log2 transformation

[utnesp]

log2 transformation

[utnesp]

I did what you suggested and here is a suggestion where I implement the use of BiomaRt. You would need to initialise a mart before you begin:

mart = useMart("ENSEMBL_MART_ENSEMBL", dataset='hsapiens_gene_ensembl')

To use other marts, you can check here:
https://github.com/utnesp/Easy-bioMart

Here is the code:

ext_name2wiggleplotr.metadata <- function(ext_name, biomart = mart){
    df <- getBM(
        filters= "external_gene_name",
        attributes= c("ensembl_transcript_id", "ensembl_gene_id", "external_gene_name", "strand", "transcript_appris"),
        values= ext_name,
        mart = biomart)
    colnames(df) <- c("transcript_id", "gene_id", "gene_name", "strand", "transcript_appris")
    return(df)
}

metadata2wiggleplotr.exons <- function(ensembl_transcript_id, biomart = mart){
    df <- getBM(
        filters= "ensembl_transcript_id",
        attributes= c("ensembl_transcript_id", "chromosome_name", "exon_chrom_start", "exon_chrom_end", "strand"),
        values= ensembl_transcript_id,
        mart = biomart)
    
    df.GRanges <- GRanges(seqnames = df$chromosome_name, ranges = IRanges(start = df$exon_chrom_start, end = df$exon_chrom_end), strand = df$strand)
    mcols(df.GRanges)$transcript_id <- df$ensembl_transcript_id

    return(df.GRanges)
}

metadata2wiggleplotr.utr <- function(ensembl_transcript_id, biomart = mart){
df <- getBM(
        filters= "ensembl_transcript_id",
        attributes= c("ensembl_transcript_id", "chromosome_name", "5_utr_start", "5_utr_end", "strand"),
        values= ensembl_transcript_id,
        mart = biomart)

df2 <- getBM(
        filters= "ensembl_transcript_id",
        attributes= c("ensembl_transcript_id", "chromosome_name", "3_utr_start", "3_utr_end", "strand"),
        values= ensembl_transcript_id,
        mart = biomart)

colnames(df) <- c("transcript_id", "chromosome_name", "utr_start", "utr_end", "strand")
colnames(df2) <- colnames(df)

df <- rbind(df, df2)
df <- df[!is.na(df$utr_start) & !is.na(df$utr_end), ]

df.GRanges <- GRanges(seqnames = df$chromosome_name, ranges = IRanges(start = df$utr_start, end = df$utr_end), strand = df$strand)
mcols(df.GRanges)$transcript_id <- df$transcript_id

return(df.GRanges)
}

plotWiggle <- function(gene, log2.transform = F, prior.count = NULL, plot.principal.only = F, heights = c(0.75, 0.25), ...) {
    metadata <- ext_name2wiggleplotr.metadata(gene)
    
    suppressWarnings( if(plot.principal.only == T & !is.na(metadata$transcript_appris)) metadata <- metadata[!metadata$transcript_appris == "", ] )
    
    metadata$transcript_appris <- NULL
    
    # create empty GRangesList for used in creation of exonic and utr regions (for loop)
    all.exons <- GRangesList()
    all.ccds <- GRangesList()
    
    for (i in 1:nrow(metadata)) {
        # get utr regions
        utrs <- metadata2wiggleplotr.utr(metadata$transcript_id[i])
        # get exon regions including utr regions
        exons <- metadata2wiggleplotr.exons(metadata$transcript_id[i])
        
        # convert to GRangesList (required for setdiff)
        utrs <- GRangesList(test = utrs); names(utrs) <- metadata$transcript_id[i]
        exons <- GRangesList(test = exons); names(exons) <- metadata$transcript_id[i]
        
        # find differences beteen exonic and utr regions giving rising to coding regions
        ccds <- setdiff(exons,utrs, ignore.strand = F)
        
        # append transcript exon and utr regions to GRangesList
        all.exons <- append(all.exons, exons, length(all.exons))
        all.ccds <- append(all.ccds, ccds, length(all.ccds))
    }
        
    # Plot coverage 
    if(log2.transform == T) {
        p <- plotCoverage(all.exons, all.ccds, transcript_annotations = metadata, return_subplots_list = T, ...)
        if(is.null(prior.count)) prior.count = 1-min(p$coverage_plot$data$coverage) # add prior count to avoid negative values
        p$coverage_plot$data$coverage <- p$coverage_plot$data$coverage + prior.count
        suppressMessages(p$coverage_plot <- p$coverage_plot + scale_y_continuous(trans=log2_trans())) #scale_y_continuous(trans='log2'))
        suppressMessages(cowplot::plot_grid(p$coverage_plot, p$tx_structure, align = "v", rel_heights = heights, ncol = 1))
    } else {
        plotCoverage(all.exons, all.ccds, transcript_annotations = metadata, ...)    
    }
}

[kauralasoo]

Hi Peter,

Thanks! That's a great suggestion. I will think about the best way of implementing it. I would like to avoid adding too many custom parameters to the plotCoverage function if possible (but log2 might well be one of those that should be in there).

Alternative option would be to refactor the code so that plotCoverage also return the original data frame that is then used by ggplot2 to make the coverage plots. This way users could look at the wiggplotr ggplot code and customise it to fit their needs (e.g. log2 transform the y-axis). Would that be useful for you?

Unfortunately I have a lot of other obligations at the moment so not sure when I will be able to look at the code.

Best,
Kaur

[kauralasoo]

Thanks, Peter!
Looks like you managed to solve your problem?

I did not realise that you could also do this via the return_subplots_list = T option, which is not documented at all in the vignette at the moment. For the next release, I will update the vignette to describe the return_subplots_list option and how it can be used to modify different aspects of the plot. I usually use this to construct more complex composite plots.

Best,
Kaur