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.

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 :

0. Intro­duc­tion
   • C'est quoi déjà, un RNA-seq ?
   • Et ces don­nées, elles res­semblent à quoi ?
1. La pré­pa­ra­tion
2. Le Contrôle de Qua­li­té ("QC")
3. L'alignement ("map­ping")
   • 3.1. En par­tant d'une anno­ta­tion exis­tante
   • 3.2. En recons­trui­sant l'annotation
   • 3.3. Le résul­tat
   • 3.4. Cuf­flinks & autres
4. Le comp­tage
   • 4.1. La pré­pa­ra­tion d'une table
   • 4.2. Le pro­blème des jonc­tions
   • 4.3. Le pro­blème des trans­crits alter­na­tifs
   • 4.4. Les dupli­cats PCR
5. La visua­li­sa­tion
   • 5.1. Genome brow­sers
   • 5.2. MA-plot
6. Le cal­cul d'expression dif­fé­ren­tielle
   • 6.1. La nor­ma­li­sa­tion
   • 6.2. Les stats
7. L'interprétation
8. 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")

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

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 :

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é , où est le vec­teur des comptes par trans­crit, celui des comptes par exon, et 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.

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 est le nombre de reads qui mappent sur un gène dans une condi­tion, et le nombre de reads du même gène dans l'autre, sur l'axe hori­zon­tal on aura , et sur l'axe ver­ti­cal . 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 : en x , et en y.

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 . 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

[maths on]

Chaque read a un chance (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 . Qu'aucun ne pro­vienne de G, , etc. Pour reads, la pro­ba­bi­li­té que d'entre eux viennent de G et pas les autres est pro­por­tion­nelle à . 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 : , 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 (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 à par­tir du nombre de reads. Si on a des répli­cats, c'est bien : on peut esti­mer 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 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.