-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathbrowsR-vignette.Rmd
122 lines (87 loc) · 4.41 KB
/
browsR-vignette.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
---
title: "The browsR"
author: "Tobias Straub"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Genome coverage plots}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r, warning=F, message=F}
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.
```{r, eval=F}
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.
```{r, message=F, fig.width=6, fig.height=2}
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
```{r, fig.height=3, fig.width=6}
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
```{r, fig.height=2, fig.width=6}
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.
```{r, fig.height=2, fig.width=6, message=F, warning=F}
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
```{r}
sessionInfo()
```