- Le blog participatif de bioinformatique francophone depuis 2012 -

L'analyse de données RNA-seq : mode d'emploi

Un jour, un bio­lo­giste se pointe chez vous avec d'une part un disque dur externe dans la main, d'autre part l'air sou­cieux. Il veut que vous ana­ly­siez ses don­nées RNA-seq. Le disque, c'est parce qu'il a envi­ron 50Gb de don­nées à vous trans­mettre ; l'air sou­cieux, c'est parce qu'elles ont coû­té dans les 15'000 euros, et qu'il espère que pour une fois il a pas trop foi­ré ses manips. Il compte sur vous, vu qu'il n'a aucune idée de ce qu'il a entre les mains, mais il veut des p‑valeurs en sor­tie, et des petites, s'il vous plaît.

Function finder | Duncan Hull
Func­tion fin­der | Dun­can Hull

Son but peut être mul­tiple :

  • démon­trer que cer­tains gènes sont trans­crits ;
  • démon­trer que cer­tains gènes sont trans­crits dans une condi­tion, mais pas dans l'autre ;
  • démon­trer un spli­cing alter­na­tif (ou "épis­sage") de ces gènes ;
  • etc.

Le second est le plus cou­rant. Les "condi­tions" peuvent repré­sen­ter la présence/​absence d'un fac­teur de trans­crip­tion, la pré­sence ou non d'un phénotype/​maladie, des envi­ron­ne­ments de culture dif­fé­rents, des tis­sus cel­lu­laires dif­fé­rents, etc. Le but est alors de mon­trer par exemple qu'un fac­teur de trans­crip­tion favo­rise une cer­taine mala­die lorsqu'il est expri­mé dans le foie.

Il est aus­si cen­sé savoir vous dire sur quelle espèce il tra­vaille, si c'est du "single-end" ou du "pai­red-end" (voir plus bas), et ce à quoi il s'attend dans les résul­tats.

Table des matières :

  1. Intro­duc­tion
  2. La pré­pa­ra­tion
  3. Le Contrôle de Qua­li­té ("QC")
  4. L'alignement ("map­ping")
  5. Le comp­tage
  6. La visua­li­sa­tion
  7. Le cal­cul d'expression dif­fé­ren­tielle
  8. L'interprétation
  9. Conclu­sion

C'est quoi déjà, un RNA-seq ?

Rap­pe­lons la méthode expé­ri­men­tale [1]: on récolte les ARN de nom­breuses cel­lules à la fois, on les coupe en frag­ments qu'on ampli­fie par PCR et qu'on donne à man­ger à un séquen­ceur. Celui-ci nous retourne de petits bouts de séquence géno­mique, des "reads", sous forme de texte avec des A,T,G,C. Vous héri­tez de dizaines de mil­lions de ces petites phrases.

On part du prin­cipe que le nombre de reads est pro­por­tion­nel à l'abondance des ARN cor­res­pon­dants dans la cel­lule, le but étant d'estimer cette abon­dance.

Le RNA-seq est une tech­nique rela­ti­ve­ment récente et faire par­tie de ce qu'on appelle "séquen­çage de seconde géné­ra­tion" (next-gene­ra­tion sequen­cing) ou "séquen­çage à haut débit" (high-through­put sequen­cing), avec le ChIP-seq et ses contem­po­rains.

Remarque : Les ARNm que l'on étu­die avec le RNA-seq sont matures — c'est-à-dire poly-adé­ny­lés et déjà spli­cés. C'est un point impor­tant, car il se peut qu'un gène pro­duise une grande quan­ti­té d'ARN, mais que ceux-ci soient rapi­de­ment dégra­dés avant leur matu­ra­tion. Or, cer­tains de ces ARN peuvent tout de même être réac­tifs, par exemple comme fac­teurs de trans­crip­tion d'autres gènes. Les ARN imma­tures, ou "nais­sants", peuvent être étu­diés à l'aide d'autres tech­no­lo­gies, comme le NET-seq [2].

Et ces données, elles ressemblent à quoi ?

Les fichiers qu'il vous trans­met sont le plus sou­vent au for­mat FASTQ (texte). Ce pour­rait être des SRA (binaire), mais ce for­mat appa­raît plu­tôt en fin d'analyse, au moment de la mise en ligne publique des don­nées. Un FASTQ, ça res­semble à ça :

Chaque bloc de 4 lignes, com­men­çant par un "@", repré­sente un "read": la séquence d'un frag­ment d'ARN de lon­gueur fixe, la qua­li­té (fia­bi­li­té) du séquen­çage de ce mor­ceau, et diverses autres infos illi­sibles. Vous remar­quez par la même occa­sion que le séquen­ceur pro­duit des reads courts ("short read sequen­cing"), d'environ 30 à 150nt.

Remarque : Vu la taille de ces fichiers (des cen­taines de mil­lions de lignes), un échan­tillon est sou­vent divi­sé en plu­sieurs fichiers de taille égale.

Ca peut être du "single end", où chaque read est indé­pen­dant des autres, ou du "pai­red-end", où les reads vont par paires, un sur chaque brin, à une dis­tance moyenne fixe l'un de l'autre. Le pai­red-end est néces­saire pour bien détec­ter le spli­cing, et dans un tout autre contexte, les indels. Pour le reste, le single-end fait à peu près aus­si bien pour moins cher. Cer­taines tech­niques single-end per­mettent aus­si de conser­ver l'information du brin d'origine (réfé­rence ou com­plé­men­taire), et donc le sens de trans­crip­tion.

1. La préparation

  • Copier les don­nées à un empla­ce­ment qui a beau­coup d'espace dis­po­nible et le plus pos­sible de mémoire vive — l'accès à un ser­veur de sto­ckage et à un clus­ter de cal­cul est idéal.
  • Se pro­cu­rer la séquence géno­mique de l'espèce concer­née en for­mat fas­ta (.fa), par exemple depuis NCBI. Elle sert à dif­fé­rent pro­grammes pour construire des "index" qui servent au map­ping (voir ci-des­sous).

2. Le Contrôle de Qualité ("QC")

SwissMade
Swiss­Made

Il faut s'assurer tout d'abord que le séquen­ceur a bien fait son tra­vail. Les pre­mières paires de bases d'un read sont séquen­cées avec beau­coup de fia­bi­li­té, mais plus on avance dans la séquence, moins c'est pré­cis. Un pro­gramme comme Fast­QC pro­dui­ra divers graphes per­met­tant de savoir s'il faut rétré­cir les reads, et de quelle lon­gueur — le "trim­ming".

Cette étape est très impor­tante puisque des reads de mau­vaise qua­li­té s'aligneront dif­fi­ci­le­ment sur le génome, et ren­dront toute l'analyse inutile : on a déjà vu pas­ser un échan­tillon de 2% de reads map­pés à plus de 60% après trim­ming — et le doc­to­rant pas­ser des larmes au sou­la­ge­ment.

3. L'alignement ("mapping")

A prio­ri, on ne sait pas de quelle par­tie du génome pro­vient un read, mais uni­que­ment sa séquence. On va donc "ali­gner" chaque read sur le génome de l'espèce, c'est-à-dire recher­cher dans le génome la posi­tion d'une sous-séquence simi­laire à celle du read. Celle-ci devrait en théo­rie être unique pour un read suf­fi­sam­ment long (>30nt), et donc on obtient la posi­tion d'origine du trans­crit d'où il vient.

On ima­gine le génome comme une grande ligne droite, longue de quelques mil­liards de carac­tères, et les reads comme de petits seg­ments venant s'aligner le long de cette ligne, là où la séquence est simi­laire. Si un trans­crit a été séquen­cé entiè­re­ment, tous ses exons seront "recou­verts" par des reads, et même plu­sieurs fois.

Le nombre moyen de reads par posi­tion map­pée s'appelle "pro­fon­deur de séquen­çage" ("sequen­cing depth"). C'est pour garan­tir une pro­fon­deur suf­fi­sante qu'on doit séquen­cer des cen­taines de mil­lions de reads, et s'assurer de leur qua­li­té, car c'est ce qui donne la puis­sance sta­tis­tique au moment de déter­mi­ner le taux d'expression d'un trans­crit.

L'alignement n'est jamais par­fait, à cause des erreurs de séquen­çage, de la phase d'amplification, des régions répé­tées, des muta­tions, etc. On trou­ve­ra des reads dans des régions non trans­crites, mais aus­si dans des introns (par­fois un ARN n'est pas spli­cé). Le génome de réfé­rence lui-même n'est pas exact. Il est le résul­tat d'un pro­ces­sus d'assem­blage, qui fait éga­le­ment appel au séquen­çage, et qui est un consen­sus entre une dizaine d'individus par­ti­cu­liers.

3.1. En partant d'une annotation existante

Sup­po­sons qu'on a déjà des fichiers (GTF, BED, etc.) nous indi­quant les coor­don­nées de chaque exon/​gène/​transcrit. On peut trou­ver ce genre d'annotation sur Ensem­bl (Bio­Mart) ou NCBI.

Plu­sieurs pro­grammes per­mettent d'aligner des (short) reads : Bow­tie, BWA, mais aus­si une ribam­belle d'autres (sou­vent com­mer­ciaux). Ils se valent, Bow­tie et BWA étant les algo­rithmes les plus rapides et les plus répan­dus. Nous n'avons pas tes­té per­son­nel­le­ment les pro­grammes payants.

Par­mi les options les plus impor­tantes, on veut contrô­ler le nombre d'alignements repor­tés (si un read map à plu­sieurs endroits), et le nombre d'erreurs dans la séquence ("mis­matches") per­mises. En effet, le séquen­ceur pos­sède un taux d'erreur d'environ 1%, donc on s'attend à au moins 1 erreur sur un read de 100nt. De plus, bien que des séquences de 100nt soient uniques en théo­rie — si la séquence du génome était aléa­toire, — il existe de nom­breuses région répé­tées ou très simi­laires.

3.2. En reconstruisant l'annotation

Si on ne fait pas confiance aux anno­ta­tions dis­po­nibles, ou qu'on cherche de nou­veaux exons/​isoformes, on peut recons­truire com­plè­te­ment ou par­tiel­le­ment la struc­ture du génome grâce à nos reads, au moment de l'alignement. Des pro­grammes comme TopHat ou Cuf­flinks savent devi­ner la limite des exons et recons­truire les iso­formes pos­sibles à par­tir des don­nées.

Si on recons­truit ain­si de novo l'entier du génome, cette approche est en revanche très lente et requiert une grande quan­ti­té de mémoire. Un bon com­pro­mis est de leur four­nir une anno­ta­tion (sous forme de fichier GTF) qu'ils se char­ge­ront de com­plé­ter.

3.3. Le résultat

Dans tous les cas, le résul­tat de cette étape est un fichier BAM (ver­sion binaire) ou SAM (ver­sion texte). Le SAM est rare­ment utile et prend vrai­ment trop de place, c'est pour­quoi on uti­lise en prin­cipe seule­ment le BAM, qu'on mani­pule avec sam­tools. Dans les deux cas, le conte­nu est iden­tique est res­semble à ça :

Ici chaque ligne repré­sente un ali­gne­ment : l'identifiant du read, le nom de la séquence de réfé­rence, la posi­tion de l'alignement sur la réfé­rence, la séquence du read, la qua­li­té de l'alignement, et divers indi­ca­teurs qua­li­ta­tifs ("flags"). Pour plus de détails, voir la spé­ci­fi­ca­tion du for­mat SAM.

Remarque : la qua­li­té de l'alignement n'a rien à voir avec la qua­li­té du read pré­cé­dem­ment décrite. Ici elle décrit com­bien on est sûr que le read pro­vient de cette posi­tion du génome.

3.4. Cufflinks & autres

Le nom de Cuf­flinks appa­raî­tra sou­vent dans cet article. En effet, il fait à peu près tout, avec l'inconvénient d'être répu­té pour

  • avoir un out­put com­pli­qué, voire décou­ra­geant, en tout cas à pre­mière vue ;
  • pro­duire des bugs qui vous tuent un job après 3 jours sur le clus­ter (d'après cer­tains col­lègues);
  • ne pas être per­son­na­li­sable puisque ce n'est pas votre code.

A vous de voir, on vous recom­mande tou­te­fois de le tes­ter dans sa der­nière ver­sion. A noter que d'autres alter­na­tives sont publiées à l'heure de cet article, et que d'autres le seront encore. On peut citer par exemple RSEM.

4. Le comptage

Comptage
Cal­cu­la­tor | Ans­si Kos­ki­nen (CC BY-SA 3.0)

Sous l'hypothèse que le nombre de reads venant d'un cer­tain gène est pro­por­tion­nel à l'abondance de son ARN dans la cel­lule, on veut comp­ter les reads venant de chaque gène, trans­crit ou exon du génome.

4.1. La préparation d'une table

Par­tant de nos fichiers BAM, une solu­tion est d'utiliser sam­tools dans un script mai­son. En ligne de com­mande, "sam­tools mpi­leup" four­nit par exemple un for­mat exploi­table ; le wrap­per Python de sam­tools, pysam, per­met plus de flexi­bi­li­té et d'efficacité si on code dans ce lan­gage. Il existe d'autres wrap­pers pour dif­fé­rents lan­gages (R, Ruby, C++, Java, Perl, Has­kell, Lisp).

En fait, si on fait cette par­tie soi-même (ce qui est jus­ti­fiable), on se dépa­touille comme on peut, pour autant qu'on pro­duise un fichier texte lisible — et déjà utile au bio­lo­giste, — qui res­semble déjà à une table, avec un hea­der expli­cite, du genre qu'on peut char­ger dans R faci­le­ment. Du type :

[table]

GeneID,counts1,counts2,Start,End,GeneName,Strand,Chromosome
ENSMUSG00000028180,134.0,135.0,157197123,157211354,Zranb2,1,chr3
ENSMUSG00000002010,426.0,526.0,71024301,71032236,Idh3g,-1,chrX
ENSMUSG00000002017,301.0,111.0,75936425,75951286,Fam98a,-1,chr17
ENSMUSG00000028184,565.0,478.0,148478549,148652280,Lphn2,-1,chr3
ENSMUSG00000002015,497.0,923.0,70931516,70961514,Bcap31,-1,chrX
ENSMUSG00000028186,12465.0,11435.0,146260040,146295269,Uox,1,chr3
ENSMUSG00000028189,69.0,64.0,146113433,146128813,Ctbs,1,chr3
ENSMUSG00000040147,792.0,891.0,16286407,16394492,Maob,-1,chrX
ENSMUSG00000053219,22.0,22.0,21893326,21903720,Raet1c,1,chr10

[/​table]

Sinon, une librai­rie Python, HTseq, per­met de faire la même chose si on arrive à com­prendre sa docu­men­ta­tion. Cuf­flinks, à nou­veau, sait le faire tout seul aus­si, que ce soit à par­tir d'une anno­ta­tion com­plète ou avec ce qu'il construit lui-même.

4.2. Le problème des jonctions

En ali­gnant sur le génome ou l'exome, les reads qui pro­viennent d'une jonc­tion de spli­cing (à che­val sur deux exons réunis après spli­cing de l'ARN) ne retrou­ve­ront plus leur ori­gine, puisque deux moi­tiés de la séquence du read sont sépa­rées par un intron dans la réfé­rence. C'est triste, et puis avec ça on perd une par­tie des don­nées. Pour y remé­dier, on peut

  • réali­gner les reads n'ayant pas pu être ali­gnés sur le génome, sur le trans­crip­tome, et addi­tion­ner ces reads aux exons cor­res­pon­dants ;
  • uti­li­ser un pro­gramme spé­cia­li­sé, comme SOAPs­plice ou TopHat, avec les reads non ali­gnés sur le génome, pour détec­ter les jonc­tions ;
  • lan­cer un pro­gramme comme TopHat ou Cuf­flinks avec tous les reads (très long).

Les deux der­niers fonc­tion­ne­ront infi­ni­ment mieux avec du pai­red-end, puisqu'un écart exa­gé­ré entre deux reads d'une paire sera la preuve qu'une par­tie d'un trans­crit a été cou­pée, comme c'est le cas lors du spli­cing.

4.3. Le problème des transcrits alternatifs

Quand on cherche à quan­ti­fier l'expression des dif­fé­rents trans­crits alter­na­tifs (iso­formes) d'un gène, on ne peut pas sim­ple­ment comp­ter com­bien de reads tombent entre les coor­don­nées géno­miques de cha­cun, qu'on peut trou­ver sur Ensem­bl. En effet plu­sieurs d'entre eux par­tagent les mêmes exons, et si un read mappe sur un exon com­mun, il est impos­sible de savoir de quel trans­crit il pro­vient. Pour plus de détails, voir cet autre article de notre blog.

On peut tou­te­fois essayer de dis­tri­buer les reads de manière opti­male entre les iso­formes, ou d'un autre point de vue, cher­cher la dis­tri­bu­tion la plus pro­bable. Plus concrè­te­ment, des solu­tions pos­sibles sont

  • Les moindres car­rés (avec contrainte de posi­ti­vi­té, "non-nega­tive least-squares") pour résoudre de manière opti­male le sys­tème linéaire sur­dé­ter­mi­né T=CE, où T est le vec­teur des comptes par trans­crit, E celui des comptes par exon, et C la matrice de cor­res­pon­dance indi­quant (avec des 1 et des 0) quel exon se retrouve dans quel trans­crit.
  • La méthode du maxi­mum de vrai­sem­blance appli­quée à la somme des pro­ba­bi­li­tés de chaque read d'appartenir à tel ou tel trans­crit. C'est la méthode implé­men­tée dans Cuf­flinks, par exemple.
  • Un algo­rithme EM, qui fait l'objet de plu­sieurs publi­ca­tions récentes (RSEM, iRe­ckon).

Ces méthodes donnent des résul­tats assez dif­fé­rents, la pre­mière, plus équi­table, appli­quant plu­tôt "le béné­fice du doute", tan­dis que la seconde tend à pri­vi­lé­gier un trans­crit par­ti­cu­lier. Nous n'avons pas encore pu tes­ter la troi­sième. On ajoute par­fois un LASSO pour "cor­ri­ger" cer­tains défauts de la méthode.

Il est très dif­fi­cile de savoir laquelle est la plus réa­liste, n'ayant pas beau­coup de cas (pas du tout?) de gènes dont on a mesu­ré par un autre moyen l'expression rela­tive des iso­formes (qui plus est pour votre type de cel­lule, dans vos condi­tions, etc.). Les résul­tats sont donc assez hypo­thé­tiques, et la solu­tion à ce genre de ques­tions vien­dra pro­ba­ble­ment d'un autre type d'expérience que le RNA-seq. On peut cepen­dant uti­li­ser la qRT-PCR pour esti­mer plus pré­ci­sé­ment l'abondance d'un petit nombre de trans­crits, ou faire des simu­la­tions. Le sujet est encore très ouvert.

4.4. Les duplicats PCR

Durant l'étape d'amplification des frag­ments d'ARN (PCR) qui pré­cède le séquen­çage, cer­tains frag­ments peuvent être sur-ampli­fiés par erreur. Le résul­tat est un nombre de copies anor­mal de reads iden­tiques, qui peut biai­ser les résul­tats. Heu­reu­se­ment, de tels "dupli­cats PCR" sont assez faci­le­ment repé­rables après visua­li­sa­tion (sec­tion sui­vante), puisqu'on observe des colonnes dis­pro­por­tion­nées de plu­sieurs cen­taines de reads iden­tiques super­po­sés. C'est beau­coup plus dif­fi­cile à détec­ter et en amont et sur l'ensemble des gènes, ce qui pose évi­dem­ment un pro­blème sta­tis­tique au moment de la com­pa­rai­son entre échan­tillons, si un frag­ment est sur-ampli­fié dans l'un et pas dans l'autre. Il n'existe par contre pas encore à notre connais­sance de bon outil pour fil­trer ces dupli­cats de manière sys­té­ma­tique.

5. La visualisation

5.1. Genome browsers

Pour avoir déjà un aper­çu de ses don­nées, le mieux est de construire pour chaque échan­tillon ce qu'on appelle un "pro­fil" ou une "den­si­té" du signal, c'est-à-dire une sorte de "bar plot" qui à chaque posi­tion du génome indique le nombre de reads map­pés à cette posi­tion. Cela se pré­sente sous la forme d'un fichier bien for­ma­té qu'un pro­gramme de visua­li­sa­tion ("genome brow­ser") comme UCSC (en ligne) ou IGV (local) sau­ra inter­pré­ter pour pro­duire un gra­phique.

UCSC view
Un exemple de vue sur UCSC, avec de haut en bas des pro­fils RNA-seq,
des CpG-islands, des élé­ments répé­tés, les gènes Ref­seq et Ensem­bl,
le GC content et des pics de ChIP-seq.

C'est extrê­me­ment utile pour plu­sieurs rai­sons :

  • On voit les "mon­tagnes" de reads qui loca­lisent les exons et autres trans­crits.
  • On peut com­pa­rer visuel­le­ment les niveaux d'expression de plu­sieurs échan­tillons, et des gènes entre eux.
  • On peut ajou­ter des pro­fils d'autres types de signaux, comme du ChIP-seq, de la méthy­la­tion, de la conser­va­tion, des motifs, etc. pour expli­quer notre signal. C'est une approche très moderne et qui est sou­vent néces­saire à jus­ti­fier une décou­verte : si un gène est expri­mé dans une condi­tion et pas dans une autre, on s'attend à obser­ver de la Poly­mé­rase II, un pic de H3K4 méthy­lé au début, un pic de machin au milieu, un pic de truc à la fin.

Pour obte­nir ces den­si­tés, il existe une mul­ti­tude de pro­grammes que chaque labo de bioin­fo a reco­dé pour lui… La fonc­tion "geno­me­Co­ve­ra­ge­Bed" de bed­tools peut appa­rem­ment le faire, sinon télé­char­gez ou codez vous-même un n‑ème "bam2wig", sachant que ça prend du BAM en entrée (à nou­veau, on peut uti­li­ser sam­tools ou son wrap­per pysam) et que ça doit res­sor­tir un bed/​bigBed/​bedGraph/​wig/​bigWig.

Pour les visua­li­seurs en ligne, sachez tout d'abord qu'il y a moyen de voir les BAM direc­te­ment, et ensuite que les for­mats pré­fé­rés sont ceux où il n'y a pas besoin de télé­char­ger tout le fichier : on upload seule­ment un "index", et quand on fait une requête pour une nou­velle région, ça vient télé­char­ger depuis votre ordi seule­ment la par­tie qu'il faut. C'est le cas des bigBed/​big­Wig, for­mats d'ailleurs inven­tés pour le brow­ser UCSC qui four­nit les outils de conver­sion (p.ex. bed­Graph­To­Big­Wig), et des BAM. Ce sera alors énor­mé­ment plus rapide.

5.2. MA-plot

Un type de graphe est aus­si fré­quem­ment uti­li­sé pour com­pa­rer deux échan­tillons, une fois les comptes obte­nus : le MA-plot. La lettre M repré­sente les log-ratios, et la lettre A est pour "Ave­rage" (la moyenne). Comme son nom l'indique, le MA-plot est un plot de la moyenne (en x) contre le ratio (en y) de deux valeurs : l'expression d'un gène dans deux condi­tions, chaque gène étant un point du graphe.

C'est-à-dire : si x_1 est le nombre de reads qui mappent sur un gène dans une condi­tion, et x_2 le nombre de reads du même gène dans l'autre, sur l'axe hori­zon­tal on aura \frac{x_1+x_2}{2}, et sur l'axe ver­ti­cal \frac{x_1}{x_2}. En fait, avec ça le graphe serait dif­fi­cile à lire (concen­tré dans un coin), c'est pour­quoi on trans­forme l'échelle en log : \frac{1}{2}(log_{10}x_1+log_{10}x_2) = log_{10}\sqrt{x_1\cdot x_2} en x , et log_2(\frac{x_1}{x_2}) en y.

MA-plot
MA-plot : le log10 de la moyenne géo­mé­trique en x ; le log2 du ratio en y. En vio­let, dif­fé­rents quan­tiles.
(Ici des qua­si répli­cats). Source : HTSs­ta­tion

Il s'interprète comme suit : plus un point est à droite du graphe, plus le gène est glo­ba­le­ment for­te­ment expri­mé, et inver­se­ment. Plus un point est en haut du graphe, plus il est sur­ex­pri­mé dans la pre­mière condi­tion par rap­port à la seconde, et inver­se­ment. Les gènes "inté­res­sants" se trouvent donc plu­tôt vers la droite (ratios plus convain­cants si les comptes son éle­vés) et le plus pos­sible en haut ou en bas du graphe, visi­ble­ment sépa­rés de la masse des autres points.

Deux répli­cats tech­niques mon­tre­ront par exemple très peu de varia­tion autour de l'axe hori­zon­tal, alors que deux échan­tillons assez dif­fé­rents don­ne­ront au graphe une forme de gros pois­son. On peut véri­fier aus­si que la médiane est bien à zéro, sans quoi on observe un biais : les échan­tillons seraient mal nor­ma­li­sés (voir plus loin).

C'est une approche pure­ment visuelle qui ne démontre rien, mais si un gène qu'on sup­pose jouer un rôle impor­tant se démarque faci­le­ment sur ce graphe, on sait déjà qu'on a réus­si, et on peut ajou­ter une très jolie figure à la publi­ca­tion, ce qui n'est pas non plus négli­geable.

La librai­rie R lim­ma sait le faire, mais on peut faci­le­ment pro­duire son propre script, tel ce package Python, qui a l'avantage d'être inter­ac­tif : si on clique sur un point, le nom du gène asso­cié s'affiche.

6. Le calcul d'expression différentielle

Là, ça se com­plique. Pour un gène don­né, on veut savoir s'il est signi­fi­ca­ti­ve­ment plus expri­mé dans une condi­tion que dans l'autre.

Par exemple, si j'ai 100 reads dans la condi­tion 1, et 500 dans la condi­tion 2, on dirait que oui (mais en fait, on n'en sait rien). Et s'il y en a 10 et 50 ? 1 et 5 ? la pro­por­tion est la même, et pour­tant ça pour­rait être dû au hasard. Et s'il y en a 100 et 150 ? Là, on a besoin des mathé­ma­tiques pour se convaincre.

6.1. La normalisation

Quand on veut com­pa­rer des échan­tillons, il faut d'abord les rendre com­pa­rables. En effet, des biais comme le nombre total de reads ali­gnés pour un échan­tillon, des arté­facts PCR, etc. peuvent faus­ser les résul­tats. Par exemple, ima­gi­nons qu'il y a 10 fois plus de reads dans une condi­tion que dans l'autre, alors la qua­si tota­li­té des gènes sera décla­rée signi­fi­ca­ti­ve­ment enri­chie.

Pour y remé­dier, on a pen­sé à divi­ser chaque échan­tillon par le nombre total de ses reads. Mal­heu­reu­se­ment, comme en RNA-seq de petites régions ont ten­dance à accu­mu­ler une grande par­tie des reads, cela crée plus de biais que ça n'en cor­rige. Ima­gi­nez deux échan­tillons où l'un a deux fois plus de reads que l'autre, mais où la moi­tié est conte­nue dans un seul gène. Alors en divi­sant par deux, on divise tous les gènes par deux, et tous deviennent signi­fi­ca­ti­ve­ment moins expri­més…

On peut aus­si consi­dé­rer de divi­ser les comptes dans chaque gène par la taille de ce gène, pour obte­nir une "den­si­té" de reads. C'est alors un autre point de vue : cette trans­for­ma­tion est néces­saire pour com­pa­rer des gènes d'un même échan­tillon entre eux, mais pas pour com­pa­rer dif­fé­rents échan­tillons. Notons cepen­dant que cer­taines publi­ca­tions [3] montrent un biais dans la détec­tion de gènes dif­fé­ren­tiel­le­ment expri­més en faveur des gènes les plus longs. Mal­gré une perte de puis­sance au niveau sta­tis­tique, il pour­rait être pré­fé­rable de divi­ser les comptes par la lon­gueur des trans­crits (?).

En alliant les deux, on obtient des RPKM, pour Reads Per Kilo­base and per Mil­lion (of reads). La for­mule est donc 10^6\cdot 10^3 \cdot \frac{\#\_reads\_de\_G}{\#\_reads\_total \quad\cdot\quad longueur\_de\_G} . Ca, les bio­lo­gistes, ils aiment bien, parce qu'on leur a dit que c'était nor­ma­li­sé. En fait, pour la rai­son évo­quée ci-des­sus, cela ne convient pas au RNA-seq.

Une meilleure solu­tion est par exemple pré­sen­tée dans l'article lié à DESeq [4] dont nous par­le­rons plus loin. En gros, on construit une nou­velle "réfé­rence" en pre­nant la moyenne (géo­mé­trique, plus robuste) des échan­tillons et en les nor­ma­li­sant tous par rap­port à cette réfé­rence.

Remarque : en fait, le RNA-seq ne per­met pas de dif­fé­ren­cier une ampli­fi­ca­tion natu­relle de l'expression de tous les gènes d'un indi­vi­du, d'un arte­fact tech­nique d'amplification. Dans le pre­mier cas il ne fau­drait pas nor­ma­li­ser, dans le second, si. On part tou­jours de l'hypothèse que le niveau d'expression de la plu­part des gènes ne change pas entre deux échan­tillons — la plu­part du temps.

6.2. Les stats

Poisson distribution
Dis­tri­bu­tion de pro­ba­bi­li­té de Pois­son
pour dif­fé­rentes valeurs de lamb­da.
(source : Wiki­pe­dia)

[maths on]

Chaque read a un chance p (très petite) de pro­ve­nir de ce gène G par­ti­cu­lier. Pour 2 reads, la pro­ba­bi­li­té qu'ils pro­viennent tous les deux de G est p*p. Qu'aucun ne pro­vienne de G, (1-p)*(1-p), etc. Pour N reads, la pro­ba­bi­li­té que k d'entre eux viennent de G et pas les autres est pro­por­tion­nelle à p^k * (1-p)^{(N-k)} . On appelle ça une loi bino­miale. Le pro­blème avec ce modèle, c'est que les expo­sants ici sont très gros : N \approx 10^8, assez pour que même le clus­ter se mette en grève.

Heu­reu­se­ment, on peut approxi­mer une bino­miale avec une loi dite de Pois­son. Une for­mule nous donne alors la pro­ba­bi­li­té d'observer un cer­tain nombre de reads (connu), si on connaît un para­mètre \lambda (incon­nu) qui n'est autre que la moyenne de cette dis­tri­bu­tion, d'après la théo­rie. On ne s'intéresse pas vrai­ment à cette pro­ba­bi­li­té, au contraire on va sup­po­ser que notre obser­va­tion est pro­bable (puisqu'elle est s'est pro­duite) et esti­mer ce para­mètre \lambda à par­tir du nombre de reads. Si on a des répli­cats, c'est bien : on peut esti­mer \lambda grâce à la moyenne arith­mé­tique des comptes des dif­fé­rents répli­cats. Si on a un seul échan­tillon de cette condi­tion, eh bien on estime la moyenne comme étant le compte de reads unique obser­vé — ce qui est sta­tis­ti­que­ment moyen, dans tous les sens du terme. Ensuite, on va com­pa­rer les \lambda dans deux condi­tions pour savoir si sta­tis­ti­que­ment les moyennes sont les mêmes.

Remarque : en réa­li­té, on observe que la loi de Pois­son est un peu trop res­tric­tive, et on uti­lise plu­tôt une loi "néga­tive bino­miale", façon d'ajouter un terme de variance sup­plé­men­taire. Un chan­ge­ment qui amène des com­pli­ca­tions puisqu'il faut esti­mer cette variance à son tour…

[maths off]

Ras­su­rez-vous, un pro­gramme le fera aima­ble­ment à votre place, donc au pire vous pou­vez joyeu­se­ment oublier tout ça et appuyer sur le bou­ton rouge.

Actuel­le­ment, les pro­grammes le mieux à même de faire ces cal­culs sont des librai­ries R du nom de DESeq et EdgeR. Tous deux pos­sèdent un manuel très clair, un article lisible et une réfé­rence com­plète. Le pro­gramme, don­née votre table obte­nue pré­cé­dem­ment, vous retourne un nou­veau fichier tabu­lé avec une belle colonne de p‑valeurs ajus­tées, dont vous pou­vez extraire une liste des gènes les plus signi­fi­ca­ti­ve­ment sur-/sous-expri­més dans une cer­taine condi­tion. Il se charge aus­si de la nor­ma­li­sa­tion si néces­saire de votre table avant de faire quoi que ce soit.

Cuf­flinks pré­tend aus­si faire tout cela très bien dans sa nou­velle ver­sion, ce que nous n'avons pas tes­té. Le modèle mathé­ma­tique est glo­ba­le­ment le même dans tous les cas.

7. L'interprétation

A pré­sent, il faut voir si les gènes qui ont été trou­vés signi­fi­ca­tifs ont quelque chose à voir avec l'expérience. Une liste de 100 gènes, en soi, n'est pas très utile. Sou­vent le bio­lo­giste sau­ra indi­quer quels gènes il s'attend à trou­ver dans la liste, et on espère qu'ils sont dans les pre­miers. Si au contraire il ne sait pas ce qu'il cherche, il faut être conscient que la démarche est explo­ra­toire, et des argu­ments Baye­siens [5] mon­tre­ront que les résul­tats n'ont en aucun cas valeur de preuve.

Enfin, si les p‑valeurs ajus­tées sont trop grandes, (c'est nul et) on peut tou­jours fil­trer par p‑valeur simple, voire par "fold change" (ratio) pour faire croire au bio­lo­giste que non, tout espoir n'est pas per­du. L'effet est par­ti­cu­liè­re­ment obser­vé en l'absence de répli­cats.

Si c'est à peut près tout ce qu'on peut tirer du RNA-seq lui-même, l'analyse peut faire inter­ve­nir d'autres types de don­nées, qui rendent le jeu plus com­plexe et pro­ba­ble­ment plus inté­res­sant. Des contraintes par­ti­cu­lières peuvent aus­si vous obli­ger à chan­ger vos para­mètres de map­ping, vos limites de p‑valeurs, etc.

8. Conclusion

Le RNA-seq est une méthode rela­ti­ve­ment récente, et par consé­quent les outils dis­po­nibles non seule­ment ne sont pas aus­si au point que pour les microar­rays, mais qu'il évo­luent rapi­de­ment ou sont régu­liè­re­ment rem­pla­cés par d'autres — en par­ti­cu­lier ceux qui pro­posent un pipe­line com­plet. Cepen­dant c'est une tech­nique pro­met­teuse, qui se pré­sente comme le suc­ces­seur des microar­rays en livrant plus d'informations struc­tu­relles tout en ayant une sen­si­bi­li­té com­pa­rable ou supé­rieure, déjà avec les outils d'analyse exis­tants. Son prin­ci­pal défaut étant son coût, le fait que celui-ci baisse très rapi­de­ment la rend de plus en plus popu­laire, et donc rend de plus en plus néces­saire la mise au point de pipe­lines effi­caces.

Et voi­là.

Réfé­rences de l'article :
[1] A. Mor­ta­za­vi, "Map­ping and quan­ti­fying mam­ma­lian trans­crip­tomes by RNA-Seq", Nature Methods 2008

[2] L. S. Church­man, "Nascent trans­cript sequen­cing visua­lizes trans­crip­tion at nucleo­tide reso­lu­tion", Nature 2011

[3] A. Osh­lack, "Trans­cript length bias in RNA-seq data confounds sys­tems bio­lo­gy", Bio­lo­gy Direct 2009

[4] S.Anders, "Dif­fe­ren­tial expres­sion ana­ly­sis for sequence count data", Genome Bio­lo­gy 2010

[5] J.P.A. Ioan­ni­dis, "Why Most Publi­shed Research Fin­dings Are False", PLOS Medi­cine 2005

Les images ne pos­sé­dant pas de sources sont de l'auteur de l'article.

Mer­ci à Clem_​, Haut­bit, nahoy, ZaZo0o et Yoann M. pour leur relec­ture.

Vous avez aimé ? Dites-le nous !

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

Pas encore de vote pour cet article.

Partagez cet article :




Commentaires

19 réponses à “L'analyse de données RNA-seq : mode d'emploi”

  1. Excellent ! Mer­ci Julien pour ce très bon article ! 🙂

    Les lec­teurs seront peut-être aus­si inté­res­sés par ce récent Nature Pro­to­col publié par Simon & co qui s'intéresse à l'aspect dif­fé­ren­tiel :

    Count-based dif­fe­ren­tial expres­sion ana­ly­sis of RNA sequen­cing data using R and Bio­con­duc­tor

    1. Mer­ci pour le lien Maga­li 🙂

  2. You're wel­come ! 😀

  3. […] la méthode dans deux billets du site Bioin­fo-fr : Ana­lyse des don­nées de séquen­çage à ARN et L’analyse de don­nées RNA-seq : mode d’emploi. Ce deuxième fait aus­si appa­raître cer­tains des sou­cis qu’on peut […]

  4. Très bon article ! Pour infor­ma­tion, le site OMIC­tools pro­pose une liste d'outils spé­ci­fiques pour l'analyse de l'ARN-seq (http://​omic​tools​.com/​m​r​n​a​-​s​eq/).

    1. Mer­ci ! C'est pas mal, ce lien. Ce serait quand même beau­coup plus pra­tique si on pou­vait trier les outils au moins par date de publi­ca­tion, ou par appré­cia­tions d'utilisateurs. Ce qui est dif­fi­cile c'est jus­te­ment de trou­ver le bon pro­gramme !

      1. Avatar de Mitra Barzine
        Mitra Barzine

        Hum, quite à me faire pas­ser pour une spam­meuse*:

        un début de réponse pour le map­page peut être trou­vée à l'adresse sui­vante : http://​www​dev​.ebi​.ac​.uk/​f​g​/​h​t​s​_​m​a​p​p​e​rs/ (site web com­pa­gnon du papier sui­vant : http://​bio​in​for​ma​tics​.oxford​jour​nals​.org/​c​o​n​t​e​n​t​/​e​a​r​l​y​/​2​0​1​2​/​1​0​/​1​1​/​b​i​o​i​n​f​o​r​m​a​t​i​c​s​.​b​t​s​6​0​5​.​a​b​s​t​r​act)

        Je sais que Nuno est en train de tra­vailler en ce moment sur les spli­ceurs, par contre je ne sais pas quand il compte sou­mettre son papier.

        Par ailleurs, il a déve­lop­pé un pipe­line iRAP (http://​code​.google​.com/​p​/​i​r​ap/) qui per­met de plug­ger faci­le­ment un outil à l'autre. Espé­rons que ça moti­ve­ra les gens à se poser un mini­mum de ques­tions avant d'opter pour la faci­li­té avec "Tuxe­do suite"(http://​www​.nature​.com/​n​p​r​o​t​/​j​o​u​r​n​a​l​/​v​7​/​n​3​/​f​u​l​l​/​n​p​r​o​t​.​2​0​1​2​.​0​1​6​.​h​tml)

        *pre­mier, second et avant-der­nier auteurs (étaient/)sont de mon labo et le der­nier auteur est de mon ins­ti­tut aus­si.

        1. Avatar de Mitra Barzine
          Mitra Barzine

          (Zut, j'ai oublier de men­tion­ner que jus­te­ment sur le pre­mier lien, *il y a* un fil chro­no­lo­gique de publi­ca­tions pour les map­peurs)

  5. Avatar de Juliette
    Juliette

    Mer­ci beau­coup pour ces expli­ca­tions ! En tant que bio­lo­giste, elles étaient très com­pré­hen­sibles !

  6. Mer­ci pour ces expli­ca­tions ; je com­mence en Jan­vier un Stage ( M2 Bioin­fo ) sur l'analyse des don­nées de RNA-seq ; cet article m'aidera beau­coup

  7. Mer­ci a vous pour vos reac­tions, ca fait plai­sir !

  8. Pour l'analyse, mais sur­tout pour en savoir plus sur ce qui se passe avant (avant que le bio­lo­giste ne débarque avec son disque dur externe), je trouve que RNA-seq­lo­pe­dia est une res­source claire et assez com­plète.

  9. Avatar de Lenjyco

    Très bon article, qui explique bien, cela m'as per­mis de com­prendre pas mal de chose sur le séquen­çage. Il fau­drait par­ler d'Eoulsan qui est un soft(français) qui com­pile un peu les dif­fé­rentes étapes et qui peut se déployer sur clus­ter.

    1. Je ne recom­man­de­rais pas du tout d'utiliser un pipe­line com­plet, mais de choi­sir ses outils. 1. Pour la repro­duc­ti­bi­li­té, il faut savoir exac­te­ment quelle ver­sion de quoi a été uti­li­sée. 2. Parce que chaque labo de bioin­fo a ses pipe­lines, mais beau­coup font n'importe quoi.
      La vraie chose à faire, c'est de lan­cer dif­fé­rents pro­grammes qui font la même chose pour com­pa­rer les résul­tats.

  10. Vrai­ment mer­ci infi­ni­ment pour ces pré­cieuses expli­ca­tions !!!! je com­prends enfin 🙂

  11. Avatar de chettab

    Bon­jour,

    Votre article est très inté­res­sant. Mer­ci pour votre aide

    K. Chet­tab

  12. Avatar de Yoann

    "Par­tant de nos fichiers BAM, une solu­tion est d'utiliser sam­tools dans un script mai­son. En ligne de com­mande, "sam­tools mpi­leup" four­nit par exemple un for­mat exploi­table;"

    Bon­jour,
    Serait-il pos­sible d'avoir plus de détails sur ces deux phrases ? Parce que le fichier que moi j'obtiens est loin d'être "exploi­table", du coup j'aimerais savoir par quelle(s) option(s) exac­te­ment de sam­tools mpi­leup, on obtient ce fichier exploi­table.

    Celle-ci étant celle que j'ai lais­sé :
    sam­tools mpi­leup possorted_genome_bam.bam ‑o s1_samtools_res.txt

  13. Avatar de Bastien Hervé
    Bastien Hervé

    Update 2019 des outils d'analyse RNA­seq :

    Le Contrôle de Qua­li­té ("QC")

    Fast­QC est tou­jours d'actualité, cou­plé à Mul­ti­QC qui per­met de faire en rap­port HTML glo­bal de vos jeux de don­nées. Fastp, sor­ti récem­ment est éga­le­ment uti­li­sé

    L'alignement ("map­ping")

    Je ne suis pas d'accord pour uti­li­ser des outils comme BWA et Bowtie2 dans une ana­lyse RNA­seq dans le cas d'un ali­gne­ment sur génome de réfé­rence, pour la bonne rai­son que ces ali­gneurs ne sont pas "splice aware" (avec leurs para­mètres par defaut), c'est à dire qu'ils ne prennent pas en compte le fait qu'un read peut ne pas se retrou­ver en entier sur le génome de réfé­rence, car l'ARN contrai­re­ment à l'ADN est spli­cé.
    Quand l'alignement est fait sur un trans­crip­tome, celui ci est lui aus­si spli­cé on peut donc uti­li­ser des ali­gneurs "non splice aware" tel que les deux cités plus haut, mais pas avec un génome.
    Deux ali­gneurs "splice aware" répan­dus : STAR et HISAT2.
    HISAT2 est sim­ple­ment l'update de TopHat, un peu lent et gour­mand en mémoire. Il faut éga­le­ment faire le trim­ming de vos reads avant.
    STAR est assez rapide est beau­coup moins deman­deur en mémoire. Il peut éga­le­ment faire l'étape de comptes en ajou­tant l'option adé­quat.
    Les deux logi­ciels offres des résul­tats com­pa­rables.

    Nou­veau­té !

    Récem­ment des logi­ciels de pseu­do-ali­gne­ments ont vu le jour, Kal­lis­to et Sal­mon pour les plus connus, il existe aus­si Sail­fish, RNA-Skim… Ces outils font de la quan­ti­fi­ca­tion de trans­crits sans ali­gner véri­ta­ble­ment les reads sur un trans­crip­tome uni­que­ment. Ils peuvent éga­le­ment, si vous le sou­hai­tez, ali­gner les reads sur le trans­crip­tome pour géné­rer des BAM si vous en avez la néces­si­té.
    Les gros avan­tages de ces outils c'est qu'ils sont très rapides et vous libère du sto­ckage des fichiers bam qui peuvent par­fois être mons­trueux. Ils font donc 2 étapes en 1, ali­gne­ment + comp­tage.
    Sal­mon est semble t‑il plus rapide que Kal­lis­to.

    Le comp­tage

    Si vous n'avez pas choi­sit la solu­tion des pseu­dos ali­gne­ments, vous pou­vez comp­ter vos reads avec fea­tu­re­Counts ou HTseq. Le pre­mier semble plus répan­dus et je n'ai pas tes­té HTseq donc je ne ferai pas d'autre com­pa­ra­tif.

    5/​6. Visualisation/​Le cal­cul d'expression dif­fé­ren­tielle

    Les packages en vogue sont DESeq2, EdgeR et lim­ma avec cha­cun leurs types de nor­ma­li­sa­tion pour contrô­ler les effets de batch, les tailles de librai­ries, les com­po­si­tions en ARN… J'ai fait un petit post sur Bios­tars au sujet de ces nor­ma­li­sa­tions (https://​www​.bios​tars​.org/​p​/​3​4​9​8​8​1​/​#​3​5​0​181)
    Le plus répan­du est DESeq2 avec de super tuto­riels (http://​bio​con​duc​tor​.org/​p​a​c​k​a​g​e​s​/​r​e​l​e​a​s​e​/​b​i​o​c​/​v​i​g​n​e​t​t​e​s​/​D​E​S​e​q​2​/​i​n​s​t​/​d​o​c​/​D​E​S​e​q​2​.​h​tml)

    Voi­là, en espé­rant avoir été assez clair 🙂 Les étapes de ce tuto sont cor­rectes et bien expli­quées par l'auteur, mais un update des outils était néces­saire.

  14. Super article très com­plet et faci­li­tant la com­pré­hen­sion de tous ces outils ! Mer­ci beau­coup 🙂

Laisser un commentaire