Differential Modification Analysis of CD4 ChIP-Seq Dataset

For histone mark H3K4me2

Preliminary Setup

First we load the necessary libraries, along with a set of utility functions.

library(here)
library(stringr)
library(glue)
library(magrittr)
library(openxlsx)
library(SummarizedExperiment)
library(dplyr)
library(edgeR)
library(limma)
library(csaw)
library(sva)
library(ggplot2)
library(scales)
library(GGally)
library(ggalt)
library(ggthemes)
library(splines)
library(reshape2)
library(assertthat)
library(ggfortify)
library(broom)
library(ks)
library(RColorBrewer)

library(BSgenome.Hsapiens.UCSC.hg38)

library(doParallel)
ncores <- getOption("mc.cores", default=parallel::detectCores(logical = FALSE))
options(mc.cores=ncores)
registerDoParallel(cores=ncores)
library(BiocParallel)
register(DoparParam())

options(future.globals.maxSize=4 * 1024^3)
library(future)
plan(multicore)

source(here("scripts", "utilities.R"))

# Required in order to use DGEList objects with future
length.DGEList <- function(x) {
    length(unclass(x))
}

We also set the random seed for reproducibility:

set.seed(1986)

Data Loading and Preprocessing

First we load the consensus peaks called from the reads pooled from all samples. This consensus peak set is not biased toward or against any sample or condition, and therefore the peak significance is expected to be independent of any differential binding in that peak.

peakfile <- here(
    "peak_calls", "epic_hg38.analysisSet",
    glue("{params$histone_mark}_condition.ALL_donor.ALL"),
    "peaks_noBL_IDR.narrowPeak")
allpeaks <- {
    read.narrowPeak(peakfile) %>% as("GRanges") %>%
        assign_into(seqinfo(.), seqinfo(BSgenome.Hsapiens.UCSC.hg38)[seqlevels(.)]) %>%
        setNames(.$name)
}

Now we’ll load the ChIP-seq read count data set from RDS files containing SummarizedExperiment objects, and modify them to use the sample names as column names. We also ensure that the column order is identical between the two objects. Lastly, we filter out any windows with fewer than one count per sample. This is a very mild filtering criterion, but it often eliminates many windows, greatly easing the subsequent computational burden of computing the real filtering threshold.

sexpfile <-
    here("saved_data",
         glue_data(params, "csaw-counts_{window_size}-windows_{fragment_length}-reads_{histone_mark}.RDS"))
full.sexp <- readRDS(sexpfile)
colnames(full.sexp) <- colData(full.sexp)$SampleName
sexp <- full.sexp %>% .[rowSums(assay(.)) >= ncol(.),]
# Exepected number of counts per read, based on overlapping multiple windows.
# NOTE: Assumes windows exactly tile the genome (no overlaps, no gaps).
colData(sexp)$CountDupFactor <- (colData(sexp)$ext - 1) / median(width(rowRanges(sexp))) + 1

We extract the sample metadata from the SummarizedExperiment. We set all factors to use a sum-to-zero variant of the treatment-contrast coding, which will ease the subtraction of batch effects later.

sample.table <- colData(sexp) %>%
    as.data.frame %>% autoFactorize %>%
    mutate(days_after_activation=time_point %>% str_extract("\\d+$") %>% as.numeric(),
           time_point=factor(days_after_activation) %>% `levels<-`(glue("D{levels(.)}")),
           group=interaction(cell_type, time_point, sep="")) %>%
    autoFactorize %>%
    set_rownames(colnames(sexp))
for (i in names(sample.table)) {
    if (is.factor(sample.table[[i]]) && nlevels(sample.table[[i]]) > 1) {
        sample.table[[i]] %<>% C(code_control_named(levels(.)))
    }
}

Peak and Window Filtering

We first select all peaks with an IDR value of 0.2 or less. This is a quite relaxed threshold in order to ensure that even peaks present in only one condition will be included.

idr.threshold <- 0.2
genome.size <- seqlengths(seqinfo(allpeaks)) %>% as.numeric %>% sum
peaks <- allpeaks[allpeaks$qValue >= -log10(idr.threshold)]
pct.covered <- width(peaks) %>% sum %>% divide_by(genome.size) %>% multiply_by(100)
mean.pct.reads <- sexp %>% subsetByOverlaps(peaks) %>% assay("counts") %>%
    colSums %>% divide_by(colData(sexp) %$% {totals * CountDupFactor}) %>% multiply_by(100) %>%
    mean
message(glue("Selected {length(peaks)} peaks at an IDR threshold of {format(idr.threshold, digits=3)}, with an average width of {round(mean(width(peaks)))} nucleotides and covering a total of {format(pct.covered, digits=3)}% of the genome, containing on average {format(mean.pct.reads, digits=3)}% of reads"))

Selected 21631 peaks at an IDR threshold of 0.2, with an average width of 3512 nucleotides and covering a total of 2.46% of the genome, containing on average 16.6% of reads

Then we select only the windows overlapping the selected peaks.

sexp %<>% subsetByOverlaps(peaks)

Now we create a DGEList from the counts.

## Extract gene metadata and colapse lists
all.window.meta <- rowRanges(sexp) %>% as.data.frame %>%
    select(-width, -strand) %>% rename(chr=seqnames)
# Convert list columns to character vectors
all.window.meta[] %<>% lapply(function(x) if (is.list(x)) sapply(x, str_c, collapse=",") else x)
rownames(all.window.meta) <- glue_data(all.window.meta, "{chr}:{start}-{end}")
dge <- asDGEList(sexp) %>%
    assign_into(.$offset, NULL) %>%
    assign_into(.$genes, all.window.meta) %>%
    set_rownames(rownames(all.window.meta))

Then we select all windows with an average read count of at least 5.

count.threshold <- 5
filter.threshold <- aveLogCPM(count.threshold, lib.size=mean(dge$samples$lib.size))
filt <- aveLogCPM(dge) >= filter.threshold
dge <- dge[filt,]
message(glue("Excluding {sum(filt == FALSE)} out of {length(filt)} peak-overlapping windows ({format(100*(1-mean(filt)), digits=3)}%) with average count below {count.threshold}."))

Excluding 46893 out of 170399 peak-overlapping windows (27.5%) with average count below 5.

Model fitting

Now we are ready to fit a model to the data. We begin by building the design matrix, inclulding coefficients for the interaction of cell type and time point (group).

design <- model.matrix(~0 + group, sample.table, strip.prefixes = TRUE)
colnames(design)

[1] "NaiveD0"   "MemoryD0"  "NaiveD1"   "MemoryD1"  "NaiveD5"   "MemoryD5" 
[7] "NaiveD14"  "MemoryD14"

NaiveD0

MemoryD0

NaiveD1

MemoryD1

NaiveD5

MemoryD5

NaiveD14

MemoryD14

Instead of incorporating all the batch effects we think might be important into the design matrix, we let SVA infer the important confounding factors from the data, and add those to the design matrix.

dge %<>% calcNormFactors()
logcpm <- cpm(dge, log=TRUE, prior.count=1)
svobj <- sva(logcpm, design)

Number of significant surrogate variables is:  8 
Iteration (out of 5 ):1  2  3  4  5  

svmat <- add.numbered.colnames(svobj$sv, "SV")
design.sv <- cbind(design, svmat)

We plot the surrogate variables against known variables and also test for significant correlations between them.

sv.cor.tests <- expand.grid(SV=colnames(svmat),
                            Var=c("cell_type", "time_point", "group", "donor_id", "totals")) %>%
    select(Var, SV) %>%
    group_by_all() %>% 
    do(mod=list(lm(as.formula(glue("{.$SV} ~ {.$Var}")), data=cbind(sample.table, svmat)))) %>%
    mutate(aov=lapply(mod, anova)) %>%
    ungroup() %>%
    mutate(PValue=vapply(.$aov, . %$% {`Pr(>F)`[1]}, numeric(1)),
           FDR=p.adjust(PValue, method="BH")) %>%
    select(-mod, -aov)
sv.cor.tests %>% arrange(PValue)

ggduo.dataXY(sample.table[c("cell_type", "time_point", "donor_id", "totals")], svmat,
             extraData=sample.table["group"],
             types=list(comboVertical="dot_no_facet"),
             mapping=aes(color=group, shape=donor_id))

We now proceed to computing normalization factors, estimating dispersions, and fitting the quasi-likelihood GLM.

dge %<>% estimateDisp(design.sv, robust=TRUE)
plotBCV(dge)

fit <- glmQLFit(dge, design.sv, robust=TRUE)
plotQLDisp(fit)

Next, we define the differential expression tests we wish to perform as contrasts. Each contrast is an arithmetic expression in terms of the model coefficients.

celltypes <- levels(sample.table$cell_type)
all.timepoints <- levels(sample.table$time_point)
nonzero.timepoints <- setdiff(all.timepoints, "D0")

timepoint.anova.tests <- setNames(lapply(celltypes, function(ct) {
    setNames(glue("{ct}{nonzero.timepoints} - {ct}D0"),
             nm=glue("{ct}.D0v{nonzero.timepoints}"))
}), nm=str_c(celltypes, ".AllT"))
timepoint.single.tests <- as.list(unlist(unname(timepoint.anova.tests)))
celltype.singlet.tests <-
    as.list(setNames(glue("Memory{all.timepoints} - Naive{all.timepoints}"),
                     glue("NvM.{all.timepoints}")))
celltype.allt.test <- list(NvM.AllT=unlist(celltype.singlet.tests))
factorial.singlet.tests <-
    as.list(setNames(glue("(Memory{nonzero.timepoints} - MemoryD0) - (Naive{nonzero.timepoints} - NaiveD0)"),
                     glue("Fac.{nonzero.timepoints}")))
factorial.allt.test <- list(Fac.AllT=unlist(factorial.singlet.tests))
mi.vs.nf.test <- list(MD0vND14="MemoryD0 - NaiveD14")
sv.test <- list(SV=colnames(svmat))
alltests <- c(timepoint.anova.tests, timepoint.single.tests,
              celltype.allt.test, celltype.singlet.tests,
              factorial.allt.test, factorial.singlet.tests,
              mi.vs.nf.test,
              sv.test
              ) %>% lapply(as_glue)
print(alltests)

$Naive.AllT
NaiveD1 - NaiveD0
NaiveD5 - NaiveD0
NaiveD14 - NaiveD0

$Memory.AllT
MemoryD1 - MemoryD0
MemoryD5 - MemoryD0
MemoryD14 - MemoryD0

$Naive.D0vD1
NaiveD1 - NaiveD0

$Naive.D0vD5
NaiveD5 - NaiveD0

$Naive.D0vD14
NaiveD14 - NaiveD0

$Memory.D0vD1
MemoryD1 - MemoryD0

$Memory.D0vD5
MemoryD5 - MemoryD0

$Memory.D0vD14
MemoryD14 - MemoryD0

$NvM.AllT
MemoryD0 - NaiveD0
MemoryD1 - NaiveD1
MemoryD5 - NaiveD5
MemoryD14 - NaiveD14

$NvM.D0
MemoryD0 - NaiveD0

$NvM.D1
MemoryD1 - NaiveD1

$NvM.D5
MemoryD5 - NaiveD5

$NvM.D14
MemoryD14 - NaiveD14

$Fac.AllT
(MemoryD1 - MemoryD0) - (NaiveD1 - NaiveD0)
(MemoryD5 - MemoryD0) - (NaiveD5 - NaiveD0)
(MemoryD14 - MemoryD0) - (NaiveD14 - NaiveD0)

$Fac.D1
(MemoryD1 - MemoryD0) - (NaiveD1 - NaiveD0)

$Fac.D5
(MemoryD5 - MemoryD0) - (NaiveD5 - NaiveD0)

$Fac.D14
(MemoryD14 - MemoryD0) - (NaiveD14 - NaiveD0)

$MD0vND14
MemoryD0 - NaiveD14

$SV
SV1
SV2
SV3
SV4
SV5
SV6
SV7
SV8

We now perform the differential modification tests for each contrast or set of contrasts. For a single contrast, the test is analogous to a t-test. For a multi-contrast test, the global null hypothesis that all contrasts are equal to zero is used, analogous to an F-test.

window.results.tables <- bplapply(alltests, function(ct) {
    ctmat <- makeContrasts(contrasts = ct, levels=design.sv) %>% set_colnames(names(ct))
    ctest <- glmQLFTest(fit, contrast = ctmat)
    topTags(ctest, n=Inf, sort.by="none") %>% as.data.frame %>% add.qvalue
})

p <- bplapply(names(window.results.tables), function(testname) {
    pvals <- window.results.tables[[testname]]$PValue
    pi0 <- pi0est(pvals)$pi0
    p <- plotpvals(pvals, ptn = pi0) +
        labs(title=glue("P-value histogram for {testname}"),
             subtitle=glue("Est. Non-Null Prop.: {(format((1 - pi0) * 100, digits=3))}%"))
})
ggprint(p)

Next, because we wish to control the FDR at the level of peaks, not individual windows, we must combine the results for all the windows in each peak. The p-value for each peak tests the combined null hypothesis that no windows in that peak are differentially modified with respect to the specified contrast(s).

peak.results.tables <- bplapply(window.results.tables, function(tt) {
    gr <- as(tt, "GRanges")
    olap <- findOverlaps(peaks, gr)
    combineOverlaps(olap, tt) %>%
        filter(!is.na(PValue)) %>%
        arrange(PValue) %>%
        add.qvalue
})
peak.results.tables %<>% lapply(. %>% add.qvalue)

We now visualize the p-value histogram for each test. Ideally, each p-value histogram should either be uniform or biased toward zero. Any other shape is indicative of a failure of some modelling step.

p <- bplapply(names(peak.results.tables), function(testname) {
    pvals <- peak.results.tables[[testname]]$PValue
    pi0 <- pi0est(pvals)$pi0
    p <- plotpvals(pvals, ptn = pi0) +
        labs(title=glue("P-value histogram for {testname}"),
             subtitle=glue("Est. Non-Null Prop.: {(format((1 - pi0) * 100, digits=3))}%"))
})
ggprint(p)

peak.results.tables %>% lapply(function(rt) {
    summarize(rt,
              NumSignif=FDR %>% is_weakly_less_than(0.1) %>% sum,
              EstNonNull=floor(length(PValue) * (1 - propTrueNull(PValue))))
}) %>% 
    do.call(what=rbind) %>% 
    mutate(Test=rownames(.)) %>%
    select(Test, everything())

# Take all genes significant at 10% FDR or top 100 genes, whichever is more
filtered.peak.results.tables <- lapply(peak.results.tables, . %>% filter(FDR <= 0.1 | seq_len(nrow(.)) <= 100))
dir.create(here("saved_data", "ChIP-seq"), recursive = TRUE, showWarnings = FALSE)
dir.create(here("results", "ChIP-seq"), recursive = TRUE, showWarnings = FALSE)
saveRDS(window.results.tables,
        here('saved_data', 'ChIP-seq',
             glue('{params$histone_mark}-window-diffmod.RDS')))
saveRDS(peak.results.tables,
        here('saved_data', 'ChIP-seq',
             glue('{params$histone_mark}-peak-diffmod.RDS')))
save.image(here('saved_data', 'ChIP-seq',
                glue('{params$histone_mark}-diffmod.rda')))
write.xlsx(filtered.peak.results.tables,
           here("results", "ChIP-seq",
                glue('{params$histone_mark}-peak-diffmod.xlsx')))