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

Question_mark_(3534516458)
(Mar­co Bel­luc­ci, CC BY)

Lorsque l'on se lance dans l'aventure du séquen­çage haut débit de trans­crip­tome, on est ame­né à se poser LA ques­tion, oui LA, celle que l'on redoute à peu près tous quand on a un bud­get ser­ré :

À quelle pro­fon­deur dois-je séquen­cer mes échan­tillons ?

Toutes les publi­ca­tions s'accordent à le dire, plus on a de répli­cats, plus on a de puis­sance sta­tis­tique pour détec­ter les gènes dif­fé­ren­tiel­le­ment expri­més. Or, même si les coûts de séquen­çage ont rai­son­na­ble­ment bais­sé, il n'est tou­te­fois pas offert à tout le monde de séquen­cer des dizaines (voire des cen­taines) d'échantillons. Si votre pro­jet se foca­lise sur l'expression des gènes, la solu­tion pour trou­ver un bon com­pro­mis entre puis­sance sta­tis­tique et coût est d'évaluer la pro­fon­deur de séquen­çage opti­mum afin de séquen­cer le plus d'échantillons sans pour autant perdre de l'information.

Petit rap­pel sur ce qu'est la pro­fon­deur de séquen­çage

Un séquen­ceur, type Illu­mi­na HiSeq 2000 par exemple, n'est capable de séquen­cer qu'un nombre fini de reads. Dans le cas du HiSeq 2000, le séquen­çage s'effectue sur 2 flow-cells conte­nant cha­cune 8 lignes (cf. cette vidéo pour en savoir plus sur le séquen­çage Illu­mi­na). Chaque ligne peut pro­duire jusqu'à 200 mil­lion de reads dits single-end (SE, une seule extré­mi­té du cDNA est lue) ou 400 mil­lion de reads dits pai­red-end (PE, les deux extré­mi­tés du cDNA sont lues). Il est pos­sible de com­bi­ner plu­sieurs librai­ries (ou échan­tillons) sur une même ligne en uti­li­sant des sys­tèmes d'index qui per­mettent de taguer chaque frag­ment cDNA afin de retrou­ver à quel échan­tillon ils appar­tiennent. L'équation est alors simple. Plus on com­bine d'échantillons sur une ligne, moins on a de reads par échan­tillon.

Si l'on com­bine trop d'échantillons, on risque de pas­ser à côté des gènes fai­ble­ment expri­més (ARN de faible abon­dance), mais le coût du séquen­çage par échan­tillon est faible. À l'inverse, si l'on com­bine peu d'échantillons, on va peut-être séquen­cer plus de reads qu'il n'en faut pour détec­ter l'expression des gènes et le coût par échan­tillon est beau­coup plus éle­vé.

Si votre pro­jet se base sur l'expression glo­bale des gènes et non pas sur la détec­tion de trans­crits alter­na­tifs, il n'est peut-être pas néces­saire de séquen­cer 25–30 mil­lion de reads par échan­tillon comme il est habi­tuel de faire. Il est pos­sible de réduire le nombre de reads par échan­tillon et donc de mettre plus d'échantillons dans le séquen­ceur pour le même bud­get.

Éva­luer la pro­fon­deur de séquen­çage

Je vais prendre un cas concret pour expo­ser la pro­blé­ma­tique. Dans le cadre de mon pro­jet de thèse, je fais du single-cell RNA-seq sur des tis­sus embryon­naires de sou­ris pour iden­ti­fier des lignées cel­lu­laires et suivre leurs dif­fé­ren­cia­tions au cours du temps. En tout, j'ai 7 condi­tions et un total de 640 librai­ries en stock (~90 cel­lules par condi­tion). Le labo­ra­toire ne peut finan­cer qu'un seul run de séquen­çage (2 flow cells). Pro­blème : dois-je séquen­cer seule­ment une par­tie de mes échan­tillons ou bien la totale ?

Au tout début du pro­jet, nous avions réa­li­sé un test où nous avions séquen­cé en pai­red-end une tren­taine de cel­lules avec une pro­fon­deur de 25 mil­lions de reads. Grâce à ce tra­vail pré­li­mi­naire, nous avons prou­vé que nous étions capable de détec­ter tout plein de gènes nous per­met­tant d'identifier nos dif­fé­rentes lignées cel­lu­laires mais la puis­sance sta­tis­tique res­tait mal­gré tout assez faible. En effet, une des lignées étant très fai­ble­ment repré­sen­tée, les gènes spé­ci­fiques à cette popu­la­tion ne sor­taient pas comme signi­fi­ca­tifs. Nous savons donc qu'il nous faut plus de cel­lules par stade pour être plus robustes dans nos ana­lyses.

Afin de cal­cu­ler com­bien d'échantillons je peux réus­sir à caser dans le séquen­ceur sans trop perdre d'information, j'ai été ame­née à éva­luer la pro­fon­deur opti­male de séquen­çage.

Pour cela, j'ai choi­si 3 fichiers FastQ de mes pre­miers séquen­çages (fichiers bruts conte­nant les reads pas encore map­pés) repré­sen­tant 3 types cel­lu­laires dif­fé­rents et j'y ai pré­le­vé de manière aléa­toire un cer­tain pour­cen­tage de reads pour simu­ler des séquen­çages de moins en moins pro­fond.

Pour faire ce pré­lè­ve­ment aléa­toire de reads j'ai uti­li­sé le logi­ciel HTSeq (python) qui pro­pose plein d'outils per­met­tant de mani­pu­ler des fichier FastQ et Bam entre autre. Vu que j'ai mis un peu de temps à trou­ver la méthode, je par­tage le code avec vous :

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.

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.



Pour continuer la lecture :


Commentaires

4 réponses à “RNA-seq : plus de profondeur ou plus d'échantillons ?”

    1. Je ne connais­sais pas cet article, ça me ras­sure sur ma méthode 😀

  1. Avatar de mlenda

    Inté­res­sant ! Je ren­contre la même pro­blé­ma­tique. Un update post-séquen­çage ?

  2. Mer­ci pour la piqûre de rap­pel 😉
    Post séquen­çage, j'ai com­pa­ré 8 échan­tillons que j'ai séquen­cé à 25M reads et à 10M reads.
    Pour cha­cun des échan­tillons, j'ai une cor­ré­la­tion supé­rieure à 0.99 et une régres­sion linéaire proche de y=x (genre y=0.0018*0.99x) quand je com­pare les log2(RPKM) des gènes dans les deux condi­tions.
    Conclu­sion, tout à ultra bien fonc­tion­né \o/​

    Je rajoute les gra­phiques à l'article dès que j'ai remis la main des­sus 😉

Laisser un commentaire