Accessibility Tools

- Le blog participatif de bioinformatique francophone depuis 2012 -

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 :

Sachant que j'avais com­bi­né ensemble 16 échan­tillons, j'ai théo­ri­que­ment séquen­cé à : 400M/16=25 mil­lion de reads par échan­tillons. Si je veux simu­ler un séquen­çage à 10M, je vais donc pré­le­ver (10/25)*100=40% des reads de mes 3 FastQ. J'ai choi­si de simu­ler des séquen­çages allant d'une pro­fon­deur théo­rique de 0.1 à 20 mil­lion de reads et je me suis fait un petit script shell quick'n dir­ty pour lan­cer les simu­la­tions en paral­lèle sur un ser­veur de cal­cul.

Une fois les simu­la­tions faites, il faut map­per tout ça sur le génome de réfé­rence (dans mon cas, celui de  la sou­ris) et cal­cu­ler le nombre de gènes que l'on détecte pour chaque pro­fon­deur de séquen­çage simu­lée. On uti­lise le logi­ciel d'alignement de son choix (Bow­tie, BWA, Gem­Map­per, …), puis on quan­ti­fie l'expression des gènes. Pour ce faire, j'ai choi­si d'utiliser Cuf­flinks pour obte­nir des FPKM (Frag­ments Per Kilo­base Of Exon Per Mil­lion Frag­ments Map­ped). Cette mesure assez répan­due nor­ma­lise le nombre de reads map­pés sur un gène par rap­port à la lon­gueur de ce gène mais aus­si par rap­port au nombre total de reads map­pés sur le génome. Cette mesure nor­ma­li­sée per­met d'estimer l'abondance d'un ARN dans notre échan­tillon. On consi­dère qu'un FPKM de 1 cor­res­pond à la pré­sence d'une molé­cule d'ARN dans notre échan­tillon. On comp­ta­bi­lise le nombre de gènes détec­tés avec au moins 1 FPKM puis on met le tout sous forme de gra­phique pour obte­nir ce que l'on appelle des courbes de satu­ra­tion (mer­ci R et ggplot2 !). On part du prin­cipe qu'il existe une pro­fon­deur de séquen­çage à par­tir de laquelle on détecte tous les ARN, et donc tous les gènes expri­més dans notre échan­tillons de départ. Si on séquence au delà de cette pro­fon­deur, le nombre de gènes res­te­ra inchan­gé puisqu'on a atteint le maxi­mum détec­table. La courbe de satu­ra­tion per­met de défi­nir à par­tir de quelle pro­fon­deur on atteint le pla­teau.

Esti­ma­tion de la pro­fon­deur opti­mum de séquen­çage pour une ana­lyse d'expression de gène par RNA-seq (Isa­belle Sté­vant, CC BY)

Sur ce gra­phique, il est clair que j’atteins un pla­teau autour de 5 mil­lion de reads, et ce, pour les 3 cel­lules. Si je séquence plus pro­fond, je ne détec­te­rai pas plus de gènes. Dans mon cas, si je veux séquen­cer mes 640 cel­lules sur les 16 lignes dis­po­nibles du séquen­ceur, je devrai au maxi­mum com­bi­ner mes échan­tillons par 640 cel­lules /​16 lignes = 40, soit une pro­fon­deur de 10 mil­lion de reads par cel­lules (400M reads par ligne /​40 cel­lules). Selon la courbe de satu­ra­tion, je devrais détec­ter le même nombre de gènes qu'avec un séquen­çage plus pro­fond. J'ai vou­lu faire une der­nière véri­fi­ca­tion en fai­sant une régres­sion linéaire afin de m'assurer que les valeurs de FPKM pour chaque gène sont les mêmes entre la  pro­fon­deur à 20 mil­lion de reads  et celle à 10 mil­lion 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)
Consis­tance entre un séquen­çage à 20M reads vs un séquen­çage à 10M reads (don­nées simu­lées, Isa­belle Sté­vant, CC BY)

La cor­ré­la­tion étant très très proche de 1, ce qui veut dire que les don­nées sont extrê­me­ment simi­laires, je peux prendre la déci­sion de séquen­cer à une pro­fon­deur de 10 mil­lion 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 échan­tillons pour un seul run de séquen­çage !

Conclu­sion

Si tout comme moi vous vou­lez étu­dier l'expression des gènes mais pas for­cé­ment aller plus en détail (i.e. variants et trans­crits alter­na­tifs), sachez qu'il est pré­fé­rable de mul­ti­plier les répli­cats plu­tôt que de tout miser sur un séquen­çage pro­fond. La pro­fon­deur de séquen­çage peut être éva­luée en amont si vous déte­nez déjà des don­nées pré­li­mi­naires ou suf­fi­sam­ment simi­laires (même tis­su, même orga­nisme). Cette étape vous per­met­tra sans doute de mul­ti­plier vos expé­riences sans faire explo­ser le bud­get séquen­çage et ain­si d'améliorer vos résul­tats grâce à des ana­lyses d'expression dif­fé­ren­tielles plus robustes.

[NDLR] Le séquen­çage de mes échan­tillons est tou­jours en cours au moment de la rédac­tion de cet article, je n'ai donc pas encore la cer­ti­tude que ça ait fonc­tion­né…

[Mise à jour — 25/​08/​2015] Le séquen­çage est ter­mi­né depuis un moment et je suis en pleine ana­lyse des don­nées 🙂
Regression lineaire RNAseq

J'ai séquen­cé 8 échan­tillons à 25M reads et à 10M reads. Les librai­ries séquen­cées sont issues de deux pré­pa­ra­tions dif­fé­rentes (même cDNA à la base mais 2 librai­ries faites sépa­ré­ment par des gens dif­fé­rents, avec des kit issus de lot dif­fé­rents). Cela veut dire que les varia­tions obte­nues peuvent être cau­sées soit par la pré­pa­ra­tion des librai­ries, soit par la pro­fon­deur du séquen­çage, voire des deux.

Voi­ci le résul­tat de la com­pa­rai­son pour un des échan­tillons :

J'obtiens pour les 8 échan­tillons un pro­fil simi­laire à celui ci-contre. J'obtiens une cor­ré­la­tion supé­rieur à 0.99 entre les deux pro­fon­deurs de séquen­çage, et la droite de régres­sion est qua­si par­faite (y=x).

Conclu­sion : je n'ai pas per­du d'information en séquen­çant moins pro­fon­dé­ment, et j'ai pu séquen­cer plus du double d'échantillons par rap­port à ce que nous avions pré­vu.

Mer­ci à Mathu­rin, m4rsu, Nol­wenn, Syl­vain P.  et waque­teu pour la relec­ture et les conseils.

Vous avez aimé ? Dites-le nous !

Moyenne : 0 /​ 5. Nb de votes : 0

Pas encore de vote pour cet article.

Partagez cet article



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

Pour insérer du code dans vos commentaires, utilisez les balises <code> et <\code>.