library(testthat)
library(PeakSegPipeline)
library(data.table)
context("demo")
download.to <- function
(u, f, writeFun=if(grepl("bigWig", f))writeBin else writeLines){
if(!file.exists(f)){
require(httr)
f.dir <- dirname(f)
dir.create(f.dir, showWarnings=FALSE, recursive=TRUE)
request <- GET(u)
stop_for_status(request)
writeFun(content(request), f)
}
}
## Download bigWig files from github.
bigWig.part.vec <- c(
"Input/MS010302",
"bcell/MS010302",
"Input/MS002202",
"kidney/MS002202",
"Input/MS026601",
"bcell/MS026601",
"Input/MS002201",
"kidney/MS002201"
)
label.txt <- "
chr10:33,061,897-33,162,814 noPeaks
chr10:33,456,000-33,484,755 peakStart kidney
chr10:33,597,317-33,635,209 peakEnd kidney
chr10:33,662,034-33,974,942 noPeaks
chr10:35,182,820-35,261,001 noPeaks
chr10:35,261,418-35,314,654 peakStart bcell kidney
chr10:35,343,031-35,398,459 peakEnd bcell kidney
chr10:38,041,023-38,102,554 noPeaks
chr10:38,296,008-38,307,179 peakStart bcell kidney
chr10:38,379,045-38,391,967 peakStart bcell kidney
chr10:38,404,899-38,412,089 peakEnd bcell kidney
chr10:38,413,073-38,444,133 noPeaks
chr10:38,585,584-38,643,190 noPeaks
chr10:38,643,191-38,650,766 peakStart bcell kidney
chr10:38,731,066-38,750,574 peakEnd bcell kidney
chr10:38,750,960-38,790,663 noPeaks
chr10:38,807,475-38,815,200 noPeaks
chr10:38,815,201-38,816,355 peakStart bcell kidney Input
chr10:38,818,377-38,819,342 peakEnd bcell kidney Input
chr10:39,098,319-39,111,384 noPeaks
chr10:39,125,134-39,125,550 peakStart bcell kidney Input
chr10:39,125,594-39,126,266 peakEnd bcell kidney Input
chr10:39,126,866-39,140,858 noPeaks
"
test.data.dir <- file.path(Sys.getenv("HOME"), "PeakSegPipeline-test")
##test.data.dir <- file.path(tempdir(), "PeakSegPipeline-test")
non.integer.dir <- file.path(test.data.dir, "non-integer")
demo.dir <- file.path(test.data.dir, "demo")
chrom.sizes.file <- tempfile()
chrom.sizes <- data.table(chrom="chr10", bases=128616069)
fwrite(chrom.sizes, chrom.sizes.file, sep="\t", col.names=FALSE)
repos.url <- "https://raw.githubusercontent.com/tdhock/input-test-data/master/"
for(bigWig.part in bigWig.part.vec){
suffix <- ifelse(grepl("MS026601|MS002201", bigWig.part), "/", "_/")
bigWig.file <- file.path(
non.integer.dir, "samples",
sub("/", suffix, bigWig.part), "coverage.bigWig")
bigWig.url <- paste0(repos.url, bigWig.part, ".bigwig")
download.to(bigWig.url, bigWig.file)
demo.bigWig <- sub("non-integer", "demo", bigWig.file)
if(!file.exists(demo.bigWig)){
dir.create(dirname(demo.bigWig), showWarnings=FALSE, recursive=TRUE)
bw.dt <- readBigWig(bigWig.file, "chr10", 0, 128616069)
out.dt <- data.table(chrom="chr10", bw.dt)
demo.bedGraph <- sub("bigWig", "bedGraph", demo.bigWig)
fwrite(out.dt, demo.bedGraph, sep="\t", col.names=FALSE)
bedGraphToBigWig(demo.bedGraph, chrom.sizes.file, demo.bigWig)
unlink(demo.bedGraph)
}
}
for(set.dir in c(demo.dir)){
labels.file <- file.path(set.dir, "labels", "some_labels.txt")
dir.create(dirname(labels.file), showWarnings=FALSE, recursive=TRUE)
writeLines(label.txt, labels.file)
problems.bed <- file.path(set.dir, "problems.bed")
unlink(problems.bed)
cat("chr10 60000 17974675
chr10 18024675 38818835
chr10 38868835 39154935
chr10 42746000 46426964
chr10 47529169 47792476
chr10 47892476 48055707
chr10 48105707 49095536
chr10 49195536 51137410
chr10 51187410 51398845
chr10 51448845 125869472
chr10 125919472 128616069
", file=problems.bed)
}
## Pipeline should run to completion for integer count data.
system(paste("bigWigToBedGraph", demo.bigWig, "/dev/stdout|head"))
index.html <- file.path(demo.dir, "index.html")
unlink(index.html)
jobs_create_run(demo.dir)
test_that("index.html is created", {
expect_true(file.exists(index.html))
})
if(interactive()){
## Post-processing to explain the output.
joint.problems.dt <- fread(paste("cat", file.path(
demo.dir, "problems", "*", "jointProblems.bed")))
setnames(joint.problems.dt, c("chrom", "problemStart", "problemEnd"))
peaks.glob <- file.path(
demo.dir, "problems", "*", "jointProblems", "*", "peaks.bed")
##Sys.glob(peaks.glob)
joint.peaks.dt <- fread(paste("cat", peaks.glob))
setnames(
joint.peaks.dt,
c("chrom", "peakStart", "peakEnd", "sample.path", "mean"))
joint.peaks.dt[, sample.id := sub(".*/", "", sample.path)]
joint.peaks.dt[, sample.group := sub("/.*", "", sample.path)]
chunks.dt <- fread(paste("cat", file.path(
demo.dir, "problems", "*", "chunks", "*", "chunk.bed")))
setnames(chunks.dt, c("chrom", "chunkStart", "chunkEnd"))
all.problems <- fread(file.path(demo.dir, "problems.bed"))
setnames(all.problems, c("chrom", "problemStart", "problemEnd"))
data.start <- min(chunks.dt$chunkStart)
data.end <- max(chunks.dt$chunkEnd)
two.problems <- all.problems[!(
problemEnd < data.start |
data.end < problemStart),]
both.problems <- rbind(
data.table(two.problems, y="separate problems"),
data.table(joint.problems.dt, y="joint problems"))
limits.list <- list(
c(34, 35),
c(35, 36),
c(36, 37),
c(37, 38),
c(38, 39))
limits.dt <- data.table(matrix(unlist(limits.list)*1e6, ncol=2, byrow=TRUE))
setnames(limits.dt, c("plotStart", "plotEnd"))
limits.dt$y <- "plots in un-labeled regions"
labels.bed.vec <- Sys.glob(file.path(
demo.dir, "samples", "*", "*", "labels.bed"))
all.labels.list <- list()
for(labels.bed in labels.bed.vec){
sample.dir <- dirname(labels.bed)
sample.id <- basename(sample.dir)
group.dir <- dirname(sample.dir)
sample.group <- basename(group.dir)
sample.labels <- fread(labels.bed)
setnames(sample.labels, c("chrom", "labelStart", "labelEnd", "annotation"))
all.labels.list[[labels.bed]] <- data.table(
sample.id, sample.group, sample.labels)
}
all.labels <- do.call(rbind, all.labels.list)
input.labels <- all.labels[sample.group=="Input", list(
prop.noPeaks=mean(annotation=="noPeaks")
), by=.(labelStart, labelEnd)]
setkey(input.labels, labelStart, labelEnd)
input.pred <- joint.peaks.dt[, list(
n.Input=sum(sample.group=="Input")
), by=.(peakStart, peakEnd)]
setkey(input.pred, peakStart, peakEnd)
labeled.input <- foverlaps(input.pred, input.labels, nomatch=0L)
thresh.dt <- labeled.input[, data.table(WeightedROC(
n.Input, ifelse(prop.noPeaks==0, 1, -1)))]
thresh.best <- thresh.dt[which.min(FP+FN),]
## threshold is smallest n.Input that is classified as non-specific.
setkey(joint.peaks.dt, peakStart, peakEnd)
peaks.with.counts <- input.pred[joint.peaks.dt]
peaks.with.counts[, specificity := ifelse(
n.Input >= thresh.best$threshold, "non-specific", "specific")]
gg <- ggplot()+
coord_cartesian(xlim=c(data.start, data.end)/1e3, expand=TRUE)+
geom_segment(aes(
problemStart/1e3, y,
xend=problemEnd/1e3, yend=y),
color="blue",
data=both.problems)+
geom_point(aes(
problemStart/1e3, y),
color="blue",
data=both.problems)+
geom_tallrect(aes(
xmin=chunkStart/1e3, xmax=chunkEnd/1e3),
alpha=0.1,
data=chunks.dt)+
geom_segment(aes(
plotStart/1e3, y,
xend=plotEnd/1e3, yend=y),
data=limits.dt)+
geom_point(aes(
plotStart/1e3, y),
data=limits.dt)+
geom_point(aes(
peakStart/1e3, sample.path, color=specificity),
data=peaks.with.counts)+
geom_segment(aes(
peakStart/1e3, sample.path,
xend=peakEnd/1e3, yend=sample.path, color=specificity),
data=peaks.with.counts)+
ylab("")+
xlab("position on chr10 (kb = kilo bases)")+
scale_color_manual(values=c("non-specific"="red", specific="deepskyblue"))
png(file.path(demo.dir, "figure-demo-overview.png"), res=100, h=200, w=1000)
print(gg)
dev.off()
coverage.bigWig.vec <- Sys.glob(file.path(
demo.dir, "samples", "*", "*", "coverage.bigWig"))
ann.colors <-
c(noPeaks="#f6f4bf",
peakStart="#ffafaf",
peakEnd="#ff4c4c",
peaks="#a445ee")
img.list <- list()
for(limit.i in seq_along(limits.list)){
limit.vec <- limits.list[[limit.i]]*1e6
print(limit.vec)
coverage.list <- list()
for(coverage.bigWig in coverage.bigWig.vec){
sample.dir <- dirname(coverage.bigWig)
sample.id <- basename(sample.dir)
group.dir <- dirname(sample.dir)
sample.group <- basename(group.dir)
cmd <- sprintf(
"bigWigToBedGraph -chrom=chr10 -start=%d -end=%d %s /dev/stdout",
limit.vec[1], limit.vec[2],
coverage.bigWig)
sample.coverage <- fread(cmd)
setnames(sample.coverage, c("chrom", "chromStart", "chromEnd", "count"))
coverage.list[[coverage.bigWig]] <- data.table(
sample.id, sample.group, sample.coverage)
}
coverage <- do.call(rbind, coverage.list)
show.peaks <- joint.peaks.dt[
!(peakEnd < limit.vec[1] | limit.vec[2] < peakStart),]
show.labels <- all.labels[
!(labelEnd < limit.vec[1] | limit.vec[2] < labelStart),]
gg <- ggplot()+
theme_bw()+
theme(panel.margin=grid::unit(0, "lines"))+
facet_grid(sample.group + sample.id ~ ., scales="free")+
geom_tallrect(aes(
xmin=labelStart/1e3,
xmax=labelEnd/1e3,
fill=annotation),
alpha=0.5,
data=show.labels)+
scale_fill_manual(values=ann.colors)+
scale_y_continuous(
"aligned read coverage",
breaks=function(limits){
lim <- floor(limits[2])
if(lim==0){
Inf
}else{
lim
}
})+
scale_x_continuous(paste(
"position on",
coverage$chrom[1],
"(kb = kilo bases)"))+
geom_step(aes(
chromStart/1e3, count),
color="grey50",
data=coverage)+
coord_cartesian(xlim=limit.vec/1e3, expand=FALSE)
if(nrow(show.peaks)){
gg <- gg+
geom_point(aes(
peakStart/1e3, 0),
color="deepskyblue",
size=3,
data=show.peaks)+
geom_segment(aes(
peakStart/1e3, 0,
xend=peakEnd/1e3, yend=0),
color="deepskyblue",
size=3,
data=show.peaks)
}
f <- sprintf(
"%s/figure-demo-%d-%d.png",
demo.dir,
limit.vec[1]/1e6, limit.vec[2]/1e6)
png(f, res=100, width=1000, height=600)
print(gg)
dev.off()
img.list[[limit.i]] <- sprintf(
'<p><img src="%s" /></p>',
basename(f))
}
index.lines <- readLines(index.html)
new.index.lines <- c('
<h2>Explanation of results below</h2>
<ul>
<li>
PeakSegFPOP + PeakSegJoint were trained using four labeled samples
(two H3K36me3 samples, two Input samples). These labeled samples
have colored rectangular labels/annotations in the plots below,
which indicate genomic regions where a scientist has visually
determined presence or absence of peaks.
</li>
<li>
The black peak predictions represent the PeakSegFPOP model, which
predicts peaks separately for each sample. These peaks do not
occur in the exact same positions across samples, but are clustered
together to define "joint problems" where the PeakSegJoint model
is fit to all samples.
</li>
<li>
The blue peak predictions represent the PeakSegJoint model, which
predicts peaks separately for each genomic region, but jointly
across all samples. These peaks occur in the exact same positions
across samples, and are reported in the final output files
(<a href="peaks_summary.tsv">peaks_summary.tsv</a>,
<a href="peaks_matrix.tsv">peaks_matrix.tsv</a>).
</li>
<li>
The models were used to predict peaks in eight samples,
four of which have no
labels. The un-labeled samples appear without any colored
rectangles in the plots below, and can be used for evaluating the
accuracy of the learned model. For example, consider the genomic region below:</li>
<img src="problems/chr10:18024675-38818835/chunks/chr10:33061897-33974942/figure-predictions.png" />
<li>On the left, the model predicted a common peak in all bcell and kidney samples
(even though this peak was not labeled in any of the samples).</li>
<li>In the middle, the model predicted a common peak in the two kidney samples,
and no corresponding peak in the two bcell samples
(only one of the two bcell samples was labeled as having noPeaks,
and only one of the two kidney samples was labeled as having a peakStart and peakEnd).</li>
</ul>
', index.lines, '
<h2>More plots</h2>
<p>The plots below are usually not included in the output,
but we added them here, so you can see peak predictions
at all genomic regions that have coverage in this data set.</p>
<img src="figure-demo-overview.png" />
', unlist(img.list))
writeLines(new.index.lines, index.html)
browseURL(index.html)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.