Si vous analysez des données d'épigénomique telles que de l'ATAC-seq ou des ChIP-seq, vous souhaitez sûrement pouvoir représenter des exemples de pics sous forme de genomic tracks comme on en voit souvent dans les publications.
Lorsque l'on inspecte ses données de coverage (ou couverture), on charge généralement le fichier Bam ou le fichier BigWig dans notre navigateur de génome préféré (pour ma part IGV, mais aussi IGB ou UCSC). Après quelques paramétrages, on peut visualiser et comparer plusieurs échantillons ensemble comme dans l'exemple ci-dessous :
IGV, c'est super pratique pour inspecter ses données, se balader le long du génome et regarder ses signaux, mais c'est moins pratique si l'on veut en faire une figure pour un papier. Il y a bien la possibilité de sauvegarder sa vue en PNG ou en SVG, mais ensuite il faut retoucher le résultat pour pouvoir le monter en figure qualité publication, et ça quand on est bioinformaticien, on aime pas trop…
J'ai longtemps cherché un moyen d'automatiser la génération de genomic tracks pas trop moches avec R. Il existe bien des paquets pour faire cela, mais chacun a ses propres limitations et/ou son esthétique bien à lui. Je ne vais pas vous en faire une liste exhaustive, mais sachez que malgré mes recherches, je n'ai toujours pas trouvé l'outil parfait qui répond à tous mes besoins sans devoir faire des compromis.
Voici une liste des paquets que j'ai regardé (sans forcément les tester) :
Vous l'aurez compris, mon choix s'est porté sur ce dernier, Gviz, pour construire mes figures, et je vais vous partager mes paramètres pour arriver à ce résultat :
Le but de cet article n'étant pas de vous faire un tour exhaustif de Gviz, je vous conseille fortement de lire, voire même de tester l'excellente vignette du paquet qui est disponible depuis bioConductor : https://bioconductor.org/packages/release/bioc/vignettes/Gviz/inst/doc/Gviz.html
Le génome de référence
Tout d'abord, on va commencer par préparer la track d'annotation du génome de référence, celle où l'on voit les gènes.
Pour cela, vous pouvez soit utiliser les génomes déjà tous prêts présent dans bioConductor, ou bien charger le même génome que vous avez utilisé pour votre analyse (option que je préfère, pour une question d'exactitude de la présentation de résultats).
Comme je cherche à vous partager les moments où j'ai eu des difficultés afin de vous simplifier la vie (et aussi pour mon moi du futur), je vais vous montrer un cas un peu compliqué avec un GTF de chez Gencode.
Pour dessiner la track du gène Tram1 chez la souris (mm10), voici le code minimal par défaut :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 |
# Load libraries library("GenomicFeatures") library("Gviz") # Path to the GTF file my_genome <- "path/to/mouse/gencode.vM25.annotation.gtf.gz" # Tram1 gene coordinates my_gene <- GRanges( seqnames="chr1", IRanges( start=13552693, end=13601910 ) ) # Gviz gene track geneTrack <- GeneRegionTrack( my_genome, name = "Genes" ) # Plot the region containing the Tram1 gene plotTracks( geneTrack, from = start(my_gene), to = end(my_gene), chromosome = as.character(seqnames(my_gene)) ) |
Voici ce que l'on obtient :
Pas très probant… Le gène est représenté par de gros rectangles, alors qu'on s'attend à voir des exons et des introns.
Les fichiers GTF de chez Gencode contiennent des entrées pour les gènes (et pas que exon et CDS par exemple), ce qui génère de gros rectangles qui masquent la structure exon/intron des gènes.
Voici une des lignes en question dans mon fichier GTF :
|
1 |
"ENSMUSG00000025935.10"; gene_type "protein_coding"; gene_name "Tram1"; level 2 ; mgi_id "MGI:1919515"; havana_gene "OTTMUSG00000049113.1"; |
Pour s'en débarrasser, l'astuce consiste à charger le génome comme un objet TxDb, lisible par Gviz. Cela va enlever les gros rectangles qui masquent la structure des gènes :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
# Load GTF file as a TxDb object TxDb <- GenomicFeatures::makeTxDbFromGFF(my_genome) # Gviz gene track geneTrack <- GeneRegionTrack( TxDb, name = "Genes" ) # Plot the region containing the Tram1 gene plotTracks( geneTrack, from = start(my_gene), to = end(my_gene), chromosome = as.character(seqnames(my_gene)) ) |
C'est beaucoup mieux, mais il manque encore quelque chose, le nom des gènes ! En temps normal, voici à quoi devrait ressembler le code pour afficher les noms des gènes :
|
1 2 3 4 5 6 7 |
# Gviz gene track geneTrack <- GeneRegionTrack( myGTF, # Show gene names transcriptAnnotation="symbol", name = "Genes" ) |
Mais voilà ce que ça donne avec un objet TxDb, on se retrouve avec les ID Ensembl à la place des noms.
Là encore, il va falloir avoir recours à une astuce. Une fois transformé en TxDb, nous avons perdu l'information du nom des gènes car les symboles ne font pas parti du format TxDb. Il va alors falloir aller les récupérer depuis le fichier GTF :
|
1 2 3 4 5 6 7 8 |
# Load the GTF file with rtracklayer genome_gtf <- rtracklayer::import(my_genome) # Extract gene_id and gene_symbol and remove duplicates gene2symbol <- unique(mcols(genome_gtf)[,c("gene_id","gene_name")]) # Define gene_id as rownames rownames(gene2symbol) <- gene2symbol$gene_id |
Nous obtenons alors une table de correspondance des ID des gènes et de leurs noms respectifs :
|
1 2 3 4 5 6 7 8 9 10 |
> head(gene2symbol) DataFrame with 6 rows and 2 columns gene_id gene_name <character> <character> ENSMUSG00000102693.1 ENSMUSG00000102693.1 4933401J01Rik ENSMUSG00000064842.1 ENSMUSG00000064842.1 Gm26206 ENSMUSG00000051951.5 ENSMUSG00000051951.5 Xkr4 ENSMUSG00000102851.1 ENSMUSG00000102851.1 Gm18956 ENSMUSG00000103377.1 ENSMUSG00000103377.1 Gm37180 ENSMUSG00000104017.1 ENSMUSG00000104017.1 Gm37363 |
Pour afficher les noms des gènes sur la track, il va falloir remplacer les ID par les symboles, comme ceci :
|
1 2 3 4 5 6 7 |
geneTrack <- GeneRegionTrack( TxDb, # Show gene names transcriptAnnotation="symbol", name = "Genes" ) ranges(geneTrack)$symbol <- gene2symbol[ranges(geneTrack)$gene, "gene_name"] |
Enfin, ici nous voyons les transcrits alternatifs des gènes. Ce n'est pas trop gênant quand il y en a peu, mais ça le devient quand il y en a vraiment plein. Nous allons alors réduire l'affichage pour ne montrer qu'un seul modèle de gène avec tous les exons alternatifs :
|
1 2 3 4 5 6 7 8 9 10 |
# Gviz gene track geneTrack <- GeneRegionTrack( TxDb, # Show gene names transcriptAnnotation="symbol", # Collapse all alternative transcripts collapseTranscripts = "meta", name = "Genes" ) ranges(geneTrack)$symbol <- gene2symbol[ranges(geneTrack)$gene, "gene_name"] |
Vous noterez au passage que l'on a perdu les boîtes plus fines qui montrent les UTR des gènes. Comme il existe plusieurs UTR possibles, Gviz fait le choix de ne pas les distinguer et de dessiner les UTR aussi épais que des exons.
Travaillons maintenant un peu l'esthétique de la track, parce que le diable se cache dans les détails. Premièrement, je n'aime pas le jaune pour les boîtes des gènes, je n'aime pas non plus qu'il y ait un contours autours des boîtes, parce qu'ils font apparaître les exons plus gros qu'ils ne le sont réellement, je n'aime pas le fond grisé du titre de la track qui manque cruellement de contraste, et enfin chez la souris les gènes doivent être notés en italique. Voilà à quoi ressemble ma version "pimpée" :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 |
# Gviz gene track geneTrack <- GeneRegionTrack( TxDb, # Collapse all alternative transcripts collapseTranscripts = "meta", # Print gene symbols transcriptAnnotation="symbol", # Name of the track name = "Genes", # Gene name in italic fontface.group="italic", # Remove borders around exons col = 0, # Color of the exons fill = "#585858", # Apply the exon color to the thin line in introns col.line = NULL, # Color of the gene names fontcolor.group= "#333333", # Font size fontsize.group=18 ) ranges(geneTrack)$symbol <- gene2symbol[ranges(geneTrack)$gene, "gene_name"] # Plot the region containing the Tram1 gene plotTracks( geneTrack, from = start(my_gene), to = end(my_gene), chromosome = as.character(seqnames(my_gene)), # Remove the grey background and the white borders of the track name background.title = "transparent", col.border.title="transparent", # Track name color col.title = "#333333" ) |
Le résultat est plus sobre, et c'est ce que l'on veut pour une figure de papier, parce que l'on veut que l'attention du lecteur se fixe sur les données que l'on va présenter avec.
Les fichiers Bed
Les analyses d'ATAC ou ChIP-seq consistent à appeler des pics, c'est à dire à détecter les signaux enrichis d'ouverture de la chromatine ou de site de fixation de protéines sur l'ADN. Le résultats de cette analyse se concrétise en un ficher Bed contenant les coordonnées des pics. Nous pouvons représenter ces régions génomique autours de nos gènes préféré dans Gviz comme ceci :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 |
my_peaks <- rtracklayer::import("example.bed") # Get peaks in the Tram1 gene area my_peak_locus <- subsetByOverlaps(my_peaks, my_gene) peak_track <- AnnotationTrack( my_peak_locus, fill = "#5f4780", col.line = NULL, col = 0, name = "Peaks" ) # Plot the region containing the Tram1 gene plotTracks( # Plot peak track on top of the gene track c(peak_track, geneTrack), from = start(my_gene), to = end(my_gene), chromosome = as.character(seqnames(my_gene)), # Remove the grey background and the white borders of the track name background.title = "transparent", col.border.title="transparent", # Track name color col.title = "#333333" ) |
Les bigWig
Charger les bigWig dans R
Les fichiers bigWig sont souvent lourds (autours de 100 Mb, voire beaucoup plus selon le type d'expérience), alors charger ces fichiers en entier dans R peut vite devenir limitant en terme de RAM. Pour ne pas exploser les ressources de votre machine, il est possible de n'importer que les données d'un locus précis grâce à
|
1 |
rtracklayer |
. Voyons comment charger de façon relativement automatisée des fichiers bigWig.
Dans un premier temps, on va lister les fichiers à notre disposition qui sont rangés ensemble dans un seul dossier. Les fichiers bigWig sont tous nommés comme suit : "Condition_Replicat.bw" afin de faciliter l'extraction des informations depuis leur nom (c'est une habitude que j'ai prise, vous faites peut-être différemment) :
|
1 2 3 4 |
# path to the bigWig files bw_folder <- "path/to/my_bg" # List all the bigWig files that have a ".bw" extention bw_files <- list.files(path=bw_folder, pattern = ".bw") |
|
1 2 3 4 |
> bw_files [1] "Cond-1_REP1.bw" "Cond-1_REP2.bw" "Cond-1_REP3.bw" "Cond-2_REP1.bw" [5] "Cond-2_REP2.bw" "Cond-2_REP3.bw" |
Nous allons ensuite extraire les conditions des noms de fichiers, définir une couleur pour chacune d'entre elles, et importer les données relative à notre gène d'intérêt, ici Tram1 :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
# Extract the condition from the bw file name , i.e. the string before the first "_" conditions <- sapply(strsplit(bw_files, "_"), `[`, 1) colors <- c("#fdb049", "#94d574") names(colors) <- unique(conditions) # Import the signal around the gene of interest from the bigwig files import_bw_file <- function(folder, file, locus){ gr <- rtracklayer::import( paste(folder, file, sep="/"), # Specify the region you want to import as a GRanges object which=locus ) return(gr) } # For each bigWig file, make a GRanges object of the region around our gene of interest gr_list <- lapply(bw_files, function(file){ gr <- import_bw_file(bw_folder, file, my_gene) }) names(gr_list) <- conditions |
Visualiser les bigWig avec Gviz
Nous avons maintenant une liste d'objets GRanges avec les données de couvertures autours de notre gène préféré. Nous pouvons maintenant les représenter avec Gviz :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 |
# Generate the bigWig signal track bigwig_tracks <- function(gr, locus, colors){ # My annotation contains non-regular chromosome names # The following option avoids an error because of the chr names options(ucscChromosomeNames=FALSE) # Get the condition of the sample condition <- unique(gr$condition) # Get the corresponding color color <- colors[condition] # Prepare the track gTrack <- DataTrack( # Granges object range = gr, # Type of graph : histogram type = "hist", # No borders around histogram bars col.histogram=0, # Histogram color fill.histogram=color, # Draw a line at 0 baseline=0, # Baseline color col.baseline=color, # Size of the baseline lwd.baseline=1, # Number of bins for the histogram window=1000, chromosome = as.character(seqnames(locus)), name = unique(gr$condition), # The scale starts at 0 and finish at the maximum of the bw values ylim=c(0, trunc(max(gr$score), digit=4)), # Show a tick only for the maximum value yTicksAt=c(0,trunc(max(gr$score), digit=4)), # Height of the track size=1 ) return(gTrack) } bw_track_list <- lapply(gr_list, function(gr){ bw_track <- bigwig_tracks(gr, my_gene, colors) }) # Plot the region containing the Tram1 gene plotTracks( # Plot bigWig tracks, with the peak and the gene track c(bw_track_list, peak_track, geneTrack), from = start(my_gene), to = end(my_gene), chromosome = as.character(seqnames(my_gene)), # Remove the grey background and the white borders of the track name background.title = "transparent", col.border.title="transparent", # Track name color col.title = "#333333", # Axis color col.axis = "#333333" ) |
Et voilà le résultat :
Superposer les réplicats
Parce que je suis toujours plus exigeante, j'aimerais superposer mes réplicats, pour n'avoir que 2 tracks à montrer, ce qui fera gagner beaucoup de place. Je vais reproduire la fonction "Overlay" de IGV, avec mes réplicats qui consiste à superposer les tracks, et je vais appliquer une légère transparence pour voir les différences entre les réplicats :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 |
# For each condition, make a list of GRanges object for the replicates gr_cond_list <- lapply(unique(conditions), function(cond){ # Get the file names corresponding to the current condition cond_bw_files <- grep(cond, bw_files, value=TRUE) # For each bigWig file, make a GRanges object of the region around our gene of interest lapply(cond_bw_files, function(file){ condition <- sapply(strsplit(file, "_"), `[`, 1) gr <- import_bw_file(bw_folder, file, condition, my_gene) }) }) # The result looks like this : # gr_cond_list : # |-Cond1 : # |-Rep1 # |-Rep2 # |-Rep3 # |-Cond2 : # |-Rep1 # |-Rep2 # |-Rep3 # Generate the bigWig signal track bigwig_tracks_overlay <- function(gr_list, locus, colors){ options(ucscChromosomeNames=FALSE) track_list <- lapply(gr_list, function(gr){ condition <- unique(gr$condition) color <- colors[condition] gTrack <- DataTrack( # Granges object range = gr, # Type of graph : histogram type = "hist", # No borders around histogram bars col.histogram=0, # Histogram color fill.histogram=color, # Make the histogram transparent alpha=0.8, # Prevent the track title to also be transparent alpha.title = 1, # Draw a line at 0 baseline=0, # Baseline color col.baseline=color, # Size of the baseline lwd.baseline=1, # Number of bins for the histogram window=1000, chromosome = as.character(seqnames(locus)), name = unique(gr$condition), # The scale starts at 0 and finish at the maximum of the bw values ylim=c(0, trunc(max(gr$score), digit=4)), # Show a tick only for the maximum value yTicksAt=c(0,trunc(max(gr$score), digit=4)), # Height of the track size=2 ) }) # Overlay the tracks from the same condition gTrack <- OverlayTrack(track_list) return(gTrack) } # Run the bigwig_tracks_overlay function bw_track_list <- lapply(gr_cond_list, function(cond_gr){ bw_track <- bigwig_tracks_overlay(cond_gr, my_gene, colors) }) # Plot the region containing the Tram1 gene plotTracks( # Plot peak track on top of the gene track c(bw_track_list, peak_track, geneTrack), from = start(my_gene), to = end(my_gene), chromosome = as.character(seqnames(my_gene)), # Remove the grey background and the white borders of the track name background.title = "transparent", col.border.title="transparent", # Track name color col.title = "#333333", # Axis color col.axis = "#333333" ) |
Normaliser les valeurs des tracks
Enfin, pour aller encore plus dans les détails, vous noterez que la hauteur des tracks est réglée sur les valeurs maximales des données de couvertures, mais que ces valeurs ne sont pas les mêmes d'une condition à l'autre. Pour faciliter la comparaison des deux conditions, on va fixer la hauteur des tracks à la valeur maximale toutes conditions confondues, pour pouvoir apprécier visuellement les différences entre les deux :
|
1 2 3 4 5 |
# Merge all the GRanges objects as one max_score <- do.call("c", do.call("c", gr_cond_list)) # Get the max score value max_score <- max(max_score$score) |
On va maintenant pouvoir passer cette valeur à notre fonction
|
1 |
bigwig_tracks_overlay |
:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
# Generate the bigWig signal track bigwig_tracks_overlay <- function(gr_list, locus, max_score, colors){ [...] # The scale starts at 0 and finish at the maximum of the bw values ylim=c(0, trunc(max_score, digit=4)), # Show a tick only for the maximum value yTicksAt=c(0,trunc(max_score, digit=4)), [...] ) }) # Overlay the tracks from the same condition gTrack <- OverlayTrack(track_list) return(gTrack) } |
Mettre en évidence des régions ou des pics
On y est presque, réglons les derniers détails. Pour mettre en valeur les régions où il y a des pics, on peut dessiner des rectangles en arrière plan qui traversent les tracks afin de mieux voir leur délimitations. Pour cela, on utilise la fonction
|
1 |
HighlightTrack |
qui prend une liste de positions
|
1 |
start |
et
|
1 |
end |
pour délimiter les régions à mettre en évidence, et qui englobe les tracks sous lesquelles on veut faire passer les rectangles :
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 |
ht <- HighlightTrack( # The highlights will be dranw across all the tracks trackList = c(bw_track_list, peak_track, geneTrack), # Start position of the peaks start = start(my_peak_locus), # End position of the peaks end = end(my_peak_locus), chromosome = as.character(seqnames(my_peak_locus)), # No borders col = 0, # Highlight in grey fill = "#666666", # Make it very transparent alpha=0.1 ) # Plot the region containing the Tram1 gene plotTracks( # Plot highlight track that contains all the tracks ht, from = start(my_gene), to = end(my_gene), chromosome = as.character(seqnames(my_gene)), # Remove the grey background and the white borders of the track name background.title = "transparent", col.border.title="transparent", # Track name color col.title = "#333333", # Axis color col.axis = "#333333" ) |
Et voilà ce que ça donne :
On voit maintenant bien où sont les pics par rapport aux signaux des bigWigs et aussi par rapport au gène. Notez toutefois qu'il y a une limitation à l'usage des highlights si vous exportez vos figures en PNG avec des unités autres que des pixels. Voici ce que ça donne quand on l'exporte en donnant des tailles en cm :
|
1 2 3 4 5 6 7 8 9 10 11 12 |
png("test.png", width=14, height=8, unit="cm", res=200) plotTracks( ht, from = start(my_gene), to = end(my_gene), chromosome = as.character(seqnames(my_gene)), background.title = "transparent", col.border.title="transparent", col.title = "#333333", col.axis = "#333333" ) dev.off() |
Les rectangles gris vers la gauche ne sont plus alignés correctement avec les pics ! Je ne compred pas d'où vient ce problème, alors il faut trouver un autre moyen d'exporter ses figures. Je n'ai pas testé toutes les options, mais l’export en PDF fonctionne bien.
Limitations et conclusion
Comme vous l'avez vu, pour générer des figures de genomic tracks de bonne qualité visuelle, il faut en suer un petit peu, et c'est pour ça que j'ai voulu partager mon code avec vous, afin de vous épargner du temps et des tâtonnements. Bien que très modulable, Gviz réclame pas mal de petits bidouillages à gauche à droite pour obtenir le résultat que l'on souhaite, surtout quand on est un poil maniaque des graphiques comme moi… Je n'ai pas présenté toutes mes astuces ici mais si vous voulez par exemple avoir les titres des tracks à l'horizontale, sachez que l'option
|
1 |
rotation.title=0 |
existe bien, mais n'est que partiellement implémentée et il faudra bidouiller avec la largeur de la colonne de titre et rajouter des espaces dans les noms des titres pour que ça passe…
Néanmoins, je trouve quand même que l'effort en vaut la peine, le résultat est visuellement satisfaisant, et je ne pense pas que les autres paquets R puissent faire mieux à l'heure de l'écriture de cet article.
Laisser un commentaire