Astuce :
RNA-seq : plus de profondeur ou plus d'échantillons ?

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 :

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.

Estimation de la profondeur optimum de séquençage pour une analyse d'expression de gène par RNA-seq (Isabelle Stévant, CC BY)

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 :

Consistence entre un séquençage à 20M reads vs un séquençage à 10M reads (données simulées, Isabelle Stévant, CC BY)

Consistance entre un séquençage à 20M reads vs un séquençage à 10M reads (données simulées, Isabelle Stévant, CC BY)

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.Regression lineaire RNAseq

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.

  • À propos de
  • Passionnée d'informatique, de logiciels libres, de graphisme, touche à tout, curieuse mais grognon.Licence de biologie et master de bioinfo à l'université de Rennes, et PhD en bioinformatique à l'université de Genève.

4 commentaires sur “RNA-seq : plus de profondeur ou plus d'échantillons ?

  1. Très intéressant ! Ça m'a rappelé l'article paru dans Bioinformatics, où les auteurs posent la même question, et arrivent à une conclusion similaire : http://bioinformatics.oxfordjournals.org/content/early/2013/12/06/bioinformatics.btt688.full.pdf

    • Je ne connaissais pas cet article, ça me rassure sur ma méthode 😀

  2. Intéressant ! Je rencontre la même problématique. Un update post-séquençage ?

  3. Merci pour la piqûre de rappel 😉
    Post séquençage, j'ai comparé 8 échantillons que j'ai séquencé à 25M reads et à 10M reads.
    Pour chacun des échantillons, j'ai une corrélation supérieure à 0.99 et une régression linéaire proche de y=x (genre y=0.0018*0.99x) quand je compare les log2(RPKM) des gènes dans les deux conditions.
    Conclusion, tout à ultra bien fonctionné \o/

    Je rajoute les graphiques à l'article dès que j'ai remis la main dessus 😉

Laisser un commentaire