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 :

GeneID counts1 counts2 Start End Gene­Name Strand Chro­mo­some
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

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