Introduction

Users want to provide here background information about the design of their ChIP-Seq project.

Background and objectives

This report describes the analysis of several ChIP-Seq experiments studying the DNA binding patterns of the transcriptions factors … from organism ….

Experimental design

Typically, users want to specify here all information relevant for the analysis of their NGS study. This includes detailed descriptions of FASTQ files, experimental design, reference genome, gene annotations, etc.

Workflow environment

Load packages

The systemPipeR package needs to be loaded to perform the analysis steps shown in this report (H Backman and Girke 2016). The package allows users to run the entire analysis workflow interactively or with a single command while also generating the corresponding analysis report. For details see systemPipeR's main vignette.

library(systemPipeR)

Generate workflow environment

systemPipeRdata package is a helper package to generate a fully populated systemPipeR workflow environment in the current working directory with a single command. All the instruction for generating the workflow template are provide in the systemPipeRdata vignette here.

After building and loading the workflow environment generated by genWorkenvir from systemPipeRdata all data inputs are stored in a data/ directory and all analysis results will be written to a separate results/ directory, while the systemPipeChIPseq.Rmd script and the targets file are expected to be located in the parent directory. The R session is expected to run from this parent directory. Additional parameter files are stored under param/.

To work with real data, users want to organize their own data similarly and substitute all test data for their own data. To rerun an established workflow on new data, the initial targets file along with the corresponding FASTQ files are usually the only inputs the user needs to provide.

For more details, please consult the documentation here. More information about the targets files from systemPipeR can be found here.

Run workflow

Now open the R markdown script systemPipeRIBOseq.Rmdin your R IDE (e.g. vim-r or RStudio) and run the workflow as outlined below.

Here pair-end workflow example is provided. Please refer to the main vignette systemPipeR.Rmd for running the workflow with single-end data.

If you are running on a single machine, use following code as an example to check if some tools used in this workflow are in your environment PATH. No warning message should be shown if all tools are installed.

Read preprocessing

Experiment definition provided by targets file

The targets file defines all FASTQ files and sample comparisons of the analysis workflow.

targetspath <- system.file("extdata", "targetsPE_chip.txt", package = "systemPipeR")
targets <- read.delim(targetspath, comment.char = "#")
targets[1:4, -c(5, 6)]
##                     FileName1                   FileName2
## 1 ./data/SRR446027_1.fastq.gz ./data/SRR446027_2.fastq.gz
## 2 ./data/SRR446028_1.fastq.gz ./data/SRR446028_2.fastq.gz
## 3 ./data/SRR446029_1.fastq.gz ./data/SRR446029_2.fastq.gz
## 4 ./data/SRR446030_1.fastq.gz ./data/SRR446030_2.fastq.gz
##   SampleName Factor        Date SampleReference
## 1        M1A     M1 23-Mar-2012                
## 2        M1B     M1 23-Mar-2012                
## 3        A1A     A1 23-Mar-2012             M1A
## 4        A1B     A1 23-Mar-2012             M1B

Read quality filtering and trimming

The following example shows how one can design a custom read preprocessing function using utilities provided by the ShortRead package, and then apply it with preprocessReads in batch mode to all FASTQ samples referenced in the corresponding SYSargs2 instance (trim object below). More detailed information on read preprocessing is provided in systemPipeR's main vignette.

First, we construct SYSargs2 object from cwl and yml param and targets files.

dir_path <- system.file("extdata/cwl/preprocessReads/trim-pe",
    package = "systemPipeR")
trim <- loadWF(targets = targetspath, wf_file = "trim-pe.cwl",
    input_file = "trim-pe.yml", dir_path = dir_path)
trim <- renderWF(trim, inputvars = c(FileName1 = "_FASTQ_PATH1_",
    FileName2 = "_FASTQ_PATH2_", SampleName = "_SampleName_"))
trim
output(trim)[1:2]

Next, we execute the code for trimming all the raw data.

filterFct <- function(fq, cutoff = 20, Nexceptions = 0) {
    qcount <- rowSums(as(quality(fq), "matrix") <= cutoff, na.rm = TRUE)
    fq[qcount <= Nexceptions]
    # Retains reads where Phred scores are >= cutoff with N
    # exceptions
}
preprocessReads(args = trim, Fct = "filterFct(fq, cutoff=20, Nexceptions=0)",
    batchsize = 1e+05)
writeTargetsout(x = trim, file = "targets_chip_trimPE.txt", step = 1,
    new_col = c("FileName1", "FileName2"), new_col_output_index = c(1,
        2), overwrite = TRUE)

FASTQ quality report

The following seeFastq and seeFastqPlot functions generate and plot a series of useful quality statistics for a set of FASTQ files including per cycle quality box plots, base proportions, base-level quality trends, relative k-mer diversity, length and occurrence distribution of reads, number of reads above quality cutoffs and mean quality distribution. The results are written to a PDF file named fastqReport.pdf. Parallelization of FASTQ quality report via scheduler (e.g. Slurm) across several compute nodes.

library(BiocParallel)
library(batchtools)
f <- function(x) {
    library(systemPipeR)
    targets <- system.file("extdata", "targetsPE_chip.txt", package = "systemPipeR")
    dir_path <- system.file("extdata/cwl/preprocessReads/trim-pe",
        package = "systemPipeR")
    trim <- loadWorkflow(targets = targets, wf_file = "trim-pe.cwl",
        input_file = "trim-pe.yml", dir_path = dir_path)
    trim <- renderWF(trim, inputvars = c(FileName1 = "_FASTQ_PATH1_",
        FileName2 = "_FASTQ_PATH2_", SampleName = "_SampleName_"))
    seeFastq(fastq = infile1(trim)[x], batchsize = 1e+05, klength = 8)
}

resources <- list(walltime = 120, ntasks = 1, ncpus = 4, memory = 1024)
param <- BatchtoolsParam(workers = 4, cluster = "slurm", template = "batchtools.slurm.tmpl",
    resources = resources)
fqlist <- bplapply(seq(along = trim), f, BPPARAM = param)

pdf("./results/fastqReport.pdf", height = 18, width = 4 * length(fqlist))
seeFastqPlot(unlist(fqlist, recursive = FALSE))
dev.off()
Figure 1: FASTQ quality report for 18 samples

Alignments

Read mapping with Bowtie2

The NGS reads of this project will be aligned with Bowtie2 against the reference genome sequence (Langmead and Salzberg 2012). The parameter settings of the aligner are defined in the bowtie2-index.cwl and bowtie2-index.yml files. In ChIP-Seq experiments it is usually more appropriate to eliminate reads mapping to multiple locations. To achieve this, users want to remove the argument setting -k 50 non-deterministic in the configuration files.

Building the index:

dir_path <- system.file("extdata/cwl/bowtie2/bowtie2-idx", package = "systemPipeR")
idx <- loadWorkflow(targets = NULL, wf_file = "bowtie2-index.cwl",
    input_file = "bowtie2-index.yml", dir_path = dir_path)
idx <- renderWF(idx)
idx
cmdlist(idx)

## Run in single machine
runCommandline(idx, make_bam = FALSE)

The following submits 18 alignment jobs via a scheduler to a computer cluster.

targets <- system.file("extdata", "targetsPE_chip.txt", package = "systemPipeR")
dir_path <- system.file("extdata/cwl/bowtie2/bowtie2-pe", package = "systemPipeR")
args <- loadWF(targets = targets, wf_file = "bowtie2-mapping-pe.cwl",
    input_file = "bowtie2-mapping-pe.yml", dir_path = dir_path)
args <- renderWF(args, inputvars = c(FileName1 = "_FASTQ_PATH1_",
    FileName2 = "_FASTQ_PATH2_", SampleName = "_SampleName_"))
args
cmdlist(args)[1:2]
output(args)[1:2]
moduleload(modules(args))  # Skip if a module system is not used
resources <- list(walltime = 120, ntasks = 1, ncpus = 4, memory = 1024)
reg <- clusterRun(args, FUN = runCommandline, more.args = list(args = args,
    dir = FALSE), conffile = ".batchtools.conf.R", template = "batchtools.slurm.tmpl",
    Njobs = 18, runid = "01", resourceList = resources)
getStatus(reg = reg)
waitForJobs(reg = reg)
args <- output_update(args, dir = FALSE, replace = TRUE, extension = c(".sam",
    ".bam"))  ## Updates the output(args) to the right location in the subfolders
output(args)

Alternatively, one can run the alignments sequentially on a single system.

args <- runCommandline(args, force = F)

Check whether all BAM files have been created and write out the new targets file.

writeTargetsout(x = args, file = "targets_bam.txt", step = 1,
    new_col = "FileName", new_col_output_index = 1, overwrite = TRUE)
outpaths <- subsetWF(args, slot = "output", subset = 1, index = 1)
file.exists(outpaths)

Read and alignment stats

The following provides an overview of the number of reads in each sample and how many of them aligned to the reference.

read_statsDF <- alignStats(args = args)
write.table(read_statsDF, "results/alignStats.xls", row.names = FALSE,
    quote = FALSE, sep = "\t")
read.delim("results/alignStats.xls")

Utilities for coverage data

The following introduces several utilities useful for ChIP-Seq data. They are not part of the actual workflow.

Rle object stores coverage information

library(rtracklayer)
library(GenomicRanges)
library(Rsamtools)
library(GenomicAlignments)
outpaths <- subsetWF(args, slot = "output", subset = 1, index = 1)
aligns <- readGAlignments(outpaths[1])
cov <- coverage(aligns)
cov

Resizing aligned reads

trim(resize(as(aligns, "GRanges"), width = 200))

Naive peak calling

islands <- slice(cov, lower = 15)
islands[[1]]

Plot coverage for defined region

library(ggbio)
myloc <- c("Chr1", 1, 1e+05)
ga <- readGAlignments(outpaths[1], use.names = TRUE, param = ScanBamParam(which = GRanges(myloc[1],
    IRanges(as.numeric(myloc[2]), as.numeric(myloc[3])))))
autoplot(ga, aes(color = strand, fill = strand), facets = strand ~
    seqnames, stat = "coverage")

Peak calling with MACS2

Merge BAM files of replicates prior to peak calling

Merging BAM files of technical and/or biological replicates can improve the sensitivity of the peak calling by increasing the depth of read coverage. The mergeBamByFactor function merges BAM files based on grouping information specified by a factor, here the Factor column of the imported targets file. It also returns an updated SYSargs2 object containing the paths to the merged BAM files as well as to any unmerged files without replicates. This step can be skipped if merging of BAM files is not desired.

dir_path <- system.file("extdata/cwl/mergeBamByFactor", package = "systemPipeR")
args <- loadWF(targets = "targets_bam.txt", wf_file = "merge-bam.cwl",
    input_file = "merge-bam.yml", dir_path = dir_path)
args <- renderWF(args, inputvars = c(FileName = "_BAM_PATH_",
    SampleName = "_SampleName_"))

args_merge <- mergeBamByFactor(args = args, overwrite = TRUE)
writeTargetsout(x = args_merge, file = "targets_mergeBamByFactor.txt",
    step = 1, new_col = "FileName", new_col_output_index = 1,
    overwrite = TRUE)
# Skip if a module system is not used
module("list")
module("unload", "miniconda2")
module("load", "python/2.7.14")  # Make sure to set up your enviroment variable for MACS2 

Peak calling without input/reference sample

MACS2 can perform peak calling on ChIP-Seq data with and without input samples (Zhang et al. 2008). The following performs peak calling without input on all samples specified in the corresponding args object. Note, due to the small size of the sample data, MACS2 needs to be run here with the nomodel setting. For real data sets, users want to remove this parameter in the corresponding *.param file(s).

dir_path <- system.file("extdata/cwl/MACS2/MACS2-noinput/", package = "systemPipeR")
args <- loadWF(targets = "targets_mergeBamByFactor.txt", wf_file = "macs2.cwl",
    input_file = "macs2.yml", dir_path = dir_path)
args <- renderWF(args, inputvars = c(FileName = "_FASTQ_PATH1_",
    SampleName = "_SampleName_"))

runCommandline(args, make_bam = FALSE, force = T)
outpaths <- subsetWF(args, slot = "output", subset = 1, index = 1)
file.exists(outpaths)
writeTargetsout(x = args, file = "targets_macs.txt", step = 1,
    new_col = "FileName", new_col_output_index = 1, overwrite = TRUE)

Peak calling with input/reference sample

To perform peak calling with input samples, they can be most conveniently specified in the SampleReference column of the initial targets file. The writeTargetsRef function uses this information to create a targets file intermediate for running MACS2 with the corresponding input samples.

writeTargetsRef(infile = "targets_mergeBamByFactor.txt", outfile = "targets_bam_ref.txt",
    silent = FALSE, overwrite = TRUE)
dir_path <- system.file("extdata/cwl/MACS2/MACS2-input/", package = "systemPipeR")
args_input <- loadWF(targets = "targets_bam_ref.txt", wf_file = "macs2-input.cwl",
    input_file = "macs2.yml", dir_path = dir_path)
args_input <- renderWF(args_input, inputvars = c(FileName1 = "_FASTQ_PATH1_",
    FileName2 = "_FASTQ_PATH2_", SampleName = "_SampleName_"))
cmdlist(args_input)[1]
## Run
args_input <- runCommandline(args_input, make_bam = FALSE, force = T)
outpaths_input <- subsetWF(args_input, slot = "output", subset = 1,
    index = 1)
file.exists(outpaths_input)
writeTargetsout(x = args_input, file = "targets_macs_input.txt",
    step = 1, new_col = "FileName", new_col_output_index = 1,
    overwrite = TRUE)

The peak calling results from MACS2 are written for each sample to separate files in the results directory. They are named after the corresponding files with extensions used by MACS2.

Identify consensus peaks

The following example shows how one can identify consensus preaks among two peak sets sharing either a minimum absolute overlap and/or minimum relative overlap using the subsetByOverlaps or olRanges functions, respectively. Note, the latter is a custom function imported below by sourcing it.

# source('http://faculty.ucr.edu/~tgirke/Documents/R_BioCond/My_R_Scripts/rangeoverlapper.R')
outpaths <- subsetWF(args, slot = "output", subset = 1, index = 1)  ## escolher um dos outputs index
peak_M1A <- outpaths["M1A"]
peak_M1A <- as(read.delim(peak_M1A, comment = "#")[, 1:3], "GRanges")
peak_A1A <- outpaths["A1A"]
peak_A1A <- as(read.delim(peak_A1A, comment = "#")[, 1:3], "GRanges")
(myol1 <- subsetByOverlaps(peak_M1A, peak_A1A, minoverlap = 1))
# Returns any overlap
myol2 <- olRanges(query = peak_M1A, subject = peak_A1A, output = "gr")
# Returns any overlap with OL length information
myol2[values(myol2)["OLpercQ"][, 1] >= 50]
# Returns only query peaks with a minimum overlap of 50%

Annotate peaks with genomic context

Annotation with ChIPpeakAnno package

The following annotates the identified peaks with genomic context information using the ChIPpeakAnno and ChIPseeker packages, respectively (Zhu et al. 2010; Yu, Wang, and He 2015).

library(ChIPpeakAnno)
library(GenomicFeatures)
dir_path <- system.file("extdata/cwl/annotate_peaks", package = "systemPipeR")
args <- loadWF(targets = "targets_macs.txt", wf_file = "annotate-peaks.cwl",
    input_file = "annotate-peaks.yml", dir_path = dir_path)
args <- renderWF(args, inputvars = c(FileName = "_FASTQ_PATH1_",
    SampleName = "_SampleName_"))

txdb <- makeTxDbFromGFF(file = "data/tair10.gff", format = "gff",
    dataSource = "TAIR", organism = "Arabidopsis thaliana")
ge <- genes(txdb, columns = c("tx_name", "gene_id", "tx_type"))
for (i in seq(along = args)) {
    peaksGR <- as(read.delim(infile1(args)[i], comment = "#"),
        "GRanges")
    annotatedPeak <- annotatePeakInBatch(peaksGR, AnnotationData = genes(txdb))
    df <- data.frame(as.data.frame(annotatedPeak), as.data.frame(values(ge[values(annotatedPeak)$feature,
        ])))
    outpaths <- subsetWF(args, slot = "output", subset = 1, index = 1)
    write.table(df, outpaths[i], quote = FALSE, row.names = FALSE,
        sep = "\t")
}
writeTargetsout(x = args, file = "targets_peakanno.txt", step = 1,
    new_col = "FileName", new_col_output_index = 1, overwrite = TRUE)

The peak annotation results are written for each peak set to separate files in the results directory. They are named after the corresponding peak files with extensions specified in the annotate_peaks.param file, here *.peaks.annotated.xls.

Annotation with ChIPseeker package

Same as in previous step but using the ChIPseeker package for annotating the peaks.

library(ChIPseeker)
for (i in seq(along = args)) {
    peakAnno <- annotatePeak(infile1(args)[i], TxDb = txdb, verbose = FALSE)
    df <- as.data.frame(peakAnno)
    outpaths <- subsetWF(args, slot = "output", subset = 1, index = 1)
    write.table(df, outpaths[i], quote = FALSE, row.names = FALSE,
        sep = "\t")
}
writeTargetsout(x = args, file = "targets_peakanno.txt", step = 1,
    new_col = "FileName", new_col_output_index = 1, overwrite = TRUE)

Summary plots provided by the ChIPseeker package. Here applied only to one sample for demonstration purposes.

peak <- readPeakFile(infile1(args)[1])
covplot(peak, weightCol = "X.log10.pvalue.")
outpaths <- subsetWF(args, slot = "output", subset = 1, index = 1)
peakHeatmap(outpaths[1], TxDb = txdb, upstream = 1000, downstream = 1000,
    color = "red")
plotAvgProf2(outpaths[1], TxDb = txdb, upstream = 1000, downstream = 1000,
    xlab = "Genomic Region (5'->3')", ylab = "Read Count Frequency")

Count reads overlapping peaks

The countRangeset function is a convenience wrapper to perform read counting iteratively over serveral range sets, here peak range sets. Internally, the read counting is performed with the summarizeOverlaps function from the GenomicAlignments package. The resulting count tables are directly saved to files, one for each peak set.

library(GenomicRanges)
dir_path <- system.file("extdata/cwl/count_rangesets", package = "systemPipeR")
args <- loadWF(targets = "targets_macs.txt", wf_file = "count_rangesets.cwl",
    input_file = "count_rangesets.yml", dir_path = dir_path)
args <- renderWF(args, inputvars = c(FileName = "_FASTQ_PATH1_",
    SampleName = "_SampleName_"))

## Bam Files
targets <- system.file("extdata", "targetsPE_chip.txt", package = "systemPipeR")
dir_path <- system.file("extdata/cwl/bowtie2/bowtie2-pe", package = "systemPipeR")
args_bam <- loadWF(targets = targets, wf_file = "bowtie2-mapping-pe.cwl",
    input_file = "bowtie2-mapping-pe.yml", dir_path = dir_path)
args_bam <- renderWF(args_bam, inputvars = c(FileName = "_FASTQ_PATH1_",
    SampleName = "_SampleName_"))
args_bam <- output_update(args_bam, dir = FALSE, replace = TRUE,
    extension = c(".sam", ".bam"))
outpaths <- subsetWF(args_bam, slot = "output", subset = 1, index = 1)

bfl <- BamFileList(outpaths, yieldSize = 50000, index = character())
countDFnames <- countRangeset(bfl, args, mode = "Union", ignore.strand = TRUE)
writeTargetsout(x = args, file = "targets_countDF.txt", step = 1,
    new_col = "FileName", new_col_output_index = 1, overwrite = TRUE)

Differential binding analysis

The runDiff function performs differential binding analysis in batch mode for several count tables using edgeR or DESeq2 (Robinson, McCarthy, and Smyth 2010; Love, Huber, and Anders 2014). Internally, it calls the functions run_edgeR and run_DESeq2. It also returns the filtering results and plots from the downstream filterDEGs function using the fold change and FDR cutoffs provided under the dbrfilter argument.

dir_path <- system.file("extdata/cwl/rundiff", package = "systemPipeR")
args_diff <- loadWF(targets = "targets_countDF.txt", wf_file = "rundiff.cwl",
    input_file = "rundiff.yml", dir_path = dir_path)
args_diff <- renderWF(args_diff, inputvars = c(FileName = "_FASTQ_PATH1_",
    SampleName = "_SampleName_"))

cmp <- readComp(file = args_bam, format = "matrix")
dbrlist <- runDiff(args = args_diff, diffFct = run_edgeR, targets = targets.as.df(targets(args_bam)),
    cmp = cmp[[1]], independent = TRUE, dbrfilter = c(Fold = 2,
        FDR = 1))
writeTargetsout(x = args_diff, file = "targets_rundiff.txt",
    step = 1, new_col = "FileName", new_col_output_index = 1,
    overwrite = TRUE)

GO term enrichment analysis

The following performs GO term enrichment analysis for each annotated peak set.

dir_path <- system.file("extdata/cwl/annotate_peaks", package = "systemPipeR")
args <- loadWF(targets = "targets_bam_ref.txt", wf_file = "annotate-peaks.cwl",
    input_file = "annotate-peaks.yml", dir_path = dir_path)
args <- renderWF(args, inputvars = c(FileName1 = "_FASTQ_PATH1_",
    FileName2 = "_FASTQ_PATH2_", SampleName = "_SampleName_"))

args_anno <- loadWF(targets = "targets_macs.txt", wf_file = "annotate-peaks.cwl",
    input_file = "annotate-peaks.yml", dir_path = dir_path)
args_anno <- renderWF(args_anno, inputvars = c(FileName = "_FASTQ_PATH1_",
    SampleName = "_SampleName_"))
annofiles <- subsetWF(args_anno, slot = "output", subset = 1,
    index = 1)
gene_ids <- sapply(names(annofiles), function(x) unique(as.character(read.delim(annofiles[x])[,
    "geneId"])), simplify = FALSE)
load("data/GO/catdb.RData")
BatchResult <- GOCluster_Report(catdb = catdb, setlist = gene_ids,
    method = "all", id_type = "gene", CLSZ = 2, cutoff = 0.9,
    gocats = c("MF", "BP", "CC"), recordSpecGO = NULL)

Motif analysis

Parse DNA sequences of peak regions from genome

Enrichment analysis of known DNA binding motifs or de novo discovery of novel motifs requires the DNA sequences of the identified peak regions. To parse the corresponding sequences from the reference genome, the getSeq function from the Biostrings package can be used. The following example parses the sequences for each peak set and saves the results to separate FASTA files, one for each peak set. In addition, the sequences in the FASTA files are ranked (sorted) by increasing p-values as expected by some motif discovery tools, such as BCRANK.

library(Biostrings)
library(seqLogo)
library(BCRANK)
dir_path <- system.file("extdata/cwl/annotate_peaks", package = "systemPipeR")
args <- loadWF(targets = "targets_macs.txt", wf_file = "annotate-peaks.cwl",
    input_file = "annotate-peaks.yml", dir_path = dir_path)
args <- renderWF(args, inputvars = c(FileName = "_FASTQ_PATH1_",
    SampleName = "_SampleName_"))

rangefiles <- infile1(args)
for (i in seq(along = rangefiles)) {
    df <- read.delim(rangefiles[i], comment = "#")
    peaks <- as(df, "GRanges")
    names(peaks) <- paste0(as.character(seqnames(peaks)), "_",
        start(peaks), "-", end(peaks))
    peaks <- peaks[order(values(peaks)$X.log10.pvalue., decreasing = TRUE)]
    pseq <- getSeq(FaFile("./data/tair10.fasta"), peaks)
    names(pseq) <- names(peaks)
    writeXStringSet(pseq, paste0(rangefiles[i], ".fasta"))
}

Motif discovery with BCRANK

The Bioconductor package BCRANK is one of the many tools available for de novo discovery of DNA binding motifs in peak regions of ChIP-Seq experiments. The given example applies this method on the first peak sample set and plots the sequence logo of the highest ranking motif.

set.seed(0)
BCRANKout <- bcrank(paste0(rangefiles[1], ".fasta"), restarts = 25,
    use.P1 = TRUE, use.P2 = TRUE)
toptable(BCRANKout)
topMotif <- toptable(BCRANKout, 1)
weightMatrix <- pwm(topMotif, normalize = FALSE)
weightMatrixNormalized <- pwm(topMotif, normalize = TRUE)
pdf("results/seqlogo.pdf")
seqLogo(weightMatrixNormalized)
dev.off()
Figure 2: One of the motifs identified by BCRANK

Version Information

## R version 4.0.0 (2020-04-24)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/local/bin/R/lib/R/lib/libRblas.so
## LAPACK: /usr/local/bin/R/lib/R/lib/libRlapack.so
## 
## 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] stats4    parallel  stats     graphics  grDevices
## [6] utils     datasets  methods   base     
## 
## other attached packages:
##  [1] systemPipeR_1.22.0          ShortRead_1.46.0           
##  [3] GenomicAlignments_1.24.0    SummarizedExperiment_1.18.2
##  [5] DelayedArray_0.14.1         matrixStats_0.56.0         
##  [7] Biobase_2.48.0              BiocParallel_1.22.0        
##  [9] Rsamtools_2.4.0             Biostrings_2.56.0          
## [11] XVector_0.28.0              GenomicRanges_1.40.0       
## [13] GenomeInfoDb_1.24.2         IRanges_2.22.2             
## [15] S4Vectors_0.26.1            BiocGenerics_0.34.0        
## [17] BiocStyle_2.16.0           
## 
## loaded via a namespace (and not attached):
##   [1] colorspace_1.4-1         rjson_0.2.20            
##   [3] hwriter_1.3.2            ellipsis_0.3.1          
##   [5] rprojroot_1.3-2          fs_1.5.0                
##   [7] rstudioapi_0.11          bit64_4.0.2             
##   [9] AnnotationDbi_1.50.3     codetools_0.2-16        
##  [11] splines_4.0.0            knitr_1.29              
##  [13] jsonlite_1.7.0           annotate_1.66.0         
##  [15] GO.db_3.11.4             dbplyr_1.4.4            
##  [17] png_0.1-7                pheatmap_1.0.12         
##  [19] graph_1.66.0             BiocManager_1.30.10     
##  [21] compiler_4.0.0           httr_1.4.2              
##  [23] GOstats_2.54.0           backports_1.1.9         
##  [25] assertthat_0.2.1         Matrix_1.2-18           
##  [27] limma_3.44.3             formatR_1.7             
##  [29] htmltools_0.5.0          prettyunits_1.1.1       
##  [31] tools_4.0.0              gtable_0.3.0            
##  [33] glue_1.4.1               GenomeInfoDbData_1.2.3  
##  [35] Category_2.54.0          dplyr_1.0.2             
##  [37] rsvg_2.1                 batchtools_0.9.13       
##  [39] rappdirs_0.3.1           V8_3.2.0                
##  [41] Rcpp_1.0.5               pkgdown_1.5.1           
##  [43] vctrs_0.3.2              rtracklayer_1.48.0      
##  [45] xfun_0.16                stringr_1.4.0           
##  [47] lifecycle_0.2.0          XML_3.99-0.5            
##  [49] edgeR_3.30.3             zlibbioc_1.34.0         
##  [51] MASS_7.3-52              scales_1.1.1            
##  [53] BSgenome_1.56.0          VariantAnnotation_1.34.0
##  [55] hms_0.5.3                RBGL_1.64.0             
##  [57] RColorBrewer_1.1-2       yaml_2.2.1              
##  [59] curl_4.3                 memoise_1.1.0           
##  [61] ggplot2_3.3.2            biomaRt_2.44.1          
##  [63] latticeExtra_0.6-29      stringi_1.4.6           
##  [65] RSQLite_2.2.0            genefilter_1.70.0       
##  [67] desc_1.2.0               checkmate_2.0.0         
##  [69] GenomicFeatures_1.40.1   DOT_0.1                 
##  [71] rlang_0.4.7              pkgconfig_2.0.3         
##  [73] bitops_1.0-6             evaluate_0.14           
##  [75] lattice_0.20-41          purrr_0.3.4             
##  [77] bit_4.0.4                tidyselect_1.1.0        
##  [79] GSEABase_1.50.1          AnnotationForge_1.30.1  
##  [81] magrittr_1.5             bookdown_0.20           
##  [83] R6_2.4.1                 generics_0.0.2          
##  [85] base64url_1.4            DBI_1.1.0               
##  [87] withr_2.2.0              pillar_1.4.6            
##  [89] survival_3.2-3           RCurl_1.98-1.2          
##  [91] tibble_3.0.3             crayon_1.3.4            
##  [93] BiocFileCache_1.12.1     rmarkdown_2.3           
##  [95] jpeg_0.1-8.1             progress_1.2.2          
##  [97] locfit_1.5-9.4           grid_4.0.0              
##  [99] data.table_1.13.0        blob_1.2.1              
## [101] Rgraphviz_2.32.0         digest_0.6.25           
## [103] xtable_1.8-4             brew_1.0-6              
## [105] openssl_1.4.2            munsell_0.5.0           
## [107] askpass_1.1

Funding

This project was supported by funds from the National Institutes of Health (NIH) and the National Science Foundation (NSF).

References

H Backman, Tyler W, and Thomas Girke. 2016. “systemPipeR: NGS workflow and report generation environment.” BMC Bioinformatics 17 (1): 388. https://doi.org/10.1186/s12859-016-1241-0.

Langmead, Ben, and Steven L Salzberg. 2012. “Fast Gapped-Read Alignment with Bowtie 2.” Nat. Methods 9 (4): 357–59. https://doi.org/10.1038/nmeth.1923.

Love, Michael, Wolfgang Huber, and Simon Anders. 2014. “Moderated Estimation of Fold Change and Dispersion for RNA-seq Data with DESeq2.” Genome Biol. 15 (12): 550. https://doi.org/10.1186/s13059-014-0550-8.

Robinson, M D, D J McCarthy, and G K Smyth. 2010. “EdgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data.” Bioinformatics 26 (1): 139–40. https://doi.org/10.1093/bioinformatics/btp616.

Yu, Guangchuang, Li-Gen Wang, and Qing-Yu He. 2015. “ChIPseeker: An R/Bioconductor Package for ChIP Peak Annotation, Comparison and Visualization.” Bioinformatics 31 (14): 2382–3. https://doi.org/10.1093/bioinformatics/btv145.

Zhang, Y, T Liu, C A Meyer, J Eeckhoute, D S Johnson, B E Bernstein, C Nussbaum, et al. 2008. “Model-Based Analysis of ChIP-Seq (MACS).” Genome Biol. 9 (9). https://doi.org/10.1186/gb-2008-9-9-r137.

Zhu, Lihua J, Claude Gazin, Nathan D Lawson, Hervé Pagès, Simon M Lin, David S Lapointe, and Michael R Green. 2010. “ChIPpeakAnno: A Bioconductor Package to Annotate ChIP-seq and ChIP-chip Data.” BMC Bioinformatics 11: 237. https://doi.org/10.1186/1471-2105-11-237.