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 |
[crayon-67300eecc53fa061714115 ]# 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 |
[crayon-67300eecc53ff881986879 ]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
1 |
c() |
permet de combiner des éléments dans un vecteur (ou une colonne). La fonction
1 |
rep(x, n) |
permet de répéter 'n' fois l'objet x. Ainsi, avec le code
1 |
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
1 |
strand = "*" |
lors de la création de l'objet granges. L'argument
1 |
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
1 |
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
1 |
mcols() |
.
1 2 3 4 5 |
[crayon-67300eecc540e566130627 ]# 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
1 |
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 |
[crayon-67300eecc5413156533287 ]# 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
1 |
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 |
[crayon-67300eecc5418904376623 ]# 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 |
[crayon-67300eecc541b830369881 ]# 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
1
granges
1 |
granges |
Nous pouvons effectuer des opérations sur notre objet :
-
1reduce(gr)
-
1gaps(gr)
-
1coverage(gr)
D'autres fonctions existent :
-
1flank(gr, width = n, start = TRUE)1start = FALSE
-
1shift(gr, shift = n)1shift = -n
-
1resize(gr, width = n, fix = "start")1fix = "start"1fix = "end"
-
1promoters(gr, upstream = x, downstream = y)
Plusieurs objets
1
granges
1 |
granges |
Après avoir vu différentes opérations sur un seul objet
1 |
granges |
, voyons comment comparer plusieurs objets.
Commençons par créer deux objets :
1 |
gr1 |
et
1 |
gr2 |
.
1 2 3 4 5 6 7 8 9 |
[crayon-67300eecc5446766680573 ]# 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
1 |
granges |
:
-
1union(gr1, gr2)
-
1intersect(gr1, gr2)
Nous pouvons également rechercher et compter les régions où les intervalles se recouvrent :
-
1findOverlaps(gr1, gr2)1<strong>Hits</strong>
-
1countOverlaps(gr1, gr2)
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 |
[crayon-67300eecc5453279481038 ]gr3 <- gr2[1, ] gr4 <- gr1[2, ] |
-
1distanceToNearest(gr3, gr4)
Importer et exporter des fichiers
La bibliothèque rtracklayer permet d'importer et d'exporter des données sous forme d'objet
1 |
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 |
[crayon-67300eecc5459141613525 ]# 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.
Laisser un commentaire