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.

# Installer BiocManager
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

# Installer la bibliothèque GenomicRanges
BiocManager::install("GenomicRanges")

# Charger la bibliothèque GenomicRanges
library(GenomicRanges)

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

gr - GRanges(
  seqnames = c(rep("chr1", 4), rep("chr2", 4)),
  ranges = IRanges(start = c(1, 8, 3, 6, 1, 2, 1, 8),
                   end = c(5, 9, 6, 8, 2, 8, 4, 10)),
  strand = c(rep("+", 2), rep("-", 2), rep("+", 2), rep("-", 2)))

Dans le code ci-des­sus, la fonc­tion

per­met de com­bi­ner des élé­ments dans un vec­teur (ou une colonne). La fonc­tion

per­met de répé­ter 'n' fois l'objet x. Ain­si, avec le code

, 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

lors de la créa­tion de l'objet granges. L'argument

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

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

.

# Ajouter des identifiants
mcols(gr)$id - paste0(seqnames(gr), "_",
                       c(rep("plus", 2), rep("minus", 2),
                         rep("plus", 2), rep("minus", 2)), "_",
                         start(gr), "_", end(gr))

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

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 :

# Longueurs des chromosomes
seqlengths(gr) - c("chr1" = 10, "chr2" = 10)

# Chromosomes circulaires ou non
isCircular(gr) - c(FALSE, FALSE)

# Nom du génome
genome(gr) - "mon.genome"

Affi­chons à nou­veau les infor­ma­tions à l'aide de

:

Extraire des données

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

# extraire les chromosomes
seqnames(gr)

# extraire les intervalles génomiques
granges(gr)

# extraire les brins
strand(gr)

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

# Extraire les intervalles sur le chr1
gr[seqnames(gr) == "chr1", ]

# Extraire les intervalles sur les brins positifs
gr[strand(gr) == "+", ]

# Extraire les intervalles sur le brin négatif du chr2
gr[seqnames(gr) == "chr2" & strand(gr) == "-", ]

# Extraire l'intervalle n°3
gr[mcols(gr)$id == "interval3"]

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

, voyons com­ment com­pa­rer plu­sieurs objets.

Com­men­çons par créer deux objets :

et

.

# gr1 : intervalles du chr1 sur le brin positif
gr1 - gr[seqnames(gr) == "chr1" & strand(gr) == "+", ]

# créer l'objet gr2
gr2 - GRanges(
  seqnames = "chr1",
  ranges = IRanges(start = c(2, 6),
                   end = c(3, 10)),
  strand = "+")

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

:

• 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 :

gr3 - gr2[1, ]
gr4 - gr1[2, ]

• 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

.

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 !

# Installer la bibliothèque GenomicRanges
BiocManager::install("GenomicRanges")

# Charger la bibliothèque
library(rtracklayer)

# Ajouter une colonne "name"
mcols(gr)$name - mcols(gr)$id

# Exporter l'objet granges au format bed
export(gr, "gr.bed")

# Importer un fichier bed
fichier_bed - rtracklayer::import("gr.bed")

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.