Lorsque l'on se lance dans l'aventure du séquençage haut débit de transcriptome, on est amené à se poser LA question, oui LA, celle que l'on redoute à peu près tous quand on a un budget serré :
À quelle profondeur dois-je séquencer mes échantillons ?
Toutes les publications s'accordent à le dire, plus on a de réplicats, plus on a de puissance statistique pour détecter les gènes différentiellement exprimés. Or, même si les coûts de séquençage ont raisonnablement baissé, il n'est toutefois pas offert à tout le monde de séquencer des dizaines (voire des centaines) d'échantillons. Si votre projet se focalise sur l'expression des gènes, la solution pour trouver un bon compromis entre puissance statistique et coût est d'évaluer la profondeur de séquençage optimum afin de séquencer le plus d'échantillons sans pour autant perdre de l'information.
Petit rappel sur ce qu'est la profondeur de séquençage
Un séquenceur, type Illumina HiSeq 2000 par exemple, n'est capable de séquencer qu'un nombre fini de reads. Dans le cas du HiSeq 2000, le séquençage s'effectue sur 2 flow-cells contenant chacune 8 lignes (cf. cette vidéo pour en savoir plus sur le séquençage Illumina). Chaque ligne peut produire jusqu'à 200 million de reads dits single-end (SE, une seule extrémité du cDNA est lue) ou 400 million de reads dits paired-end (PE, les deux extrémités du cDNA sont lues). Il est possible de combiner plusieurs librairies (ou échantillons) sur une même ligne en utilisant des systèmes d'index qui permettent de taguer chaque fragment cDNA afin de retrouver à quel échantillon ils appartiennent. L'équation est alors simple. Plus on combine d'échantillons sur une ligne, moins on a de reads par échantillon.
Si l'on combine trop d'échantillons, on risque de passer à côté des gènes faiblement exprimés (ARN de faible abondance), mais le coût du séquençage par échantillon est faible. À l'inverse, si l'on combine peu d'échantillons, on va peut-être séquencer plus de reads qu'il n'en faut pour détecter l'expression des gènes et le coût par échantillon est beaucoup plus élevé.
Si votre projet se base sur l'expression globale des gènes et non pas sur la détection de transcrits alternatifs, il n'est peut-être pas nécessaire de séquencer 25–30 million de reads par échantillon comme il est habituel de faire. Il est possible de réduire le nombre de reads par échantillon et donc de mettre plus d'échantillons dans le séquenceur pour le même budget.
Évaluer la profondeur de séquençage
Je vais prendre un cas concret pour exposer la problématique. Dans le cadre de mon projet de thèse, je fais du single-cell RNA-seq sur des tissus embryonnaires de souris pour identifier des lignées cellulaires et suivre leurs différenciations au cours du temps. En tout, j'ai 7 conditions et un total de 640 librairies en stock (~90 cellules par condition). Le laboratoire ne peut financer qu'un seul run de séquençage (2 flow cells). Problème : dois-je séquencer seulement une partie de mes échantillons ou bien la totale ?
Au tout début du projet, nous avions réalisé un test où nous avions séquencé en paired-end une trentaine de cellules avec une profondeur de 25 millions de reads. Grâce à ce travail préliminaire, nous avons prouvé que nous étions capable de détecter tout plein de gènes nous permettant d'identifier nos différentes lignées cellulaires mais la puissance statistique restait malgré tout assez faible. En effet, une des lignées étant très faiblement représentée, les gènes spécifiques à cette population ne sortaient pas comme significatifs. Nous savons donc qu'il nous faut plus de cellules par stade pour être plus robustes dans nos analyses.
Afin de calculer combien d'échantillons je peux réussir à caser dans le séquenceur sans trop perdre d'information, j'ai été amenée à évaluer la profondeur optimale de séquençage.
Pour cela, j'ai choisi 3 fichiers FastQ de mes premiers séquençages (fichiers bruts contenant les reads pas encore mappés) représentant 3 types cellulaires différents et j'y ai prélevé de manière aléatoire un certain pourcentage de reads pour simuler des séquençages de moins en moins profond.
Pour faire ce prélèvement aléatoire de reads j'ai utilisé le logiciel HTSeq (python) qui propose plein d'outils permettant de manipuler des fichier FastQ et Bam entre autre. Vu que j'ai mis un peu de temps à trouver la méthode, je partage le code avec vous :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 |
import sys, random import HTSeq """ Prélève aléatoirement une fraction de reads dans des fichiers FastQ Exemple pour des RNA-seq paired-end reads Usage : python getRandomReads.py 0.5 input_1.fastq.gz input_2.fastq.gz output_1.fastq.gz output_2.fastq.gz """ # Pourcentage de read à prélever fraction = float( sys.argv[1] ) # Lire fichier FastQ read-pair 1 in1 = iter( HTSeq.FastqReader( sys.argv[2] ) ) # Lire fichier FastQ read-pair 2 in2 = iter( HTSeq.FastqReader( sys.argv[3] ) ) # Fichier de sortie FastQ read-pair 1 out1 = open( sys.argv[4], "w" ) # Fichier de sortie FastQ read-pair 2 out2 = open( sys.argv[5], "w" ) # La magie se passe ici while True : read1 = next( in1 ) read2 = next( in2 ) if random.random() < fraction : read1.write_to_fastq_file( out1 ) read2.write_to_fastq_file( out2 ) # Fermer les fichier créés out1.close() out2.close() |
Sachant que j'avais combiné ensemble 16 échantillons, j'ai théoriquement séquencé à : 400M/16=25 million de reads par échantillons. Si je veux simuler un séquençage à 10M, je vais donc prélever (10/25)*100=40% des reads de mes 3 FastQ. J'ai choisi de simuler des séquençages allant d'une profondeur théorique de 0.1 à 20 million de reads et je me suis fait un petit script shell quick'n dirty pour lancer les simulations en parallèle sur un serveur de calcul.
Une fois les simulations faites, il faut mapper tout ça sur le génome de référence (dans mon cas, celui de la souris) et calculer le nombre de gènes que l'on détecte pour chaque profondeur de séquençage simulée. On utilise le logiciel d'alignement de son choix (Bowtie, BWA, GemMapper, …), puis on quantifie l'expression des gènes. Pour ce faire, j'ai choisi d'utiliser Cufflinks pour obtenir des FPKM (Fragments Per Kilobase Of Exon Per Million Fragments Mapped). Cette mesure assez répandue normalise le nombre de reads mappés sur un gène par rapport à la longueur de ce gène mais aussi par rapport au nombre total de reads mappés sur le génome. Cette mesure normalisée permet d'estimer l'abondance d'un ARN dans notre échantillon. On considère qu'un FPKM de 1 correspond à la présence d'une molécule d'ARN dans notre échantillon. On comptabilise le nombre de gènes détectés avec au moins 1 FPKM puis on met le tout sous forme de graphique pour obtenir ce que l'on appelle des courbes de saturation (merci R et ggplot2 !). On part du principe qu'il existe une profondeur de séquençage à partir de laquelle on détecte tous les ARN, et donc tous les gènes exprimés dans notre échantillons de départ. Si on séquence au delà de cette profondeur, le nombre de gènes restera inchangé puisqu'on a atteint le maximum détectable. La courbe de saturation permet de définir à partir de quelle profondeur on atteint le plateau.
Sur ce graphique, il est clair que j’atteins un plateau autour de 5 million de reads, et ce, pour les 3 cellules. Si je séquence plus profond, je ne détecterai pas plus de gènes. Dans mon cas, si je veux séquencer mes 640 cellules sur les 16 lignes disponibles du séquenceur, je devrai au maximum combiner mes échantillons par 640 cellules /16 lignes = 40, soit une profondeur de 10 million de reads par cellules (400M reads par ligne /40 cellules). Selon la courbe de saturation, je devrais détecter le même nombre de gènes qu'avec un séquençage plus profond. J'ai voulu faire une dernière vérification en faisant une régression linéaire afin de m'assurer que les valeurs de FPKM pour chaque gène sont les mêmes entre la profondeur à 20 million de reads et celle à 10 million de reads :
La corrélation étant très très proche de 1, ce qui veut dire que les données sont extrêmement similaires, je peux prendre la décision de séquencer à une profondeur de 10 million de reads au lieu de 25 sans risque de perdre de l'information sur l'expression des gènes. Je passe donc de 256 à 640 échantillons pour un seul run de séquençage !
Conclusion
Si tout comme moi vous voulez étudier l'expression des gènes mais pas forcément aller plus en détail (i.e. variants et transcrits alternatifs), sachez qu'il est préférable de multiplier les réplicats plutôt que de tout miser sur un séquençage profond. La profondeur de séquençage peut être évaluée en amont si vous détenez déjà des données préliminaires ou suffisamment similaires (même tissu, même organisme). Cette étape vous permettra sans doute de multiplier vos expériences sans faire exploser le budget séquençage et ainsi d'améliorer vos résultats grâce à des analyses d'expression différentielles plus robustes.
[NDLR] Le séquençage de mes échantillons est toujours en cours au moment de la rédaction de cet article, je n'ai donc pas encore la certitude que ça ait fonctionné…
[Mise à jour — 25/08/2015] Le séquençage est terminé depuis un moment et je suis en pleine analyse des données 🙂
J'ai séquencé 8 échantillons à 25M reads et à 10M reads. Les librairies séquencées sont issues de deux préparations différentes (même cDNA à la base mais 2 librairies faites séparément par des gens différents, avec des kit issus de lot différents). Cela veut dire que les variations obtenues peuvent être causées soit par la préparation des librairies, soit par la profondeur du séquençage, voire des deux.
Voici le résultat de la comparaison pour un des échantillons :
J'obtiens pour les 8 échantillons un profil similaire à celui ci-contre. J'obtiens une corrélation supérieur à 0.99 entre les deux profondeurs de séquençage, et la droite de régression est quasi parfaite (y=x).
Conclusion : je n'ai pas perdu d'information en séquençant moins profondément, et j'ai pu séquencer plus du double d'échantillons par rapport à ce que nous avions prévu.
Merci à Mathurin, m4rsu, Nolwenn, Sylvain P. et waqueteu pour la relecture et les conseils.
Laisser un commentaire