1.1 Background and objectives

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

1.2 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.

2.1 Generate workflow environment

Load workflow environment with sample data into your current working directory. The sample data are described here.

library(systemPipeRdata)
genWorkenvir(workflow = "chipseq")
setwd("chipseq")

Alternatively, this can be done from the command-line as follows:

Rscript -e "systemPipeRdata::genWorkenvir(workflow='chipseq')"

In the workflow environments generated by genWorkenvir 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.

2.2 Run workflow

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

q("no")  # closes R session on head node

srun --x11 --partition=short --mem=2gb --cpus-per-task 4 --ntasks 1 --time 2:00:00 --pty bash -l
module load R/3.4.2
R

system("hostname")  # should return name of a compute node starting with i or c 
getwd()  # checks current working directory of R session
dir()  # returns content of current working directory

library(systemPipeR)

## Loading required package: Rsamtools

## Loading required package: GenomeInfoDb

## 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':
## 
##     Filter, Find, Map, Position, Reduce,
##     anyDuplicated, append, as.data.frame, basename,
##     cbind, colnames, dirname, do.call, duplicated,
##     eval, evalq, get, grep, grepl, intersect,
##     is.unsorted, lapply, mapply, match, mget, order,
##     paste, pmax, pmax.int, pmin, pmin.int, rank,
##     rbind, rownames, sapply, setdiff, sort, table,
##     tapply, union, unique, unsplit, which, which.max,
##     which.min

## Loading required package: S4Vectors

## Loading required package: stats4

## 
## Attaching package: 'S4Vectors'

## The following object is masked from 'package:base':
## 
##     expand.grid

## Loading required package: IRanges

## Loading required package: GenomicRanges

## Loading required package: Biostrings

## Loading required package: XVector

## 
## Attaching package: 'Biostrings'

## The following object is masked from 'package:base':
## 
##     strsplit

## Loading required package: ShortRead

## Loading required package: BiocParallel

## Loading required package: 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:Biostrings':
## 
##     type

## The following objects are masked from 'package:base':
## 
##     aperm, apply, rowsum

## Registered S3 methods overwritten by 'ggplot2':
##   method         from 
##   [.quosures     rlang
##   c.quosures     rlang
##   print.quosures rlang

## 

## 
## Attaching package: 'systemPipeR'

## The following object is masked from 'package:BiocStyle':
## 
##     output

source("systemPipeChIPseq_Fct.R")

3.1 Experiment definition provided by targets file

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

targetspath <- system.file("extdata", "targets_chip.txt", package = "systemPipeR")
targets <- read.delim(targetspath, comment.char = "#")
targets[1:4, -c(5, 6)]

##                   FileName SampleName Factor SampleLong
## 1 ./data/SRR446027_1.fastq        M1A     M1  Mock.1h.A
## 2 ./data/SRR446028_1.fastq        M1B     M1  Mock.1h.B
## 3 ./data/SRR446029_1.fastq        A1A     A1   Avr.1h.A
## 4 ./data/SRR446030_1.fastq        A1B     A1   Avr.1h.B
##   SampleReference
## 1                
## 2                
## 3             M1A
## 4             M1B

3.2 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 SYSargs instance (args object below). More detailed information on read preprocessing is provided in systemPipeR's main vignette.

args <- systemArgs(sysma = "param/trim.param", mytargets = "targets_chip.txt")
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 = args, Fct = "filterFct(fq, cutoff=20, Nexceptions=0)", 
    batchsize = 1e+05)
writeTargetsout(x = args, file = "targets_chip_trim.txt", overwrite = TRUE)

3.3 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.

args <- systemArgs(sysma = "param/tophat.param", mytargets = "targets_chip.txt")
library(BiocParallel)
library(batchtools)
f <- function(x) {
    library(systemPipeR)
    args <- systemArgs(sysma = "param/tophat.param", mytargets = "targets_chip.txt")
    seeFastq(fastq = infile1(args)[x], batchsize = 1e+05, klength = 8)
}
moduleload(modules(args))  # Skip if a module system is not used
resources <- list(walltime = 120, ntasks = 1, ncpus = cores(args), 
    memory = 1024)
param <- BatchtoolsParam(workers = 4, cluster = "slurm", template = "batchtools.slurm.tmpl", 
    resources = resources)
fqlist <- bplapply(seq(along = args), f, BPPARAM = param)

pdf("./results/fastqReport.pdf", height = 18, width = 4 * length(fqlist))
seeFastqPlot(unlist(fqlist, recursive = FALSE))
dev.off()