Replicate statistics
Tom
2018-10-02
Overview
This document provides quality control figures for DamID replicates.
Set-up
# Libraries used
library(GenomicRanges)## Loading required package: stats4## Loading required package: BiocGenerics## Loading required package: parallel##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, cbind, colMeans,
## colnames, colSums, do.call, duplicated, eval, evalq, Filter,
## Find, get, grep, grepl, intersect, is.unsorted, lapply,
## lengths, Map, mapply, match, mget, order, paste, pmax,
## pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce,
## rowMeans, rownames, rowSums, sapply, setdiff, sort, table,
## tapply, union, unique, unsplit, which, which.max, which.min## Loading required package: S4Vectors##
## Attaching package: 'S4Vectors'## The following object is masked from 'package:base':
##
## expand.grid## Loading required package: IRanges## Loading required package: GenomeInfoDblibrary(GenomicAlignments)## Loading required package: SummarizedExperiment## Loading required package: Biobase## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.## Loading required package: DelayedArray## Loading required package: matrixStats##
## Attaching package: 'matrixStats'## The following objects are masked from 'package:Biobase':
##
## anyMissing, rowMedians##
## Attaching package: 'DelayedArray'## The following objects are masked from 'package:matrixStats':
##
## colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges## The following object is masked from 'package:base':
##
## apply## Loading required package: Biostrings## Loading required package: XVector##
## Attaching package: 'Biostrings'## The following object is masked from 'package:DelayedArray':
##
## type## The following object is masked from 'package:base':
##
## strsplit## Loading required package: Rsamtoolslibrary(rtracklayer)
library(ggplot2)
library(reshape2)
library(GGally)
library(RColorBrewer)
# Basenames
basename <- snakemake@params[["basename"]]
samples <- snakemake@config[["replicates"]][[basename]]
if (length(samples) < 2) {
stop("Not enough samples!")
}
# Bin sizes
bins <- snakemake@config[["bins"]]
bins <- strsplit(bins, ",")[[1]]
bins.kb <- paste0(bins, "kb")
bins.with_gatc <- c("gatc", paste0(bins, "kb"))
bins.count <- length(bins.with_gatc)
# Input directories - note that the .Rmd file is the working directory!
map_dir <- file.path(snakemake@config[["report_relative_location"]],
snakemake@config[["output_dir"]],
snakemake@config["out_map"])
counts_dir <- file.path(snakemake@config[["report_relative_location"]],
snakemake@config[["output_dir"]],
snakemake@config[["out_count"]])
norm_dir <- file.path(snakemake@config[["report_relative_location"]],
snakemake@config[["output_dir"]],
snakemake@config[["out_norm"]])
hmm_dir <- file.path(snakemake@config[["report_relative_location"]],
snakemake@config[["output_dir"]],
snakemake@config[["out_hmm"]])Samples to be analyzed:
print(basename)## [1] "H1_4xAP3"print(samples)## [1] "H1_r1_4xAP3" "H1_r2_4xAP3"Read-in data
I will read in the data
- Normalized scores
- HMM model
FixGATCsequences <- function(df) {
# The GATC fragments have the downside that they overlap, as they end with
# GATC and also start with GATC. This function quickly fixes this issue.
df.gr <- as(df, "GRanges")
# Fix the start / end
start(df.gr) <- start(df.gr) + 2
end(df.gr) <- end(df.gr) - 2
# Also trim the object to be within range
df.gr <- trim(df.gr)
as(df.gr, "data.frame")[, c("seqnames", "start", "end", "score")]
}
# 1) Normalized values
replicate.norm <- list()
for (s in samples) {
target.norm <- list()
for (bin in bins.with_gatc) {
df.name <- file.path(norm_dir,
paste0("bin-", bin),
paste0(s,
"-",
bin,
".norm.txt.gz"))
df <- read.table(df.name,
sep = "\t", stringsAsFactors = FALSE,
col.names = c("seqnames", "start", "end", "score"))
if (bin == "gatc") {
df <- FixGATCsequences(df)
}
target.norm <- c(target.norm, list(df))
}
names(target.norm) <- bins.with_gatc
replicate.norm <- c(replicate.norm, list(target.norm))
}
names(replicate.norm) <- samples
# 1b) Normalized combined values
replicate.norm.combined <- list()
for (bin in bins.with_gatc) {
df.name <- file.path(norm_dir,
paste0("bin-", bin),
paste0(basename,
"-",
bin,
"-combined.norm.txt.gz"))
df <- read.table(df.name,
sep = "\t", stringsAsFactors = FALSE,
col.names = c("seqnames", "start", "end", "score"))
if (bin == "gatc") {
df <- FixGATCsequences(df)
}
replicate.norm.combined <- c(replicate.norm.combined, list(df))
}
names(replicate.norm.combined) <- bins.with_gatc
# 2) HMM
replicate.hmm <- list()
for (s in samples) {
target.hmm <- list()
for (bin in bins.kb) {
df.name <- file.path(hmm_dir,
paste0("bin-", bin),
paste0(s,
"-",
bin,
"_HMM.txt.gz"))
df <- read.table(df.name,
sep = "\t", stringsAsFactors = FALSE,
col.names = c("seqnames", "start", "end", "score"))
if (bin == "gatc") {
df <- FixGATCsequences(df)
}
target.hmm <- c(target.hmm, list(df))
}
names(target.hmm) <- bins.kb
replicate.hmm <- c(replicate.hmm, list(target.hmm))
}
names(replicate.hmm) <- samplesTop of the normalized counts and HMM model for sample #1:
head(replicate.norm[[samples[1]]][[1]])## seqnames start end score
## 1 chr1 2 11161 NA
## 2 chr1 11161 12412 NA
## 3 chr1 12412 12462 NA
## 4 chr1 12462 12687 NA
## 5 chr1 12687 12830 NA
## 6 chr1 12830 13316 NAhead(replicate.hmm[[samples[1]]][[1]])## seqnames start end score
## 1 chr1 0 1000 <NA>
## 2 chr1 1000 2000 <NA>
## 3 chr1 2000 3000 <NA>
## 4 chr1 3000 4000 <NA>
## 5 chr1 4000 5000 <NA>
## 6 chr1 5000 6000 <NA>Quality controls
Various quality plots will be shown below. This list is still under construction.
1) Bins called
First, how many bins have reads in replicate #1 and replicate #2?
for (i in 1:length(samples)) {
print(samples[i])
print(paste(sum(! is.na(replicate.norm[[samples[i]]][[1]][, 4])), # NA bins
"/",
nrow(replicate.norm[[samples[i]]][[1]]),
"GATC bins have reads.",
sep = " "))
print("")
}## [1] "H1_r1_4xAP3"
## [1] "1228764 / 7180359 GATC bins have reads."
## [1] ""
## [1] "H1_r2_4xAP3"
## [1] "1019712 / 7180359 GATC bins have reads."
## [1] ""Or the same thing in a plot
library(Vennerable)
plotGATCOverlap <- function(replicate.norm, samples, n_bin = 1) {
x <- ! is.na(replicate.norm[[samples[1]]][[n_bin]][, 4])
y <- ! is.na(replicate.norm[[samples[2]]][[n_bin]][, 4])
ab <- sum(x & y)
a <- sum(x & !y)
b <- sum(!x & y)
if (length(samples) == 2) {
w <- Venn(SetNames = samples)
Weights(w) <- c(0, a, b, ab)
plot(w, show = list(Faces = F))
} else if (length(samples >= 3)) {
z <- ! is.na(replicate.norm[[samples[3]]][[n_bin]][, 4])
abc <- sum(x & y & z)
ac <- sum(x & !y & z)
bc <- sum(!x & y & z)
c <- sum(!x & !y & z)
if (length(samples) == 3) {
w <- Venn(SetNames = samples)
Weights(w) <- c(0, a, b, ab, c, ac, bc, abc)
plot(w, show = list(Faces = F))
} else if (length(samples) == 4) {
q <- ! is.na(replicate.norm[[samples[4]]][[n_bin]][, 4])
abcd <- sum(x & y & z & q)
abd <- sum(x & y & !z & q)
acd <- sum(x & !y & z & q)
bcd <- sum(!x & y & z & q)
ad <- sum(x & !y & !z & q)
bd <- sum(!x & y & !z & q)
cd <- sum(!x & !y & z & q)
d <- sum(!x & !y & !z & q)
w <- Venn(SetNames = samples)
Weights(w) <- c(0, a, b, ab, c, ac, bc, abc, d, ad, bd, abd, cd, acd, bcd, abcd)
plot(w, show = list(Faces = F))
} else {
stop("Sorry, this number of names is not supported yet")
}
}
}
for (i in 1:length(bins.with_gatc)) {
print(paste0("Current bin size: ", bins.with_gatc[i]))
tryCatch(plotGATCOverlap(replicate.norm, samples, n_bin = i),
error = function(e) print(e))
}## [1] "Current bin size: gatc"
## [1] "Current bin size: 1kb"
## [1] "Current bin size: 2kb"
## [1] "Current bin size: 5kb"
## [1] "Current bin size: 10kb"
## [1] "Current bin size: 20kb"
## [1] "Current bin size: 25kb"
## [1] "Current bin size: 50kb"
## [1] "Current bin size: 80kb"
## [1] "Current bin size: 100kb"
## [1] "Current bin size: 250kb"
2) HMM overlap
How well does the HMM calling overlap between the replicates for different bin sizes?
library(Vennerable)
plotVennFromList <- function(df.list) {
names.plot <- paste0("bin_", names(df.list))
# Convert HMM list into one GRanges object
# For now, this assumes that:
# 1) Model go from small to large
hmm.gr <- as(df.list[[1]],
"GRanges")
names(mcols(hmm.gr)) <- names.plot[1]
for (i in 2:length(df.list)) {
# Determine overlap
ovl <- findOverlaps(hmm.gr,
as(df.list[[i]],
"GRanges"),
type = "within")
mcols(hmm.gr)[, names.plot[i]] <- df.list[[i]][subjectHits(ovl), "score"]
}
# Convert into data frame
hmm.df <- as(hmm.gr,
"data.frame")
# Only select the interesting data, and only the complete cases
df.tmp <- hmm.df[, c("seqnames", "start", "end", names.plot)]
df.tmp <- df.tmp[complete.cases(df.tmp), ]
# Next, create a vector of overlap
if (length(names.plot) == 2) {
x <- (df.tmp[, names.plot[1]] == "AD")
y <- (df.tmp[, names.plot[2]] == "AD")
ab <- sum(x & y)
a <- sum(x & !y)
b <- sum(!x & y)
w <- Venn(SetNames = names.plot)
Weights(w) <- c(0, a, b, ab)
plot(w, show = list(Faces = F))
} else if (length(names.plot) == 3) {
x <- (df.tmp[, names.plot[1]] == "AD")
y <- (df.tmp[, names.plot[2]] == "AD")
z <- (df.tmp[, names.plot[3]] == "AD")
abc <- sum(x & y & z)
ab <- sum(x & y & !z)
ac <- sum(x & !y & z)
bc <- sum(!x & y & z)
a <- sum(x & !y & !z)
b <- sum(!x & y & !z)
c <- sum(!x & !y & z)
w <- Venn(SetNames = names.plot)
Weights(w) <- c(0, a, b, ab, c, ac, bc, abc)
plot(w, show = list(Faces = F))
} else if (length(names.plot) == 4) {
x <- (df.tmp[, names.plot[1]] == "AD")
y <- (df.tmp[, names.plot[2]] == "AD")
z <- (df.tmp[, names.plot[3]] == "AD")
q <- (df.tmp[, names.plot[4]] == "AD")
abcd <- sum(x & y & z & q)
abc <- sum(x & y & z & !q)
abd <- sum(x & y & !z & q)
acd <- sum(x & !y & z & q)
bcd <- sum(!x & y & z & q)
ab <- sum(x & y & !z & !q)
ac <- sum(x & !y & z & !q)
ad <- sum(x & !y & !z & q)
bc <- sum(!x & y & z & !q)
bd <- sum(!x & y & !z & q)
cd <- sum(!x & !y & z & q)
a <- sum(x & !y & !z & !q)
b <- sum(!x & y & !z & !q)
c <- sum(!x & !y & z & !q)
d <- sum(!x & !y & !z & q)
w <- Venn(SetNames = names.plot)
Weights(w) <- c(0, a, b, ab, c, ac, bc, abc, d, ad, bd, abd, cd, acd, bcd, abcd)
plot(w, show = list(Faces = F))
} else {
stop("Sorry, this number of names is not supported yet")
}
}
# Unfortunately, there is a bug in Vennerable. Try it, but don't bother if
# it doesn't work.
for (b in bins.kb) {
print(paste0("Current bin size: ", b))
hmm.list <- lapply(samples, function(x) replicate.hmm[[x]][[b]])
names(hmm.list) <- samples
tryCatch(plotVennFromList(hmm.list),
error = function(e) print(e))
}## [1] "Current bin size: 1kb"
## <simpleError in xy[1:3, ]: subscript out of bounds>
## [1] "Current bin size: 2kb"
## [1] "Current bin size: 5kb"
## [1] "Current bin size: 10kb"
## [1] "Current bin size: 20kb"
## [1] "Current bin size: 25kb"
## [1] "Current bin size: 50kb"
## [1] "Current bin size: 80kb"
## [1] "Current bin size: 100kb"
## [1] "Current bin size: 250kb"
3) Correlation between samples
First, the correlation + spearman number.
op <- par(pty="s")
n <- 10000
for (i in 1:length(bins.with_gatc)) {
# Get a data frame with the observations
df <- do.call(cbind, lapply(samples, function(x) replicate.norm[[x]][[i]][, 4]))
df <- df[complete.cases(df), ]
df <- data.frame(df)
names(df) <- samples
# Take a samples of this data frame
s <- sample(1:nrow(df), min(n, nrow(df)), replace = F)
df.reduced <- df[s, ]
# Get the limits
limits <- c(min(df), max(df))
# Plot
my_dens <- function(data, mapping, ...) {
ggplot(data = data, mapping=mapping) +
geom_point(..., alpha = 0.2) +
geom_abline(slope = 1, lty = "dashed", col = "red")
}
print(ggpairs(df.reduced,
lower = list(continuous = my_dens)) +
ggtitle(paste(bins.with_gatc[i],
"|", nrow(df.reduced), "points plotted")) +
xlab("Normalized score (log2)") +
ylab("Normalized score (log2)") +
theme_bw())
# plot(df,
# pch = 19, col = alpha("black", 0.2),
# xlab = paste(samples[1], "(log 2)"),
# ylab = paste(samples[2], "(log 2)"),
# xlim = limits, ylim = limits,
# main = paste(bins.with_gatc[i],
# "|", nrow(df), "points plotted"))
# abline(a = 0, b = 1, lty = 2, col = "red")
}
par(op)Also, let’s plot the cross correlation. Note, this is only for replicate #1 and #2.
# op <- par(pty="s")
#
# for (i in 1:length(bins.with_gatc)) {
# x <- replicate.norm[[samples[1]]][[i]][, 4]
# y <- replicate.norm[[samples[2]]][[i]][, 4]
#
# ccf(x, y, na.action = na.pass,
# lag.max = 20,
# main = paste(basename,
# bins.with_gatc[i],
# sep = " - "))
# }
#
# par(op)4) Autocorrelation
Multiple replicates should filter out the (technical) noise present in samples. If we use the autocorrelation as quality measure, can we see this?
# Set-up for the ACF
acf.df <- c()
for (b in bins.with_gatc) {
acf.vec <- c()
for (s in samples) {
# Where -b is the bin size and -s is the sample
# Calculate the ACF - with the notion that NAs are simply skipped
a <- acf(replicate.norm[[s]][[b]][, 4], lag.max = 2, na.action = na.pass, plot = F)
acf.vec <- c(acf.vec, a$acf[2])
}
# And for the combined replicate
a <- acf(replicate.norm.combined[[b]][, 4], lag.max = 2, na.action = na.pass, plot = F)
acf.vec <- c(acf.vec, a$acf[2])
acf.df <- cbind(acf.df, acf.vec)
}
acf.df <- data.frame(acf.df)
names(acf.df) <- bins.with_gatc
acf.df$sample <- c(samples, paste0(basename, "-combined"))
# Plot all of this
acf.df.melt <- melt(acf.df, id.vars = "sample")
ggplot(acf.df.melt, aes(x = variable, y = value, col = sample)) +
geom_point() +
ggtitle("Combined replicate ACF") +
xlab("Bin size") +
ylab("ACF") +
scale_color_brewer(palette = "Set1") +
theme_bw()
SessionInfo
sessionInfo()## R version 3.4.4 (2018-03-15)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.5 LTS
##
## Matrix products: default
## BLAS: /usr/lib/libblas/libblas.so.3.6.0
## LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 methods stats graphics grDevices utils
## [8] datasets base
##
## other attached packages:
## [1] Vennerable_3.1.0.9000 RColorBrewer_1.1-2
## [3] GGally_1.4.0 reshape2_1.4.3
## [5] ggplot2_2.2.1 rtracklayer_1.38.3
## [7] GenomicAlignments_1.14.1 Rsamtools_1.30.0
## [9] Biostrings_2.46.0 XVector_0.18.0
## [11] SummarizedExperiment_1.8.1 DelayedArray_0.4.1
## [13] matrixStats_0.53.1 Biobase_2.38.0
## [15] GenomicRanges_1.30.1 GenomeInfoDb_1.14.0
## [17] IRanges_2.12.0 S4Vectors_0.16.0
## [19] BiocGenerics_0.24.0
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.14 pillar_1.1.0 compiler_3.4.4
## [4] plyr_1.8.4 bitops_1.0-6 tools_3.4.4
## [7] zlibbioc_1.24.0 digest_0.6.15 tibble_1.4.1
## [10] evaluate_0.10.1 gtable_0.2.0 lattice_0.20-35
## [13] rlang_0.1.6 graph_1.56.0 Matrix_1.2-14
## [16] yaml_2.1.19 GenomeInfoDbData_1.0.0 stringr_1.3.0
## [19] knitr_1.18 rprojroot_1.3-2 grid_3.4.4
## [22] reshape_0.8.7 RBGL_1.54.0 XML_3.98-1.11
## [25] BiocParallel_1.12.0 rmarkdown_1.8 magrittr_1.5
## [28] backports_1.1.2 scales_0.5.0 htmltools_0.3.6
## [31] colorspace_1.3-2 labeling_0.3 stringi_1.1.6
## [34] lazyeval_0.2.1 RCurl_1.95-4.10 munsell_0.4.3