Replicate statistics
Tom
2019-04-15
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, basename, cbind,
## colMeans, colnames, colSums, dirname, 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## Loading required package: BiocParallel##
## Attaching package: 'DelayedArray'## The following objects are masked from 'package:matrixStats':
##
## colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges## The following objects are masked from 'package:base':
##
## aperm, 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] "pADamID-Hap1_S_LMNB2"print(samples)## [1] "pADamID-Hap1_S_r1_LMNB2" "pADamID-Hap1_S_r2_LMNB2"
## [3] "pADamID-Hap1_S_r3_LMNB2"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 10000 <NA>
## 2 chr1 10000 20000 <NA>
## 3 chr1 20000 30000 <NA>
## 4 chr1 30000 40000 <NA>
## 5 chr1 40000 50000 <NA>
## 6 chr1 50000 60000 <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] "pADamID-Hap1_S_r1_LMNB2"
## [1] "31242 / 7180359 GATC bins have reads."
## [1] ""
## [1] "pADamID-Hap1_S_r2_LMNB2"
## [1] "30685 / 7180359 GATC bins have reads."
## [1] ""
## [1] "pADamID-Hap1_S_r3_LMNB2"
## [1] "57871 / 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: 10kb"## Warning in min(which(!fequal(fromdist, 0))): no non-missing arguments to
## min; returning Inf
## <simpleError in npoints[nextix, , drop = FALSE]: incorrect number of dimensions>
## [1] "Current bin size: 20kb"## [1] "Current bin size: 80kb"## Warning in min(which(!fequal(fromdist, 0))): no non-missing arguments to
## min; returning Inf## Warning in .find.triangle.within.face(drawing, faceName): NAs introduced by
## coercion to integer range
## <simpleError in if (nointersect) { fix <- ix break}: missing value where TRUE/FALSE needed>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: 10kb"
## <simpleError in `[.data.frame`(hmm.df, , c("seqnames", "start", "end", names.plot)): undefined columns selected>
## [1] "Current bin size: 20kb"
## <simpleError in `[.data.frame`(hmm.df, , c("seqnames", "start", "end", names.plot)): undefined columns selected>
## [1] "Current bin size: 80kb"
## <simpleError in `[.data.frame`(hmm.df, , c("seqnames", "start", "end", names.plot)): undefined columns selected>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.5.3 (2019-03-11)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.6 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 stats graphics grDevices utils datasets
## [8] methods 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_3.1.0 rtracklayer_1.42.1
## [7] GenomicAlignments_1.18.1 Rsamtools_1.34.1
## [9] Biostrings_2.50.2 XVector_0.22.0
## [11] SummarizedExperiment_1.12.0 DelayedArray_0.8.0
## [13] BiocParallel_1.16.6 matrixStats_0.54.0
## [15] Biobase_2.42.0 GenomicRanges_1.34.0
## [17] GenomeInfoDb_1.18.2 IRanges_2.16.0
## [19] S4Vectors_0.20.1 BiocGenerics_0.28.0
##
## loaded via a namespace (and not attached):
## [1] tidyselect_0.2.5 xfun_0.5 purrr_0.3.0
## [4] lattice_0.20-38 colorspace_1.4-0 htmltools_0.3.6
## [7] yaml_2.2.0 XML_3.98-1.17 RBGL_1.58.1
## [10] rlang_0.3.1 pillar_1.3.1 glue_1.3.0
## [13] withr_2.1.2 GenomeInfoDbData_1.2.0 plyr_1.8.4
## [16] stringr_1.4.0 zlibbioc_1.28.0 munsell_0.5.0
## [19] gtable_0.2.0 evaluate_0.13 labeling_0.3
## [22] knitr_1.21 Rcpp_1.0.0 scales_1.0.0
## [25] graph_1.60.0 digest_0.6.18 stringi_1.3.1
## [28] dplyr_0.8.0.1 grid_3.5.3 tools_3.5.3
## [31] bitops_1.0-6 magrittr_1.5 lazyeval_0.2.1
## [34] RCurl_1.95-4.11 tibble_2.0.1 crayon_1.3.4
## [37] pkgconfig_2.0.2 Matrix_1.2-15 assertthat_0.2.0
## [40] rmarkdown_1.11 reshape_0.8.8 R6_2.4.0
## [43] compiler_3.5.3