Search
Question: ggbio: examples of overlaying multiple assay datasets
0
gravatar for Abhishek Pratap
3.7 years ago by
United States
Abhishek Pratap170 wrote:

Hi Tengfei and others

Just wondering if you know of some code based examples where people
have used ggbio based circular plots for overlaying different
expression based datasets.

I have RNA, miRNA and methylation datasets that I would love to
present it users in a single plot using some sort of a circular plot
as my first attempt.

pointers appreciated.

(i did send this earlier on the mailing list but assuming this forum is the new way to go).. I do believe this would be a very useful resource ..look fwd to reusable BioC help.


Thanks!
-Abhi

ADD COMMENTlink modified 3.7 years ago by Panagiotis Moulos30 • written 3.7 years ago by Abhishek Pratap170
2
gravatar for Panagiotis Moulos
3.7 years ago by
Greece
Panagiotis Moulos30 wrote:

Hello Abhi. ggbio plots are highly customizable and maybe is hard to find an example that does exactly what you want. A very useful starting point is the vignette of ggbio.

Now, the following example does not display data in circular plots but in a karyogram (which may also fit your needs). However, it illustrates the logic behind ggbio and it should not be very hard to adjust it for your data and circular plots.

# Load requried libraries
library(ggbio)
library(GenomicRanges)

# The custom ggbiosave function
ggbiosave <- function (filename = default_name(plot), plot = last_plot(), 
    device = default_device(filename), path = NULL, scale = 1, 
    width = par("din")[1], height = par("din")[2], units = c("in", 
        "cm", "mm"), dpi = 300, limitsize = TRUE, ...) 
{
    #if (!inherits(plot, "ggplot") & !is(plot, "Tracks")) 
    #    stop("plot should be a ggplot2 plot or tracks object")
    eps <- ps <- function(..., width, height) grDevices::postscript(..., 
        width = width, height = height, onefile = FALSE, horizontal = FALSE, 
        paper = "special")
    tex <- function(..., width, height) grDevices::pictex(..., 
        width = width, height = height)
    pdf <- function(..., version = "1.4") grDevices::pdf(..., 
        version = version)
    svg <- function(...) grDevices::svg(...)
    wmf <- function(..., width, height) grDevices::win.metafile(..., 
        width = width, height = height)
    emf <- function(..., width, height) grDevices::win.metafile(..., 
        width = width, height = height)
    png <- function(..., width, height) grDevices::png(..., width = width, 
        height = height, res = dpi, units = "in")
    jpg <- jpeg <- function(..., width, height) grDevices::jpeg(..., 
        width = width, height = height, res = dpi, units = "in")
    bmp <- function(..., width, height) grDevices::bmp(..., width = width, 
        height = height, res = dpi, units = "in")
    tiff <- function(..., width, height) grDevices::tiff(..., 
        width = width, height = height, res = dpi, units = "in")
    default_name <- function(plot) {
        paste(digest.ggplot(plot), ".pdf", sep = "")
    }
    default_device <- function(filename) {
        pieces <- strsplit(filename, "\\.")[[1]]
        ext <- tolower(pieces[length(pieces)])
        match.fun(ext)
    }
    units <- match.arg(units)
    convert_to_inches <- function(x, units) {
        x <- switch(units, `in` = x, cm = x/2.54, mm = x/2.54/10)
    }
    convert_from_inches <- function(x, units) {
        x <- switch(units, `in` = x, cm = x * 2.54, mm = x * 
            2.54 * 10)
    }
    if (!missing(width)) {
        width <- convert_to_inches(width, units)
    }
    if (!missing(height)) {
        height <- convert_to_inches(height, units)
    }
    if (missing(width) || missing(height)) {
        message("Saving ", prettyNum(convert_from_inches(width * 
            scale, units), digits = 3), " x ", prettyNum(convert_from_inches(height * 
            scale, units), digits = 3), " ", units, " image")
    }
    width <- width * scale
    height <- height * scale
    if (limitsize && (width >= 50 || height >= 50)) {
        stop("Dimensions exceed 50 inches (height and width are specified in inches/cm/mm, not pixels).", 
            " If you are sure you want these dimensions, use 'limitsize=FALSE'.")
    }
    if (!is.null(path)) {
        filename <- file.path(path, filename)
    }
    device(file = filename, width = width, height = height, ...)
    on.exit(capture.output(dev.off()))
    print(plot)
    invisible()
}

# Firstly, we retrieve chromosome information for the organism of initerest to be used
# as seqlengths in our custom datasets. We will come back to this data frame holding
# chromosomal lengths later.
download.file(paste("http://hgdownload.cse.ucsc.edu/goldenPath/mm9/database/chromInfo.txt.gz",sep=""),
  file.path(tempdir(),"chromInfo.txt.gz"))
if (.Platform$OS.type == "unix")
    system(paste("gzip -df",file.path(tempdir(),"chromInfo.txt.gz")))
else
    unzip(file.path(tempdir(),"chromInfo.txt.gz"))
chrom.info <- read.delim(file.path(tempdir(),"chromInfo.txt"),header=FALSE)
chrom.info <- chrom.info[-grep("rand|chrM|hap",chrom.info[,1]),1:2]
rownames(chrom.info) <- chrom.info[,1]

# Then, download cytoband information for the organism of interest to be used for coloring
# the bands of each chromosome
download.file(paste("http://hgdownload.soe.ucsc.edu/goldenPath/mm9/database/cytoBand.txt.gz",sep=""),
  file.path(tempdir(),"cytoBand.txt.gz"))
if (.Platform$OS.type == "unix")
    system(paste("gzip -df",file.path(tempdir(),"cytoBand.txt.gz")))
else
    unzip(file.path(tempdir(),"chromInfo.txt.gz"))

# Create the cytoband information to be used for the respective layer in ggbio
cyto.text <- read.delim(file.path(tempdir(),"cytoBand.txt"),header=FALSE)
names(cyto.text) <- c("chromosome","start","end","cytoband","gieStain")
cyto.info <- makeGRangesFromDataFrame(
    df=cyto.text,
    keep.extra.columns=TRUE,
    seqnames.field="chromosome"
)
seqlengths(cyto.info) <- chrom.info[names(seqlengths(cyto.info)),2]
cyto.info <- keepSeqlevels(cyto.info, paste("chr",c(1:19,"X","Y"),sep=""))

# Now, you have to find a way to format your data in a GRanges object and think
# how you can represent what you want to show in an axes-system suitable for layered
# visualization in ggbio. In the following code, I use a hypothetical file which
# contains CGH aberration data. The first three columns contain genomic co-ordinates
# (chromosome, start, end), the fourth column contains 1 if the respective line
# represents an amplification, or -1 for a deletion. Finally, the fifth column
# contain the name of the aCGH array into which this aberration was discovered.
cgh.data <- read.delim("cgh_aberrations.txt")

# This variable makes sure that there is some space across points in the vertical
# axis that will be created above each chromosome in the karyogram
add <- 0

# Loop through the arrays of our aCGH experiment
for (arr in unique(as.character(cgh.data$array))) {
    ind <- which(cgh.data$array==arr)
    for (i in ind) {
        # Transform 1s and -1s so as to be distinctly visible in the final graphic
        # Maybe not clear now, you have to play with it to understand
        if (cgh.data$aberration[i]==1)
            cgh.data$aberration[i] <- cgh.data$aberration[i] + (3*add + 1)
        else
            cgh.data$aberration[i] <- cgh.data$aberration[i] - (3*add + 1)
    }
    add <- add + 1
}
# To separate amplifications from deletions, we have to draw a horizontal line
# somewhere, which adds one more layer in our final plot.
cgh.data$aberration <- cgh.data$aberration + 55

# After all these transformations, create a GRanges object with the data to display
cgh.data.gr <- makeGRangesFromDataFrame(
    df=cgh.data,
    keep.extra.columns=TRUE,
    seqnames.field="chromosome"
)

# Back to chromosomal lengths. We will create a GRanges object with chromosomal
# lengths and additional data so that we can draw the horizontal line that will
# separate amplifications from deletions.
seqls <- as.numeric(chrom.info[,2])
names(seqls) <- rownames(chrom.info)
chrom.info.gr <- GRanges(
    seqnames=rep(as.character(rownames(chrom.info)),each=2),
    IRanges(
        start=c(
            1,chrom.info[1,2],
            1,chrom.info[2,2],
            1,chrom.info[3,2],
            1,chrom.info[4,2],
            1,chrom.info[5,2],
            1,chrom.info[6,2],
            1,chrom.info[7,2],
            1,chrom.info[8,2],
            1,chrom.info[9,2],
            1,chrom.info[10,2],
            1,chrom.info[11,2],
            1,chrom.info[12,2],
            1,chrom.info[13,2],
            1,chrom.info[14,2],
            1,chrom.info[15,2],
            1,chrom.info[16,2],
            1,chrom.info[17,2],
            1,chrom.info[18,2],
            1,chrom.info[19,2],
            1,chrom.info[20,2],
            1,chrom.info[21,2]
        ),
        end=rep(as.numeric(chrom.info[,2]),each=2)
    ),
    seqlengths=seqls
)
elementMetadata(chrom.info.gr) <- DataFrame(zero=rep(55,length(chrom.info.gr)))

# Now, create the final GRanges object which will hold the scores which you want
# to display in a karyogram. Notice that the main difference with cgh.data.gr is
# a change in the coordinates which place the start and end in the same location
# (middle of aberration), otherwise the result will be very blurry. We also add
# seqlengths which is necessary for the plot.
score.layer <- GRanges(
    seqnames=seqnames(cgh.data.gr),
    IRanges(
        start=round(start(cgh.data.gr) + (end(cgh.data.gr) - start(cgh.data.gr) + 1)/2),
        end=round(start(cgh.data.gr) + (end(cgh.data.gr) - start(cgh.data.gr) + 1)/2),
    ),
    seqlengths=seqlengths(cyto.info)
)
elementMetadata(score.layer) <- elementMetadata(cgh.data.gr)

# and finally, create the karyogram!
p <- ggplot(cyto.info) + 
    theme_bw() +
    layout_karyogram(cytoband=TRUE,rect.height=8) +
    guides(fill=FALSE) +
    layout_karyogram(chrom.info.gr,aes(x=start,y=zero),geom="line",ylim=c(15,100),size=0.1,guide="none") +
    layout_karyogram(score.layer,aes(x=start,y=aberration,color=array),geom="point",size=0.75,ylim=c(15,100)) +
    facet_wrap(~seqnames,nrow=7) +
    theme(axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.x=element_text(angle=45,hjust=1,vjust=1.2),
        legend.key=element_blank(),
        legend.title=element_blank(),
        legend.position="bottom") +
    guides(colour=guide_legend(override.aes=list(size=2),nrow=2,byrow=TRUE,keyheight=unit(0.5,"lines")))
ggbiosave(filename=file.path(the.path,"cgh_karyogram_multi.pdf"),plot=p,width=210,height=297,
    dpi=600,units="mm")

The above example will produce something similar to the following plot (also here):

Example karyogram with data

ADD COMMENTlink modified 3.7 years ago • written 3.7 years ago by Panagiotis Moulos30
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 142 users visited in the last hour