ggbio - Visualize genomic data

Gene model from OrganismDb object

OrganismDb object has a simpler API to retrieve data from different annotation resources, so we could label our transcripts in different ways.

library(ggbio)
library(Homo.sapiens)
#load gene symbol : GRanges, one gene/row
data(genesymbol, package = "biovizBase")
#retrieve information of the gene of interest
wh <- genesymbol[c("BRCA1", "NBR1")]
wh <- range(wh, ignore.strand = TRUE)
#Plot the different transcripts  for our gene of interest
p.txdb <- autoplot(Homo.sapiens, which = wh)
p.txdb
#Change inton geometry, use gap.geom
autoplot(Homo.sapiens, which = wh, gap.geom = "chevron")

Different arguments to change colors : label.color (color of the label), color(line color) and fill (fill color of exons):

autoplot(Homo.sapiens, which = wh, label.color = "black", color = "brown",
fill = "brown")

Label could be turned off by setting it to FALSE, you could also use expression to make a flexible label combination from column names.

columns(Homo.sapiens)

##  [1] "GOID"         "TERM"         "ONTOLOGY"     "DEFINITION"  
##  [5] "ENTREZID"     "PFAM"         "IPI"          "PROSITE"     
##  [9] "ACCNUM"       "ALIAS"        "CHR"          "CHRLOC"      
## [13] "CHRLOCEND"    "ENZYME"       "MAP"          "PATH"        
## [17] "PMID"         "REFSEQ"       "SYMBOL"       "UNIGENE"     
## [21] "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS" "GENENAME"    
## [25] "UNIPROT"      "GO"           "EVIDENCE"     "GOALL"       
## [29] "EVIDENCEALL"  "ONTOLOGYALL"  "OMIM"         "UCSCKG"      
## [33] "CDSID"        "CDSNAME"      "CDSCHROM"     "CDSSTRAND"   
## [37] "CDSSTART"     "CDSEND"       "EXONID"       "EXONNAME"    
## [41] "EXONCHROM"    "EXONSTRAND"   "EXONSTART"    "EXONEND"     
## [45] "GENEID"       "TXID"         "EXONRANK"     "TXNAME"      
## [49] "TXCHROM"      "TXSTRAND"     "TXSTART"      "TXEND"

#Flexible label
autoplot(Homo.sapiens, which = wh, columns = c("GENENAME", "GO"), names.expr = "GENENAME::GO")

Gene model from TranscriptDb object

TranscriptDb doesn't contain any gene symbol information, so we use tx id as default for label.

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
autoplot(txdb, which = wh)

Create a bam object

RSamtools package is required. The following code is just an example to create a bam object. This bam object can be used in autoplot function

bam<-BamFile(file="file.bam", index="file.bai")
#use it for autoplot
autoplot(bam, which = wh)

Visualize bam file

ggbio supports visualization of alignments file stored in bam, autoplot method accepts :

• bam file path (indexed)
• BamFile object

It's simple to just pass a file path to autoplot function, you can stream a chunk of region by providing 'which' parameter. Otherwise please use method 'estiamte' to show overall estiamted coverage.

fl.bam <- system.file("extdata", "wg-brca1.sorted.bam", package = "biovizBase")
#keeps only the seqlevels in value and removes all others
wh <- keepSeqlevels(wh, "chr17")
autoplot(fl.bam, which = wh)

Mismatch proportion

To show mismatch proportion, you have to provide reference sequence, the mismatched proportion is color coded in the bar chart.

library(BSgenome.Hsapiens.UCSC.hg19)
bg <- BSgenome.Hsapiens.UCSC.hg19
p.mis <- autoplot(fl.bam, bsgenome = bg, which = wh, stat = "mismatch")
p.mis

View all coverage distribution

To view overall estimated coverage distribution, please use method ‘estiamte’. ‘which’ parameter also accept characters. And there is a hidden value called ‘..coverage..’ to let you do simple transformation in aes().

#View all coverage distribution
autoplot(fl.bam, method = "estimate")
#Select chromosomes of interest
#Log transformation of coverage
autoplot(fl.bam, method = "estimate", 
         which = paste0("chr", 17:18),
         aes(y = log(..coverage..)))

Load data

#Load the data
data("CRC", package = "biovizBase")
head(hg19sub)

## GRanges with 6 ranges and 0 metadata columns:
##       seqnames         ranges strand
##                  
##   [1]        1 [1, 249250621]      *
##   [2]        2 [1, 243199373]      *
##   [3]        3 [1, 198022430]      *
##   [4]        4 [1, 191154276]      *
##   [5]        5 [1, 180915260]      *
##   [6]        6 [1, 171115067]      *
##   ---
##   seqlengths:
##            1         2         3 ...        20        21        22
##    249250621 243199373 198022430 ...  63025520  48129895  51304566

Create ideogram, label and scale track

The function layouts the circle by the order you created from inside to outside.

p <- ggbio() + circle(hg19sub, geom = "ideo", fill = "gray70") + #Ideogram
circle(hg19sub, geom = "scale", size = 2) + #Scale
circle(hg19sub, geom = "text", aes(label = seqnames), vjust = 0, size = 3) # label
p #print plot

Show somatic mutation

We add a “rectangle” track to show somatic mutation.

head(mut.gr)

## GRanges with 6 ranges and 10 metadata columns:
##       seqnames                 ranges strand | Hugo_Symbol Entrez_Gene_Id
##                           |          
##   [1]        1 [ 11003085,  11003085]      + |      TARDBP          23435
##   [2]        1 [ 62352395,  62352395]      + |       INADL          10207
##   [3]        1 [194960885, 194960885]      + |         CFH           3075
##   [4]        2 [ 10116508,  10116508]      - |        CYS1         192668
##   [5]        2 [ 33617747,  33617747]      + |     RASGRP3          25780
##   [6]        2 [ 73894280,  73894280]      + |     C2orf78         388960
##         Center NCBI_Build   Strand Variant_Classification Variant_Type
##                              
##   [1]    Broad         36        +               Missense          SNP
##   [2]    Broad         36        +               Missense          SNP
##   [3]    Broad         36        +               Missense          SNP
##   [4]    Broad         36        -               Missense          SNP
##   [5]    Broad         36        +               Missense          SNP
##   [6]    Broad         36        +               Missense          SNP
##       Reference_Allele Tumor_Seq_Allele1 Tumor_Seq_Allele2
##                                   
##   [1]                G                 G                 A
##   [2]                T                 T                 G
##   [3]                G                 G                 A
##   [4]                C                 C                 T
##   [5]                C                 C                 T
##   [6]                T                 T                 C
##   ---
##   seqlengths:
##            1         2         3 ...        20        21        22
##    249250621 243199373 198022430 ...  63025520  48129895  51304566

p <- ggbio() + circle(mut.gr, geom = "rect", color = "steelblue") + #somatic mutation
circle(hg19sub, geom = "ideo", fill = "gray70") +#Ideogram
circle(hg19sub, geom = "scale", size = 2) +#Scale
circle(hg19sub, geom = "text", aes(label = seqnames), vjust = 0, size = 3)#label
p

Many other examples are available in the ggbio package vignette