Comments (3)
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, ...)
}
}
from wiggleplotr.
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
from wiggleplotr.
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
from wiggleplotr.
Related Issues (12)
- Allow free scales for the read coverage plots. HOT 2
- Warn if the region that is being plotted is too wide to avoid running out of memory and crashing R.
- Filter classes were removed from ensembldb HOT 1
- exon labeling HOT 2
- Plot transcripts from GTF HOT 2
- Allow region_coords in plotCoverage while rescale_introns = TRUE
- Set order of transcripts HOT 1
- plotCoverage can't open BW files HOT 16
- plotCoverage for non cds region HOT 1
- Consider support for the D4 file format
- Set the Y-axis max value a little bit more than free-scale in coveragePlot
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from wiggleplotr.