Manipulation d'intervalles génomiques dans R

Introduction

Nous avons abor­dé, dans le pré­cé­dent article de cette série, les bases de la mani­pu­la­tion d'intervalles dans R.

Ce deuxième article a pour objec­tif de mon­trer com­ment mani­pu­ler des inter­valles géno­miques. Au niveau le plus basique, un inter­valle est défi­ni par deux nombres entiers posi­tifs déli­mi­tant son début et sa fin. Nous allons pou­voir ajou­ter des couches sup­plé­men­taires d'informations, telles que les chro­mo­somes et leurs tailles ain­si que le brin.

Nous allons pour cela uti­li­ser une biblio­thèque de la suite Bio­con­duc­tor, Geno­mi­cRanges.

Cette biblio­thèque est en quelque sorte l'équivalent des fameux bed­tools, très uti­li­sés pour les ana­lyses en géno­mique.

Installer et charger la bibliothèque GenomicRanges

Com­men­çons par ins­tal­ler et char­ger la biblio­thèque.

Créer un objet GRanges

Coordonnées des intervalles

Nous allons créer un objet conte­nant des inter­valles, ain­si que des infor­ma­tions sur les chro­mo­somes et les brins de chro­mo­somes qui portent ces inter­valles. Un objet granges est l'équivalent d'un fichier au for­mat .bed (Brow­ser Exten­sible Data) conte­nant au mini­mum trois colonnes : le nom de la séquence, un inter­valle (un objet qui com­prend une posi­tion de début et une posi­tion de fin) et le brin (posi­tif "+" ou néga­tif "-").

Dans le code ci-des­sus, la fonc­tion c() per­met de com­bi­ner des élé­ments dans un vec­teur (ou une colonne). La fonc­tion rep(x, n) per­met de répé­ter 'n' fois l'objet x. Ain­si, avec le code seqnames = c(rep("chr1", 4), rep("chr2", 4)), nous assi­gnons les quatre pre­miers inter­valles au chro­mo­some 1 (chr1), et les quatre inter­valles sui­vants au chro­mo­some 2 (chr2).

On indique aus­si que la moi­tié de nos inter­valles sont sur le brin '+' et le reste sur le brin '-'. Si le brin n'a pas d'importance, on uti­lise alors strand = "*" lors de la créa­tion de l'objet granges. L'argument ignore.strand = TRUE per­met d'ignorer le brin dans la plu­part des fonc­tions de la biblio­thèque Geno­mi­cRanges.

Voyons à quoi res­semble notre nou­vel objet :

Métadonnées

Notre objet gr ne contient pour le moment que des infor­ma­tions néces­saire pour la créa­tion des inter­valles ("GRanges object with 8 ranges and 0 meta­da­ta columns"). Nous pou­vons ajou­ter des infor­ma­tions sup­plé­men­taires comme des noms de gènes, le taux en GC, etc. à l'aide du para­mètre mcols().

Notre objet res­semble main­te­nant à ceci :

Informations sur les chromosomes

Pour l'instant, les inter­valles de notre objet sont situés sur des brins de chro­mo­somes, mais nous n'avons pas pré­ci­sé leurs tailles, ni le génome ("2 sequences from an uns­pe­ci­fied genome ; no seq­lengths"). Sans ces infor­ma­tions, nous ne pou­vons faire des opé­ra­tions que sur les régions de chro­mo­somes cou­vertes par ces inter­valles.

La com­mande seqinfo() per­met d'afficher et modi­fier les infor­ma­tions liées aux séquences :

Les infor­ma­tions concer­nant la taille des chro­mo­somes de l'espèce ou du génome de réfé­rence peuvent être trou­vées sur le site de l'UCSC ou celui d'Ensem­bl. La biblio­thèque BSge­nomes, per­met aus­si de récu­pé­rer ces infor­ma­tions pour cer­tains génomes.

Nous pou­vons aus­si ajou­ter ces infor­ma­tions de la façon sui­vante :

Affi­chons à nou­veau les infor­ma­tions à l'aide de seqinfo() :

Extraire des données

Nous pou­vons accé­der aux infor­ma­tions à l'aide des com­mandes de base sui­vantes :

Nous pou­vons éga­le­ment uti­li­ser des opé­ra­teurs logiques :

Opérations sur des intervalles génomiques

Un seul objet granges

Nous pou­vons effec­tuer des opé­ra­tions sur notre objet :

• reduce(gr) : fusion­ner les inter­valles che­vau­chants ou adja­cents
• gaps(gr) : extraire les posi­tions non cou­vertes par des inter­valles
• coverage(gr) : cal­cu­ler la cou­ver­ture (nombre de fois que chaque posi­tion est cou­verte par un inter­valle)

D'autres fonc­tions existent :

• flank(gr, width = n, start = TRUE) : extraire une région de lon­gueur 'n' en amont. Pour extraire une région de lon­gueur 'n' en aval, uti­li­ser start = FALSE.
• shift(gr, shift = n) : dépla­cer les inter­valle de 'n' posi­tions vers la droite (ne tient pas compte du brin). Pour déca­ler les inter­valles vers la gauche, uti­li­ser shift = -n.
• resize(gr, width = n, fix = "start") : modi­fier la taille des inter­valles de lon­gueur 'n', en dépla­çant la fin de l'intervalle (fix = "start") ou le début de l'intervalle (fix = "end").
• promoters(gr, upstream = x, downstream = y) : extraire les régions pro­mo­trices, entre x posi­tions en amont (par défaut 2000) et y posi­tions en aval (par défaut 200)

Plusieurs objets granges

Après avoir vu dif­fé­rentes opé­ra­tions sur un seul objet granges, voyons com­ment com­pa­rer plu­sieurs objets.

Com­men­çons par créer deux objets : gr1 et gr2.

Les fonc­tions sui­vantes per­mettent d'extraire l'union et l'intersection de deux objets granges :

• union(gr1, gr2) : union des deux objets.
• intersect(gr1, gr2) : inter­sec­tion des deux objets.

Nous pou­vons éga­le­ment recher­cher et comp­ter les régions où les inter­valles se recouvrent :

• findOverlaps(gr1, gr2) : détec­ter les régions où les inter­valles se recouvrent. Cette fonc­tion retourne un objet de type Hits qui indique pour chaque inter­valle du pre­mier objet (que­ry) les inter­valles du deuxième objet qui le recouvrent (sub­ject).
• countOverlaps(gr1, gr2) pour chaque inter­valle du pre­mier objet, indique le nombre d'intervalles du deuxième objet qui le recouvrent.

Pour finir, voyons une der­nière fonc­tion qui per­met de cal­cu­ler la dis­tance entre deux objets. Créons tout d'abord deux nou­veaux objets :

• distanceToNearest(gr3, gr4) : cal­cule la dis­tance entre chaque inter­valle du pre­mier objet et les inter­valles du deuxième objet.

Importer et exporter des fichiers

La biblio­thèque rtra­ck­layer per­met d'importer et d'exporter des don­nées sous forme d'objet granges.

Atten­tion ! Si vous sou­hai­tez que la colonne de méta­don­nées soit expor­tée, elle doit por­ter l'en-tête name. En effet, les colonnes nom­mées name et score seront expor­tées. Pen­sez bien à véri­fier que c'est le cas avant d'exporter votre objet !

Conclusion

Nous avons vu, dans les deux pre­miers articles de cette série, com­ment créer et mani­pu­ler des inter­valles géno­miques avec R, à l'aide des biblio­thèques IRanges et Geno­mi­cRanges.

Dans le pro­chain article, nous ver­rons com­ment effec­tuer le même genre d'opérations tout en uti­li­sant une syn­taxe "tidy­verse" à l'aide de la biblio­thèque ply­ranges.

Pour aller plus loin …

• Page d'accueil de la biblio­thèque Geno­mi­cRanges (Bio­con­duc­tor)
• An intro­duc­tion to the Geno­mi­cRanges package
• Page d'accueil de la biblio­thèque BSge­nome (Bio­con­duc­tor)
• Page d'accueil de la biblio­thèque rtra­ck­layer (Bio­con­duc­tor)

---

Un grand mer­ci à Guillaume, Mat­thias et Pierre pour leurs com­men­taires, conseils et autres remarques construc­tives !

Les logos des biblio­thèques ont été télé­char­gés sur BiocS­ti­ckers.