Introduction
Nous avons abordé, dans le précédent article de cette série, les bases de la manipulation d'intervalles dans R.
Ce deuxième article a pour objectif de montrer comment manipuler des intervalles génomiques. Au niveau le plus basique, un intervalle est défini par deux nombres entiers positifs délimitant son début et sa fin. Nous allons pouvoir ajouter des couches supplémentaires d'informations, telles que les chromosomes et leurs tailles ainsi que le brin.
Nous allons pour cela utiliser une bibliothèque de la suite Bioconductor, GenomicRanges.
Cette bibliothèque est en quelque sorte l'équivalent des fameux bedtools, très utilisés pour les analyses en génomique.
Installer et charger la bibliothèque GenomicRanges
Commençons par installer et charger la bibliothèque.
|
1 2 3 4 5 6 7 8 9 |
# 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 contenant des intervalles, ainsi que des informations sur les chromosomes et les brins de chromosomes qui portent ces intervalles. Un objet granges est l'équivalent d'un fichier au format .bed (Browser Extensible Data) contenant au minimum trois colonnes : le nom de la séquence, un intervalle (un objet qui comprend une position de début et une position de fin) et le brin (positif "+" ou négatif "-").
|
1 2 3 4 5 |
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-dessus, la fonction c() permet de combiner des éléments dans un vecteur (ou une colonne). La fonction rep(x, n) permet de répéter 'n' fois l'objet x. Ainsi, avec le code seqnames = c(rep("chr1", 4), rep("chr2", 4)), nous assignons les quatre premiers intervalles au chromosome 1 (chr1), et les quatre intervalles suivants au chromosome 2 (chr2).
On indique aussi que la moitié de nos intervalles sont sur le brin '+' et le reste sur le brin '-'. Si le brin n'a pas d'importance, on utilise alors strand = "*" lors de la création de l'objet granges. L'argument ignore.strand = TRUE permet d'ignorer le brin dans la plupart des fonctions de la bibliothèque GenomicRanges.
Voyons à quoi ressemble notre nouvel objet :
Métadonnées
Notre objet gr ne contient pour le moment que des informations nécessaire pour la création des intervalles ("GRanges object with 8 ranges and 0 metadata columns"). Nous pouvons ajouter des informations supplémentaires comme des noms de gènes, le taux en GC, etc. à l'aide du paramètre mcols().
|
1 2 3 4 5 6 7 |
# 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 ressemble maintenant à ceci :
Informations sur les chromosomes
Pour l'instant, les intervalles de notre objet sont situés sur des brins de chromosomes, mais nous n'avons pas précisé leurs tailles, ni le génome ("2 sequences from an unspecified genome; no seqlengths"). Sans ces informations, nous ne pouvons faire des opérations que sur les régions de chromosomes couvertes par ces intervalles.
La commande seqinfo() permet d'afficher et modifier les informations liées aux séquences :
Les informations concernant la taille des chromosomes de l'espèce ou du génome de référence peuvent être trouvées sur le site de l'UCSC ou celui d'Ensembl. La bibliothèque BSgenomes, permet aussi de récupérer ces informations pour certains génomes.
Nous pouvons aussi ajouter ces informations de la façon suivante :
|
1 2 3 4 5 6 7 8 |
# 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" |
Affichons à nouveau les informations à l'aide de seqinfo() :
Extraire des données
Nous pouvons accéder aux informations à l'aide des commandes de base suivantes :
|
1 2 3 4 5 6 7 8 |
# extraire les chromosomes seqnames(gr) # extraire les intervalles génomiques granges(gr) # extraire les brins strand(gr) |
Nous pouvons également utiliser des opérateurs logiques :
|
1 2 3 4 5 6 7 8 9 10 11 |
# 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 pouvons effectuer des opérations sur notre objet :
reduce(gr): fusionner les intervalles chevauchants ou adjacentsgaps(gr): extraire les positions non couvertes par des intervallescoverage(gr): calculer la couverture (nombre de fois que chaque position est couverte par un intervalle)
D'autres fonctions existent :
flank(gr, width = n, start = TRUE): extraire une région de longueur 'n' en amont. Pour extraire une région de longueur 'n' en aval, utiliserstart = FALSE.shift(gr, shift = n): déplacer les intervalle de 'n' positions vers la droite (ne tient pas compte du brin). Pour décaler les intervalles vers la gauche, utilisershift = -n.resize(gr, width = n, fix = "start"): modifier la taille des intervalles de longueur '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 promotrices, entre x positions en amont (par défaut 2000) et y positions en aval (par défaut 200)
Plusieurs objets granges
Après avoir vu différentes opérations sur un seul objet granges, voyons comment comparer plusieurs objets.
Commençons par créer deux objets : gr1 et gr2.
|
1 2 3 4 5 6 7 8 9 |
# 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 fonctions suivantes permettent d'extraire l'union et l'intersection de deux objets granges :
union(gr1, gr2): union des deux objets.intersect(gr1, gr2): intersection des deux objets.
Nous pouvons également rechercher et compter les régions où les intervalles se recouvrent :
findOverlaps(gr1, gr2): détecter les régions où les intervalles se recouvrent. Cette fonction retourne un objet de typeHitsqui indique pour chaque intervalle du premier objet (query) les intervalles du deuxième objet qui le recouvrent (subject).countOverlaps(gr1, gr2)pour chaque intervalle du premier objet, indique le nombre d'intervalles du deuxième objet qui le recouvrent.
Pour finir, voyons une dernière fonction qui permet de calculer la distance entre deux objets. Créons tout d'abord deux nouveaux objets :
|
1 2 |
gr3 <- gr2[1, ] gr4 <- gr1[2, ] |
distanceToNearest(gr3, gr4): calcule la distance entre chaque intervalle du premier objet et les intervalles du deuxième objet.
Importer et exporter des fichiers
La bibliothèque rtracklayer permet d'importer et d'exporter des données sous forme d'objet granges.
Attention ! Si vous souhaitez que la colonne de métadonnées soit exportée, elle doit porter l'en-tête name. En effet, les colonnes nommées name et score seront exportées. Pensez bien à vérifier que c'est le cas avant d'exporter votre objet !
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
# 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 premiers articles de cette série, comment créer et manipuler des intervalles génomiques avec R, à l'aide des bibliothèques IRanges et GenomicRanges.
Dans le prochain article, nous verrons comment effectuer le même genre d'opérations tout en utilisant une syntaxe "tidyverse" à l'aide de la bibliothèque plyranges.
Pour aller plus loin ...
- Page d'accueil de la bibliothèque GenomicRanges (Bioconductor)
- An introduction to the GenomicRanges package
- Page d'accueil de la bibliothèque BSgenome (Bioconductor)
- Page d'accueil de la bibliothèque rtracklayer (Bioconductor)
Un grand merci à Guillaume, Matthias et Pierre pour leurs commentaires, conseils et autres remarques constructives !
Les logos des bibliothèques ont été téléchargés sur BiocStickers.
