---
title: "ChIP-Seq Workflow Template"
author: "Author: Daniela Cassol (danielac@ucr.edu) and Thomas Girke (thomas.girke@ucr.edu)"
date: "Last update: `r format(Sys.time(), '%d %B, %Y')`"
output:
BiocStyle::html_document:
toc_float: true
code_folding: show
package: systemPipeR
vignette: |
%\VignetteEncoding{UTF-8}
%\VignetteIndexEntry{WF: ChIP-Seq Workflow Template}
%\VignetteEngine{knitr::rmarkdown}
fontsize: 14pt
bibliography: bibtex.bib
---
```{css, echo=FALSE}
pre code {
white-space: pre !important;
overflow-x: scroll !important;
word-break: keep-all !important;
word-wrap: initial !important;
}
```
```{r style, echo = FALSE, results = 'asis'}
BiocStyle::markdown()
options(width=60, max.print=1000)
knitr::opts_chunk$set(
eval=as.logical(Sys.getenv("KNITR_EVAL", "TRUE")),
cache=as.logical(Sys.getenv("KNITR_CACHE", "TRUE")),
tidy.opts=list(width.cutoff=60), tidy=TRUE)
```
```{r setup, echo=FALSE, message=FALSE, wwarning=FALSE, eval=FALSE}
suppressPackageStartupMessages({
library(systemPipeR)
library(BiocParallel)
library(Biostrings)
library(Rsamtools)
library(GenomicRanges)
library(ggplot2)
library(GenomicAlignments)
library(ShortRead)
library(ape)
library(batchtools)
})
```
# 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_Backman2016-bt]. 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](http://www.bioconductor.org/packages/devel/bioc/vignettes/systemPipeR/inst/doc/systemPipeR.html).
```{r load_systempiper, eval=TRUE, message=FALSE, warning=FALSE}
library(systemPipeR)
```
## Generate workflow environment
[*systemPipeRdata*](http://bioconductor.org/packages/release/data/experiment/html/systemPipeRdata.html) package is a helper package to generate a fully populated [*systemPipeR*](http://bioconductor.org/packages/release/bioc/html/systemPipeR.html)
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](http://www.bioconductor.org/packages/devel/data/experiment/vignettes/systemPipeRdata/inst/doc/systemPipeRdata.html#1_Introduction).
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](http://www.bioconductor.org/packages/release/bioc/vignettes/systemPipeR/inst/doc/systemPipeR.html). More information about the `targets` files from *systemPipeR* can be found [here](http://www.bioconductor.org/packages/release/bioc/vignettes/systemPipeR/inst/doc/systemPipeR.html#25_structure_of_targets_file).
## Run workflow
Now open the R markdown script `systemPipeRIBOseq.Rmd`in 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.
```{r load_targets_file, eval=TRUE}
targetspath <- system.file("extdata", "targetsPE_chip.txt", package="systemPipeR")
targets <- read.delim(targetspath, comment.char = "#")
targets[1:4,-c(5,6)]
```
## 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.
```{r construct_SYSargs2_trim-se, eval=FALSE}
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.
```{r proprocess_reads, eval=FALSE, message=FALSE, warning=FALSE, cache=TRUE}
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=100000)
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.
```{r fastq_report, eval=FALSE}
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=100000, 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 [@Langmead2012-bs]. 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:
```{r bowtie2_index, eval=FALSE}
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.
```{r bowtie2_align, eval=FALSE}
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]
```
```{r bowtie2_align_cluster, eval=FALSE}
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.
```{r bowtie2_align_single, eval=FALSE}
args <- runCommandline(args, force=F)
```
Check whether all BAM files have been created and write out the new targets file.
```{r check_files_exist, eval=FALSE}
writeTargetsout(x=args, file="targets_bam.txt", step = 1,
new_col = "FileName", new_col_output_index = 1, overwrite = TRUE, remove=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.
```{r align_stats, eval=FALSE}
read_statsDF <- alignStats(args=args)
write.table(read_statsDF, "results/alignStats.xls", row.names=FALSE, quote=FALSE, sep="\t")
read.delim("results/alignStats.xls")
```
## Create symbolic links for viewing BAM files in IGV
The `symLink2bam` function creates symbolic links to view the BAM alignment files in a
genome browser such as IGV without moving these large files to a local
system. The corresponding URLs are written to a file with a path
specified under `urlfile`, here `IGVurl.txt`. Please replace the directory and the user name.
```{r symbol_links, eval=FALSE}
symLink2bam(sysargs=args, htmldir=c("~/.html/", "somedir/"),
urlbase="http://cluster.hpcc.ucr.edu/~tgirke/",
urlfile="./results/IGVurl.txt")
```
# 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
```{r rle_object, eval=FALSE}
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
```{r resize_align, eval=FALSE}
trim(resize(as(aligns, "GRanges"), width = 200))
```
## Naive peak calling
```{r rle_slice, eval=FALSE}
islands <- slice(cov, lower = 15)
islands[[1]]
```
## Plot coverage for defined region
```{r plot_coverage, eval=FALSE}
library(ggbio)
myloc <- c("Chr1", 1, 100000)
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.
```{r merge_bams, eval=FALSE}
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, remove=TRUE)
```
## Peak calling without input/reference sample
MACS2 can perform peak calling on ChIP-Seq data with and without input
samples [@Zhang2008-pc]. 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).
```{r call_peaks_macs_noref, eval=FALSE}
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.
```{r call_peaks_macs_withref, eval=FALSE}
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.
```{r consensus_peaks, eval=FALSE}
# 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 [@Zhu2010-zo; @Yu2015-xu].
```{r chip_peak_anno, eval=FALSE}
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 )
```
```{r chip_peak_anno_full_annotation, include=FALSE, eval=FALSE}
## Perform previous step with full genome annotation from Biomart
# txdb <- makeTxDbFromBiomart(biomart = "plants_mart", dataset = "athaliana_eg_gene", host="plants.ensembl.org")
# tx <- transcripts(txdb, columns=c("tx_name", "gene_id", "tx_type"))
# ge <- genes(txdb, columns=c("tx_name", "gene_id", "tx_type")) # works as well
# seqlevels(ge) <- c("Chr1", "Chr2", "Chr3", "Chr4", "Chr5", "ChrC", "ChrM")
# table(mcols(tx)$tx_type)
# tx <- tx[!duplicated(unstrsplit(values(tx)$gene_id, sep=","))] # Keeps only first transcript model for each gene]
# annotatedPeak <- annotatePeakInBatch(peaksGR, AnnotationData = tx)
```
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.
```{r chip_peak_seeker, eval=FALSE}
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.
```{r chip_peak_seeker_plots, eval=FALSE}
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.
```{r count_peak_ranges, eval=FALSE}
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(FileName1 = "_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` [@Robinson2010-uk; @Love2014-sh].
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.
```{r diff_bind_analysis, eval=FALSE}
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.
```{r go_enrich, eval=FALSE}
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`.
```{r parse_peak_sequences, eval=FALSE}
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.
```{r bcrank_enrich, eval=FALSE}
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 sessionInfo}
sessionInfo()
```
# Funding
This project was supported by funds from the National Institutes of
Health (NIH) and the National Science Foundation (NSF).
# References