In stephaniehicks/benchmarkfdrData2019: Data and Benchmarking Results from Korthauer and Kimes et al. (2019)

knitr::opts_chunk$set(echo = TRUE)

Summary

In this case study, we perform differential peak calling on ChIP-seq data of the histone H3K4Me3 for samples from two cell lines (K562 and GM12878) using publicly available data generated by the ENCODE Project. The same data set is used for all ChIP-seq differential testing case studies.

The csaw package provides a powerful and flexible approach to differential peak calling based on counting reads in sliding windows across the genome. The sliding window approach allows for the detection of differential peaks over sharper or wider regions, without being constrained by the original peak boundaries obtained from the peak calling algorithm, e.g. MACS. Differential peaks are called at each window using edgeR and p-values are combined within adjacent windows (clusters).

For this analysis, the default csaw settings are used as described in the csaw "User Guide."

Workspace Setup

library(dplyr)
library(ggplot2)
library(SummarizedBenchmark)
library(BiocParallel)
library(rtracklayer)
library(Rsamtools)
library(Rsubread)
library(csaw)
library(edgeR)
library(GenomicAlignments)

## load helper functions
for (f in list.files("../R", "\\.(r|R)$", full.names = TRUE)) {
    source(f)
}

## project data/results folders
datdir <- "data"
resdir <- "results"
sbdir <- "../../results/ChIPseq"
dir.create(datdir, showWarnings = FALSE)
dir.create(resdir, showWarnings = FALSE)
dir.create(sbdir, showWarnings = FALSE)

## intermediary files we create below
count_file <- file.path(resdir, "h3k4me3-csaw-counts.rds")
filtered_file <- file.path(resdir, "h3k4me3-csaw-counts-filtered.rds")
result_file <- file.path(resdir, "h3k4me3-csaw-results.rds")
bench_file <- file.path(sbdir, "h3k4me3-csaw-benchmark.rds")
bench_file_cov <- file.path(sbdir, "h3k4me3-csaw-cov-benchmark.rds")
bench_file_uninf <- file.path(sbdir, "h3k4me3-csaw-uninf-benchmark.rds")

## set up parallel backend
cores <- as.numeric(Sys.getenv("SLURM_NTASKS"))
multicoreParam <- MulticoreParam(workers = cores)

Data Preparation

We download the fastq files directly from the UCSC ENCODE portal.

broad_url <- "http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHistone/"
broad_fastqs <- c("wgEncodeBroadHistoneGm12878H3k4me3StdRawDataRep1.fastq.gz",
                  "wgEncodeBroadHistoneGm12878H3k04me3StdRawDataRep2V2.fastq.gz",
                  "wgEncodeBroadHistoneK562H3k4me3StdRawDataRep1.fastq.gz",
                  "wgEncodeBroadHistoneK562H3k4me3StdRawDataRep2.fastq.gz")

for (i_fq in broad_fastqs) {
    if (!file.exists(file.path(datdir, i_fq))) {
        download.file(paste0(broad_url, i_fq),
                      destfile = file.path(datdir, i_fq))
    }
}

uw_url <- "http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwHistone/"
uw_fastqs <- c("wgEncodeUwHistoneGm12878H3k4me3StdRawDataRep1.fastq.gz",
               "wgEncodeUwHistoneGm12878H3k4me3StdRawDataRep2.fastq.gz",
               "wgEncodeUwHistoneK562H3k4me3StdRawDataRep1.fastq.gz",
               "wgEncodeUwHistoneK562H3k4me3StdRawDataRep2.fastq.gz")

for (i_fq in uw_fastqs) {
    if (!file.exists(file.path(datdir, i_fq))) {
        download.file(paste0(uw_url, i_fq),
                      destfile = file.path(datdir, i_fq))
    }
}

We use the GRCh38 reference genome used in the ENCODE3 analysis. If the reference index is not already available, we must first download the reference file. Note that this is a different reference than what was originally used for the data in earlier analyses of the ENCODE project.

ref_url <- paste0("https://www.encodeproject.org/files/",
                  "GRCh38_no_alt_analysis_set_GCA_000001405.15/@@download/",
                  "GRCh38_no_alt_analysis_set_GCA_000001405.15.fasta.gz")
ref_fastagz <- file.path(datdir, basename(ref_url))
ref_fasta <- gsub("\\.gz$", "", ref_fastagz)

if (!file.exists(file.path(datdir, "ref_index.00.b.tab"))) {
    if (!file.exists(ref_fasta)) {
        download.file(ref_url, destfile = ref_fastagz)
        system(paste("gunzip", ref_fastagz))
    }
    buildindex(basename = file.path(datdir, "ref_index"),
               reference = ref_fasta)
}

We determine sample metadata from the file names.

fqfiles <- file.path(datdir, c(broad_fastqs, uw_fastqs))
labs <- gsub("wgEncode(.*)Histone.*", "\\1", basename(fqfiles))
cells <- gsub("wgEncode.*Histone(.*)H3k.*", "\\1", basename(fqfiles))

meta <- data.frame(cellline = cells, lab = labs,
                   fqfile = fqfiles,
                   bamfile = gsub("\\.fastq.gz", ".bam", fqfiles),
                   sortedfile = gsub("\\.fastq.gz", ".sorted.bam", fqfiles),
                   stringsAsFactors = FALSE)

We also download blacklisted regions for GRCh38 from ENCODE (https://www.encodeproject.org/annotations/ENCSR636HFF/).

blacklist_url <- paste0("http://mitra.stanford.edu/kundaje/akundaje/",
                        "release/blacklists/hg38-human/hg38.blacklist.bed.gz")
if (!file.exists(file.path(datdir, basename(blacklist_url)))) {
    download.file(blacklist_url, destfile = file.path(datdir, basename(blacklist_url)))
}

bl <- import(file.path(datdir, basename(blacklist_url)))

Standard set of parameters are defined for the analysis. Only the canonical set of chromosomes are used in the analysis.

std_chr <- paste0("chr", c(1:22, "X", "Y"))
param <- readParam(minq = 20, discard = bl, restrict = std_chr)

Read Counting

We count reads within sliding windows across the genome.

if (file.exists(count_file)) {
    win_cnts <- readRDS(count_file)
} else {
    ## Broad samples have quality score w/ phred offset +33 (standard)
    b_meta <- meta[meta$lab == "Broad", ]
    b_unaligned <- !file.exists(b_meta$bamfile)
    if (any(b_unaligned)) {
        align(index = file.path(datdir, "ref_index"),
              readfile1 = b_meta$fqfile[b_unaligned],
              TH1 = 2, type = 1, phredOffset = 33,
              output_file = b_meta$bamfile[b_unaligned])
    }

    ## UW samples have quality score w/ phred offset +64 (old Illumina scale)
    u_meta <- meta[meta$lab == "Uw", ]
    u_unaligned <- !file.exists(u_meta$bamfile)
    if (any(u_unaligned)) {
        align(index = file.path(datdir, "ref_index"),
              readfile1 = u_meta$fqfile[u_unaligned],
              TH1 = 2, type = 1, phredOffset = 64,
              output_file = u_meta$bamfile[u_unaligned])
    }

    ## sort bam files
    for (i in 1:nrow(meta)) {
        if (!file.exists(meta$sortedfile[i])) {
            sortBam(meta$bamfile[i], gsub("\\.bam$", "", meta$sortedfile[i]))
        }
    }

    ## mark duplicates w/ picard and index bam files
    if (any(!file.exists(paste0(meta$sortedfile, ".bai")))) {
        temp_bam <- file.path(datdir, "temp_dups.bam")
        temp_file <- file.path(datdir, "temp_mets.txt")
        temp_dir <- file.path(datdir, "temp_dups")
        dir.create(temp_dir)
        for (i_bam in meta$sortedfile) {
            code <- paste0("java -jar ${PICARD_TOOLS_HOME}/picard.jar ",
                           "MarkDuplicates I=%s O=%s M=%s ", 
                           "TMP_DIR=%s AS=true REMOVE_DUPLICATES=false ",
                           "VALIDATION_STRINGENCY=SILENT")
            code <- sprintf(code, i_bam, temp_bam, temp_file, temp_dir)
            code <- system(code)
            stopifnot(code == 0L)
            file.rename(temp_bam, i_bam)
        }
        unlink(temp_file)
        unlink(temp_dir, recursive = TRUE)

        indexBam(meta$sortedfile)
    }

    ## use correlateReads to determine fragment length (remove dups)
    b_x <- correlateReads(b_meta$sortedfile, param = reform(param, dedup = TRUE))
    u_x <- correlateReads(u_meta$sortedfile, param = reform(param, dedup = TRUE))
    b_frag_len <- which.max(b_x) - 1
    u_frag_len <- which.max(u_x) - 1

    ## count reads in sliding windows (keep dups)
    win_cnts <- windowCounts(meta$sortedfile, param = param, width = 150,
                             ext = list(rep(c(b_frag_len, u_frag_len), each = 4), 150))

    ## save unfiltered counts
    saveRDS(win_cnts, file = count_file)
}

We can apply prefiltering on windows with low abundance as described in Lun and Smyth (2016).

if (file.exists(filtered_file)) {
    filtered_cnts <- readRDS(filtered_file)
} else {
    ## filter windows by abundance
    bins <- windowCounts(meta$sortedfile, bin = TRUE, width = 2000, param = param)
    filter_stat <- filterWindows(win_cnts, bins, type = "global")
    keep <- filter_stat$filter > log2(2)
    filtered_cnts <- win_cnts[keep, ]

    ## save filtered counts
    saveRDS(filtered_cnts, file = filtered_file)
}

Data Analysis

Differential Testing

We use edgeR to test for differential binding with the filtered data.

if (file.exists(result_file)) {
    res_ranges <- readRDS(result_file)
} else {
    ## normalize for library-specific trended biases
    filtered_cnts <- normOffsets(filtered_cnts, type = "loess")

    ## create model design
    design <- model.matrix(~ lab + cellline,
                           data = as.data.frame(meta))

    ## run quasi-likelihood F-test
    y <- asDGEList(filtered_cnts)
    y <- estimateDisp(y, design)
    fit <- glmQLFit(y, design, robust = TRUE)
    res <- glmQLFTest(fit)

    ## merge p-values across regions
    merged <- mergeWindows(rowRanges(filtered_cnts), tol = 100, max.width = 5000)
    tab_comb <- combineTests(merged$id, res$table)
    tab_best <- getBestTest(merged$id, res$table)

    ## save results
    res_ranges <- merged$region
    elementMetadata(res_ranges) <-
        data.frame(tab_comb,
                   best_pos = mid(ranges(rowRanges(filtered_cnts[tab_best$best]))),
                   best_logFC = tab_best$logFC)

    ## get overall mean counts for each window
    merged_cnts <- summarizeOverlaps(res_ranges, BamFileList(meta$sortedfile))
    res_ranges$meancnt <- rowMeans(assays(merged_cnts)$counts)

    saveRDS(res_ranges, file = result_file)
}

## covert to df for downstream analysis
res_df <- as.data.frame(res_ranges)

## add random covariate 
set.seed(7778)
res_df$uninf_covar = rnorm(nrow(res_df))

Covariate Diagnostics

Number of Windows

rank_scatter(res_df, pvalue = "PValue", covariate = "nWindows") +
    ggtitle("Number of windows as independent covariate") +
    xlab("Number of Windows")

strat_hist(res_df, pvalue = "PValue", covariate = "nWindows", maxy = 35)

Region Width

rank_scatter(res_df, pvalue = "PValue", covariate = "width") +
    ggtitle("Region width as independent covariate") +
    xlab("Region Width")

strat_hist(res_df, pvalue = "PValue", covariate = "width", maxy = 35)

Mean Coverage

rank_scatter(res_df, pvalue = "PValue", covariate = "meancnt") +
    ggtitle("Mean coverage as independent covariate") +
    xlab("Mean coverage")

strat_hist(res_df, pvalue = "PValue", covariate = "meancnt", maxy = 35)

Random

rank_scatter(res_df, pvalue = "PValue", covariate = "uninf_covar") +
    ggtitle("Random independent covariate") +
    xlab("Region Width")

strat_hist(res_df, pvalue = "PValue", covariate = "uninf_covar", maxy = 35)

Multiple-Testing Correction

We use the common BenchDesign with the set of multiple testing correction methods already included. We investigate both the width of the regions and their mean coverage as the independent covariate. We won't assess the number of windows, since this is very tightly correlated with the width of the region.

cor(res_df[,c("width", "nWindows", "meancnt")])

First, we'll use the region width covariate.

if (!file.exists(bench_file)) {
    res_df$ind_covariate <- res_df$width
    res_df$pval <- res_df$PValue

    bd <- initializeBenchDesign()
    sb <- buildBench(bd, data = res_df, ftCols = "ind_covariate")
    saveRDS(sb, file = bench_file)
} else {
    sb <- readRDS(bench_file)
}

Next, we'll use the mean coverage covariate.

if (!file.exists(bench_file_cov)) {
    res_df$ind_covariate <- res_df$meancnt
    res_df$pval <- res_df$PValue

    bd <- initializeBenchDesign()
    sbC <- buildBench(bd, data = res_df, ftCols = "ind_covariate")
    saveRDS(sbC, file = bench_file_cov)
} else {
    sbC <- readRDS(bench_file_cov)
}

We'll also compare to the random covariate.

if (!file.exists(bench_file_uninf)) {
    res_df$ind_covariate <- res_df$uninf_covar
    res_df$pval <- res_df$PValue

    bd <- initializeBenchDesign()
    sbU <- buildBench(bd, data = res_df, ftCols = "ind_covariate")
    saveRDS(sbU, file = bench_file_uninf)
} else {
    sbU <- readRDS(bench_file_uninf)
}

Benchmark Metrics

Region width

Next, we'll add the default performance metric for q-value assays. First, we have to rename the assay to 'qvalue'.

assayNames(sb) <- "qvalue"
sb <- addDefaultMetrics(sb)

Now, we'll plot the results.

rejections_scatter(sb, supplementary = FALSE)
rejection_scatter_bins(sb, covariate = "ind_covariate",
                       bins = 4, supplementary = FALSE)

plotFDRMethodsOverlap(sb, alpha = 0.05, nsets = ncol(sb),
                      order.by = "freq", decreasing = TRUE,
                      supplementary = FALSE)

covariateLinePlot(sb, alpha = 0.05, covname = "ind_covariate")

Mean coverage

We'll do the same for the mean coverage covariate.

assayNames(sbC) <- "qvalue"
sbC <- addDefaultMetrics(sbC)

Now, we'll plot the results.

rejections_scatter(sbC, supplementary = FALSE)
rejection_scatter_bins(sbC, covariate = "ind_covariate",
                       bins = 4, supplementary = FALSE)

plotFDRMethodsOverlap(sbC, alpha = 0.05, nsets = ncol(sbC),
                      order.by = "freq", decreasing = TRUE,
                      supplementary = FALSE)

covariateLinePlot(sbC, alpha = 0.05, covname = "ind_covariate")

Random

We'll do the same for the random (uninformative covariate).

assayNames(sbU) <- "qvalue"
sbU <- addDefaultMetrics(sbU)

Now, we'll plot the results.

rejections_scatter(sbU, supplementary = FALSE)
rejection_scatter_bins(sbU, covariate = "ind_covariate",
                       bins = 4, supplementary = FALSE)

plotFDRMethodsOverlap(sbU, alpha = 0.05, nsets = ncol(sbU),
                      order.by = "freq", decreasing = TRUE,
                      supplementary = FALSE)

covariateLinePlot(sbU, alpha = 0.05, covname = "ind_covariate")

Covariate comparison

Here we compare the method ranks for the two covariates at alpha = 0.10.

plotMethodRanks(c(bench_file, bench_file_cov, bench_file_uninf), 
                colLabels = c("Region width", "Mean Coverage", "Random"), 
                alpha = 0.10, xlab = "Covariate", 
                excludeMethods = NULL)

Session Info

sessionInfo()

stephaniehicks/benchmarkfdrData2019 documentation built on June 20, 2021, 10 a.m.