systemPipeR 1.20.0
Users want to provide here background information about the design of their ChIP-Seq project.
This report describes the analysis of several ChIP-Seq experiments studying the DNA binding patterns of the transcriptions factors … from organism ….
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.
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.
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.
After opening the Rmd
file of this workflow in Vim and attaching a connected
R session via the F2
(or other) key, use the following command sequence to run your R
session on a computer node.
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.6.0
R
Now check whether your R session is running on a computer node of the cluster and assess your environment.
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
The systemPipeR
package needs to be loaded to perform the analysis steps shown in
this report (H Backman and Girke 2016).
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 objects are masked from 'package:base':
##
## aperm, apply, rowsum
##
##
## Attaching package: 'systemPipeR'
## The following object is masked from 'package:BiocStyle':
##
## output
If applicable users can load custom functions not provided by systemPipeR
. Skip
this step if this is not the case.
source("systemPipeChIPseq_Fct.R")
targets
fileThe 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.gz M1A M1 Mock.1h.A
## 2 ./data/SRR446028_1.fastq.gz M1B M1 Mock.1h.B
## 3 ./data/SRR446029_1.fastq.gz A1A A1 Avr.1h.A
## 4 ./data/SRR446030_1.fastq.gz A1B A1 Avr.1h.B
## SampleReference
## 1
## 2
## 3 M1A
## 4 M1B
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-se",
package = "systemPipeR")
trim <- loadWF(targets = targetspath, wf_file = "trim-se.cwl",
input_file = "trim-se.yml", dir_path = dir_path)
trim <- renderWF(trim, inputvars = c(FileName = "_FASTQ_PATH1_",
SampleName = "_SampleName_"))
trim
output(trim)[1:2]
Next, we execute the code for trimming all the raw data.
# 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 = trim, Fct = "filterFct(fq, cutoff=20, Nexceptions=0)",
batchsize = 1e+05)
writeTargetsout(x = trim, file = "targets_chip_trim.txt", step = 1,
new_col = "FileName", new_col_output_index = 1, overwrite = TRUE)
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", "targets_chip.txt", package = "systemPipeR")
dir_path <- system.file("extdata/cwl/preprocessReads/trim-se",
package = "systemPipeR")
trim <- loadWorkflow(targets = targets, wf_file = "trim-se.cwl",
input_file = "trim-se.yml", dir_path = dir_path)
trim <- renderWF(trim, inputvars = c(FileName = "_FASTQ_PATH1_",
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()