Skip to content

Commit e9fb5ee

Browse files
showteethshowteeth
committed
update doc
1 parent 68e6568 commit e9fb5ee

17 files changed

+311
-79
lines changed

‎NAMESPACE

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,8 @@ importFrom(GenomeInfoDb,seqlengths)
4040
importFrom(GenomeInfoDb,seqnames)
4141
importFrom(GenomeMatrix,matrixPlotter)
4242
importFrom(GenomicAlignments,alphabetFrequencyFromBam)
43+
importFrom(GenomicAlignments,coverage)
44+
importFrom(GenomicAlignments,readGAlignments)
4345
importFrom(GenomicRanges,GRanges)
4446
importFrom(GenomicRanges,end)
4547
importFrom(GenomicRanges,makeGRangesFromDataFrame)
@@ -50,6 +52,7 @@ importFrom(IRanges,IRanges)
5052
importFrom(IRanges,findOverlaps)
5153
importFrom(IRanges,subsetByOverlaps)
5254
importFrom(RColorBrewer,brewer.pal)
55+
importFrom(Rsamtools,ScanBamParam)
5356
importFrom(Rsamtools,indexBam)
5457
importFrom(S4Vectors,"values<-")
5558
importFrom(dplyr,arrange)
@@ -81,6 +84,7 @@ importFrom(ggplot2,geom_point)
8184
importFrom(ggplot2,geom_polygon)
8285
importFrom(ggplot2,geom_rect)
8386
importFrom(ggplot2,geom_segment)
87+
importFrom(ggplot2,geom_step)
8488
importFrom(ggplot2,geom_text)
8589
importFrom(ggplot2,geom_tile)
8690
importFrom(ggplot2,ggplot)

‎NEWS.md

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
## Minor changes
88
* `geom_coverage` and `ggcoverage` supports changing shared type of y-axis scales across facets (with `facet.y.scale` parameter).
99
* `geom_coverage` and `ggcoverage` supports joint view of multiple tracks.
10+
* `LoadTrackFile` loads bam directly instead of converting bam to bigwig and loading the bigwig.
1011

1112
-------------
1213

‎R/LoadTrack.R

Lines changed: 103 additions & 60 deletions
Original file line numberDiff line numberDiff line change
@@ -25,11 +25,11 @@
2525
#'
2626
#' @return A dataframe.
2727
#' @importFrom rtracklayer import
28-
#' @importFrom Rsamtools indexBam
28+
#' @importFrom Rsamtools indexBam ScanBamParam
2929
#' @importFrom utils read.csv
30-
#' @importFrom GenomicAlignments alphabetFrequencyFromBam
30+
#' @importFrom GenomicAlignments alphabetFrequencyFromBam readGAlignments coverage
3131
#' @importFrom GenomicRanges GRanges
32-
#' @importFrom IRanges IRanges
32+
#' @importFrom IRanges IRanges subsetByOverlaps
3333
#' @importFrom magrittr %>%
3434
#' @importFrom dplyr select filter
3535
#' @importFrom BiocParallel register MulticoreParam bplapply
@@ -170,67 +170,110 @@ LoadTrackFile <- function(track.file, track.folder = NULL, format = c("bam", "wi
170170
stop("Please provide region for visualizing single nucleotide!")
171171
}
172172
} else {
173-
# require deeptools
174-
if (is.null(bamcoverage.path)) {
175-
bamcoverage.path <- Sys.which("bamCoverage")
176-
if (bamcoverage.path == "") {
177-
stop("Can not find bamCoverage automatically, please specify the path!")
178-
}
179-
} else {
180-
bamcoverage.path <- bamcoverage.path
181-
}
182-
183173
# get used gr
184174
gr <- PrepareRegion(region = region, gtf.gr = gtf.gr, gene.name = gene.name, gene.name.type = gene.name.type, extend = extend)
185-
# read track file
186-
if (is.null(n.cores) || n.cores == 1) {
187-
track.list <- lapply(track.file, function(x) {
188-
# get basename
189-
track.file.base <- basename(x)
190-
# bigwig file
191-
out.bw.file <- tempfile(fileext = c(".bw"))
192-
# prepare bamCoverage cmd
193-
bamcoverage.cmd <- paste(
194-
bamcoverage.path, "-b", x, "-o", out.bw.file,
195-
"--binSize", bin.size, "--normalizeUsing", norm.method, bc.extra.para
196-
)
197-
# run command
198-
message(paste("Calling bamCoverage: ", bamcoverage.cmd))
199-
bamcoverage.status <- system(bamcoverage.cmd, intern = TRUE)
200-
bamcoverage.status.code <- attr(bamcoverage.status, "status")
201-
if (!is.null(bamcoverage.status.code)) {
202-
stop("Run bamCoverage error!")
203-
}
204-
# import wig, bigwig and bedgraph file
205-
single.track.df <- as.data.frame(rtracklayer::import(out.bw.file, which = gr))
206-
single.track.df$TrackFile <- track.file.base
207-
return(single.track.df)
208-
})
175+
if (norm.method == "None") {
176+
message("Calculate coverage with GenomicAlignments when norm.method is None!")
177+
if (is.null(n.cores) || n.cores == 1) {
178+
track.list <- lapply(track.file, function(x) {
179+
# get basename
180+
track.file.base <- basename(x)
181+
# load track
182+
param <- Rsamtools::ScanBamParam(which = gr)
183+
ga <- GenomicAlignments::readGAlignments(x, param = param)
184+
ga.cov <- GenomicAlignments::coverage(ga)
185+
ga.cov.gr <- GenomicRanges::GRanges(ga.cov)
186+
ga.cov.df <- GenomicRanges::subsetByOverlaps(ga.cov.gr, gr) %>% as.data.frame()
187+
# valid the region
188+
gr.df <- as.data.frame(gr)
189+
ga.cov.df[1, "start"] <- gr.df[1, "start"]
190+
ga.cov.df[nrow(ga.cov.df), "end"] <- gr.df[1, "end"]
191+
# add track file
192+
ga.cov.df$TrackFile <- track.file.base
193+
return(ga.cov.df)
194+
})
195+
} else {
196+
# register
197+
BiocParallel::register(BiocParallel::MulticoreParam(workers = n.cores), default = TRUE)
198+
track.list <- BiocParallel::bplapply(track.file, BPPARAM = BiocParallel::MulticoreParam(), FUN = function(x) {
199+
# get basename
200+
track.file.base <- basename(x)
201+
# load track
202+
param <- Rsamtools::ScanBamParam(which = gr)
203+
ga <- GenomicAlignments::readGAlignments(x, param = param)
204+
ga.cov <- GenomicAlignments::coverage(ga)
205+
ga.cov.gr <- GenomicRanges::GRanges(ga.cov)
206+
ga.cov.df <- GenomicRanges::subsetByOverlaps(ga.cov.gr, gr) %>% as.data.frame()
207+
# valid the region
208+
gr.df <- as.data.frame(gr)
209+
ga.cov.df[1, "start"] <- gr.df[1, "start"]
210+
ga.cov.df[nrow(ga.cov.df), "end"] <- gr.df[1, "end"]
211+
# add track file
212+
ga.cov.df$TrackFile <- track.file.base
213+
return(ga.cov.df)
214+
})
215+
}
209216
} else {
210-
# register
211-
BiocParallel::register(BiocParallel::MulticoreParam(workers = n.cores), default = TRUE)
212-
track.list <- BiocParallel::bplapply(track.file, BPPARAM = BiocParallel::MulticoreParam(), FUN = function(x) {
213-
# get basename
214-
track.file.base <- basename(x)
215-
# bigwig file
216-
out.bw.file <- tempfile(fileext = c(".bw"))
217-
# prepare bamCoverage cmd
218-
bamcoverage.cmd <- paste(
219-
bamcoverage.path, "-b", x, "-o", out.bw.file,
220-
"--binSize", bin.size, "--normalizeUsing", norm.method, bc.extra.para
221-
)
222-
# run command
223-
message(paste("Calling bamCoverage: ", bamcoverage.cmd))
224-
bamcoverage.status <- system(bamcoverage.cmd, intern = TRUE)
225-
bamcoverage.status.code <- attr(bamcoverage.status, "status")
226-
if (!is.null(bamcoverage.status.code)) {
227-
stop("Run bamCoverage error!")
217+
message("Calculate coverage with bamCoverage when norm.method is not None!")
218+
# require deeptools
219+
if (is.null(bamcoverage.path)) {
220+
bamcoverage.path <- Sys.which("bamCoverage")
221+
if (bamcoverage.path == "") {
222+
stop("Can not find bamCoverage automatically, please specify the path!")
228223
}
229-
# import wig, bigwig and bedgraph file
230-
single.track.df <- as.data.frame(rtracklayer::import(out.bw.file, which = gr))
231-
single.track.df$TrackFile <- track.file.base
232-
return(single.track.df)
233-
})
224+
} else {
225+
bamcoverage.path <- bamcoverage.path
226+
}
227+
# read track file
228+
if (is.null(n.cores) || n.cores == 1) {
229+
track.list <- lapply(track.file, function(x) {
230+
# get basename
231+
track.file.base <- basename(x)
232+
# bigwig file
233+
out.bw.file <- tempfile(fileext = c(".bw"))
234+
# prepare bamCoverage cmd
235+
bamcoverage.cmd <- paste(
236+
bamcoverage.path, "-b", x, "-o", out.bw.file,
237+
"--binSize", bin.size, "--normalizeUsing", norm.method, bc.extra.para
238+
)
239+
# run command
240+
message(paste("Calling bamCoverage: ", bamcoverage.cmd))
241+
bamcoverage.status <- system(bamcoverage.cmd, intern = TRUE)
242+
bamcoverage.status.code <- attr(bamcoverage.status, "status")
243+
if (!is.null(bamcoverage.status.code)) {
244+
stop("Run bamCoverage error!")
245+
}
246+
# import wig, bigwig and bedgraph file
247+
single.track.df <- as.data.frame(rtracklayer::import(out.bw.file, which = gr))
248+
single.track.df$TrackFile <- track.file.base
249+
return(single.track.df)
250+
})
251+
} else {
252+
# register
253+
BiocParallel::register(BiocParallel::MulticoreParam(workers = n.cores), default = TRUE)
254+
track.list <- BiocParallel::bplapply(track.file, BPPARAM = BiocParallel::MulticoreParam(), FUN = function(x) {
255+
# get basename
256+
track.file.base <- basename(x)
257+
# bigwig file
258+
out.bw.file <- tempfile(fileext = c(".bw"))
259+
# prepare bamCoverage cmd
260+
bamcoverage.cmd <- paste(
261+
bamcoverage.path, "-b", x, "-o", out.bw.file,
262+
"--binSize", bin.size, "--normalizeUsing", norm.method, bc.extra.para
263+
)
264+
# run command
265+
message(paste("Calling bamCoverage: ", bamcoverage.cmd))
266+
bamcoverage.status <- system(bamcoverage.cmd, intern = TRUE)
267+
bamcoverage.status.code <- attr(bamcoverage.status, "status")
268+
if (!is.null(bamcoverage.status.code)) {
269+
stop("Run bamCoverage error!")
270+
}
271+
# import wig, bigwig and bedgraph file
272+
single.track.df <- as.data.frame(rtracklayer::import(out.bw.file, which = gr))
273+
single.track.df$TrackFile <- track.file.base
274+
return(single.track.df)
275+
})
276+
}
234277
}
235278
}
236279
} else if (format == "txt") {

‎R/geom_coverage.R

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,7 @@
2323
#' @param mark.label.size The label size of mark label. Default: 4.
2424
#'
2525
#' @return Layers of ggplot2.
26-
#' @importFrom ggplot2 aes_string scale_fill_manual geom_rect geom_text aes
26+
#' @importFrom ggplot2 aes_string scale_fill_manual geom_rect geom_text aes geom_step
2727
#' @importFrom rlang .data
2828
#' @importFrom grDevices colorRampPalette
2929
#' @importFrom RColorBrewer brewer.pal

‎R/ggcoverage.R

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,7 @@
2626
#' @return A ggplot2 object.
2727
#' @importFrom magrittr %>%
2828
#' @importFrom ggplot2 ggplot aes_string scale_fill_manual geom_rect geom_text aes theme_classic theme unit
29-
#' element_blank annotate rel scale_y_continuous expansion scale_x_continuous coord_cartesian
29+
#' element_blank annotate rel scale_y_continuous expansion scale_x_continuous coord_cartesian geom_step
3030
#' @importFrom scales comma
3131
#' @importFrom grDevices colorRampPalette
3232
#' @importFrom RColorBrewer brewer.pal

‎README.Rmd

Lines changed: 53 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -113,19 +113,41 @@ gtf.gr = rtracklayer::import.gff(con = gtf.file, format = 'gtf')
113113
```
114114

115115
### Basic coverage
116-
```{r basic_coverage, warning=FALSE, fig.height = 6, fig.width = 12, fig.align = "center"}
116+
The basic coverage plot has **two types**:
117+
118+
* **facet**: Create subplot for every track (specified by `facet.key`). This is default.
119+
* **joint**: Visualize all tracks in a single plot.
120+
121+
#### joint view
122+
```{r basic_coverage_joint, warning=FALSE, fig.height = 4, fig.width = 12, fig.align = "center"}
117123
basic.coverage = ggcoverage(data = track.df, color = "auto",
124+
plot.type = "joint", group.key = "Type",
125+
mark.region = mark.region, range.position = "out")
126+
basic.coverage
127+
```
128+
129+
#### facet view
130+
```{r basic_coverage, warning=FALSE, fig.height = 6, fig.width = 12, fig.align = "center"}
131+
basic.coverage = ggcoverage(data = track.df, color = "auto", plot.type = "facet",
118132
mark.region = mark.region, range.position = "out")
119133
basic.coverage
120134
```
121135

122-
You can also change Y axis style:
136+
#### Custom Y-axis style
137+
**Change the Y-axis scale label in/out of plot region with `range.position`**:
123138
```{r basic_coverage_2, warning=FALSE, fig.height = 6, fig.width = 12, fig.align = "center"}
124-
basic.coverage = ggcoverage(data = track.df, color = "auto",
139+
basic.coverage = ggcoverage(data = track.df, color = "auto", plot.type = "facet",
125140
mark.region = mark.region, range.position = "in")
126141
basic.coverage
127142
```
128143

144+
**Shared/Free Y-axis scale with `facet.y.scale`**:
145+
```{r basic_coverage_3, warning=FALSE, fig.height = 6, fig.width = 12, fig.align = "center"}
146+
basic.coverage = ggcoverage(data = track.df, color = "auto", plot.type = "facet",
147+
mark.region = mark.region, range.position = "in",
148+
facet.y.scale = "fixed")
149+
basic.coverage
150+
```
129151

130152
### Add gene annotation
131153
```{r gene_coverage, warning=FALSE, fig.height = 8, fig.width = 12, fig.align = "center"}
@@ -313,12 +335,38 @@ graphics::par(opar)
313335
```
314336

315337
#### Add base and amino acid annotation
338+
339+
**Use twill to mark position with SNV**:
316340
```{r base_aa_coverage, warning=FALSE, fig.height = 10, fig.width = 12, fig.align = "center"}
317341
library(ggpattern)
318342
# create plot with twill mark
319-
ggcoverage(data = track.df, color = "grey", range.position = "out", single.nuc=T, rect.color = "white") +
343+
ggcoverage(data = track.df, color = "grey", range.position = "out",
344+
single.nuc=T, rect.color = "white") +
345+
geom_base(bam.file = bam.file,
346+
bs.fa.seq = BSgenome.Hsapiens.UCSC.hg19,
347+
mark.type = "twill") +
348+
geom_ideogram(genome = "hg19",plot.space = 0)
349+
```
350+
351+
**Use star to mark position with SNV**:
352+
```{r base_aa_coverage_star, warning=FALSE, fig.height = 10, fig.width = 12, fig.align = "center"}
353+
# create plot with star mark
354+
ggcoverage(data = track.df, color = "grey", range.position = "out",
355+
single.nuc=T, rect.color = "white") +
356+
geom_base(bam.file = bam.file,
357+
bs.fa.seq = BSgenome.Hsapiens.UCSC.hg19,
358+
mark.type = "star") +
359+
geom_ideogram(genome = "hg19",plot.space = 0)
360+
```
361+
362+
**Highlight position with SNV**:
363+
```{r base_aa_coverage_highlight, warning=FALSE, fig.height = 10, fig.width = 12, fig.align = "center"}
364+
# highlight
365+
ggcoverage(data = track.df, color = "grey", range.position = "out",
366+
single.nuc=T, rect.color = "white") +
320367
geom_base(bam.file = bam.file,
321-
bs.fa.seq = BSgenome.Hsapiens.UCSC.hg19) +
368+
bs.fa.seq = BSgenome.Hsapiens.UCSC.hg19,
369+
mark.type = "highlight") +
322370
geom_ideogram(genome = "hg19",plot.space = 0)
323371
```
324372

0 commit comments

Comments
 (0)