library(tsTools)
require(TxDb.Dmelanogaster.UCSC.dm6.ensGene)
library(RColorBrewer)
library(GEOquery)
library(rtracklayer)

Installation

tsTools is available on github (https://github.com/musikutiv/tsTools). Documentation is provided as Vignette.

install.packages("devtools")
library(devtools)
install_github("musikutiv/tsTools", build_vignettes=T)

library(tsTools)
browseVignettes("tsTools”)

the function plotProfiles takes the following parameters (bold = required)

| name | explanation | |------------|---------------------------------------------------------------------------------| fstart | start of the genomic window (bp) fend | end of the genomic window (bp) fchr | chromosome of the genomic window profs | named list of coverages i.e. list of SimpleRleLists one per profile to be plotted cols | colors of the profiles (= vector of length profs) ann | a named list dataframe of annotations. Dataframes with colnames chr, start, end, col. The col column specifies color if individual data points ylabel | the label of the y axis ylims | list of ylimits for plotting the coverages. list(c(min, max), c(min, max)) txdb | a transcription database used for plotting gene models (TxDb format) ftitle | title of the plot collapse | collapse gene models (default = TRUE) with.genes.highlited | vector of gene ids that should be highlighted plot.labels | plot labels for genes (default = TRUE) grid | plot grid (default = FALSE) with.average | plot average (default = FALSE)

a simple example, one profile

I load the profile from the packages. For sake of space the profile just consists of one chromosome. Therefore I have to explictly re-create a SimpleRleList object. Typically a Granges to coverage conversion would yield a SimpleRleList for direct usage.

The genomic annotations are taken from Bioconductors TxDb.Dmelanogaster.UCSC.dm6.ensGene library. Alternatively it can be created from GTF or GFF files using the GenomicFeatures package.

cov <- new("SimpleRleList")
fpath <- system.file("extdata", "covx.rds", package="tsTools")
cov[["chrX"]] <- readRDS(fpath)

plotProfiles(fstart=1660000,fend=1720000,fchr="chrX", profs=list(MSL2=cov), cols=c("red"),txdb=TxDb.Dmelanogaster.UCSC.dm6.ensGene)

three profiles

In this example I simply re-use the same profile

cov <- new("SimpleRleList")
fpath <- system.file("extdata", "covx.rds", package="tsTools")
cov[["chrX"]] <- readRDS(fpath)

plotProfiles(fstart=1660000,fend=1720000,fchr="chrX", profs=list(MSL2=cov, MSL3=cov, MSL4=cov), cols=brewer.pal(3, "Set1"),txdb=TxDb.Dmelanogaster.UCSC.dm6.ensGene)

adding annotations

fpath <- system.file("extdata", "msl2.bed", package="tsTools")
peaks <- read.delim(fpath, header=F)
peaks[,1] <- paste0("chr", peaks[,1])

ann <- list("MSL2 peaks"=data.frame(chr=peaks[,1], start=peaks[,2], end=peaks[,3], col=c("blue")))

plotProfiles(fstart=1660000,fend=1720000,fchr="chrX", profs=list(MSL2=cov), cols=c("red"),txdb=TxDb.Dmelanogaster.UCSC.dm6.ensGene, ann=ann)

getting profiles from GEO (GSM1941084)

Using the geoQuery package download a bedgraph file from GEO. rtracklayer is then used to import the track which yieals a GRanges list. Which is then converted to a SimpleRleList required for the plotProfiles function. The chromosome names of the coverage has to be adjusted as they don't match the ones of the txdb annotation.

filePaths = getGEOSuppFiles("GSM1941084")
granges <- import(rownames(filePaths)[1])
granges
cov <- coverage(granges, weight = "score")
names(cov) <- paste0("chr",names(cov))
plotProfiles(fstart=1660000,fend=1720000,fchr="chrX", profs=list(MSL2=cov), cols=c("red"),txdb=TxDb.Dmelanogaster.UCSC.dm6.ensGene, ann=ann)

Alternatives

https://github.com/rstats-gsoc/gsoc2017/wiki/Interactive-Genome-Browser-in-R http://www.bioconductor.org/packages/release/bioc/html/GenomeGraphs.html https://bioconductor.org/packages/release/bioc/vignettes/Gviz/inst/doc/Gviz.pdf https://bioconductor.org/packages/release/bioc/html/ggbio.html

sessionInfo()


musikutiv/tsTools documentation built on June 6, 2023, 1:39 a.m.