diff --git a/.gitignore b/.gitignore
index 113294a5f18d979085b4ad570d8a22a8e4eabd58..bf385c65b78082fb3581d2a2f613ae4da90d7d4b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -97,3 +97,16 @@ ENV/
# mypy
.mypy_cache/
+
+# Mac OS
+.DS_Store
+
+# nextflow analysis folders/files
+pipeline_trace*.txt*
+.nextflow*.log*
+report*.html*
+timeline*.html*
+/workflow/output/*
+/work/*
+/test_data/*
+/.nextflow/*
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 384eba4fe7f924a44609c76cfcbf0796b4390871..67829760655698ee58de138d71f224c637425807 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -1,14 +1,15 @@
before_script:
- module add python/3.6.1-2-anaconda
- - pip install --user pytest-pythonpath pytest-cov
+ - pip install --user pytest-pythonpath==0.7.1 pytest-cov==2.5.1
- module load nextflow/0.31.0
- - ln -s /work/BICF/s163035/chipseq/*fastq.gz test_data/
+ - ln -s /project/shared/bicf_workflow_ref/workflow_testdata/chipseq/*fastq.gz test_data/
stages:
- unit
- astrocyte
- single
- multiple
+ - skip
user_configuration:
stage: unit
@@ -28,32 +29,46 @@ astrocyte:
single_end_mouse:
stage: single
+ only:
+ - master
script:
- - nextflow run workflow/main.nf -resume
+ - nextflow run workflow/main.nf --astrocyte true -resume
- pytest -m singleend
- artifacts:
- expire_in: 2 days
paired_end_human:
stage: single
+ only:
+ - branches
+ except:
+ - master
script:
- - nextflow run workflow/main.nf --designFile "$CI_PROJECT_DIR/test_data/design_ENCSR729LGA_PE.txt" --genome 'GRCh38' --pairedEnd true -resume
+ - nextflow run workflow/main.nf --designFile "$CI_PROJECT_DIR/test_data/design_ENCSR729LGA_PE.txt" --genome 'GRCh38' --pairedEnd true --astrocyte false -resume
- pytest -m pairedend
- artifacts:
- expire_in: 2 days
single_end_diff:
stage: multiple
+ only:
+ - branches
+ except:
+ - master
script:
- - nextflow run workflow/main.nf --designFile "$CI_PROJECT_DIR/test_data/design_diff_SE.txt" --genome 'GRCm38' -resume
+ - nextflow run workflow/main.nf --designFile "$CI_PROJECT_DIR/test_data/design_diff_SE.txt" --genome 'GRCm38' --astrocyte false -resume
+ - pytest -m singleend
- pytest -m singlediff
- artifacts:
- expire_in: 2 days
paired_end_diff:
+ only:
+ - master
stage: multiple
script:
- - nextflow run workflow/main.nf --designFile "$CI_PROJECT_DIR/test_data/design_diff_PE.txt" --genome 'GRCh38' --pairedEnd true -resume
+ - nextflow run workflow/main.nf --designFile "$CI_PROJECT_DIR/test_data/design_diff_PE.txt" --genome 'GRCh38' --pairedEnd true --astrocyte false -resume
+ - pytest -m pairedend
- pytest -m paireddiff
- artifacts:
- expire_in: 2 days
+
+single_end_skip:
+ stage: skip
+ only:
+ - master
+ script:
+ - nextflow run workflow/main.nf --designFile "$CI_PROJECT_DIR/test_data/design_diff_SE.txt" --genome 'GRCm38' --skipDiff true --skipMotif true --skipPlotProfile true --astrocyte false -resume
+ - pytest -m singleskip_true
diff --git a/.gitlab/merge_request_templates/merge_request.md b/.gitlab/merge_request_templates/merge_request.md
new file mode 100644
index 0000000000000000000000000000000000000000..8ca7718c080cc2c047c083ed764a2051828787f6
--- /dev/null
+++ b/.gitlab/merge_request_templates/merge_request.md
@@ -0,0 +1,11 @@
+
+Please fill in the appropriate checklist below (delete whatever is not relevant).
+These are the most common things requested on pull requests (PRs).
+
+## PR checklist
+ - [ ] This comment contains a description of changes (with reason)
+ - [ ] If you've fixed a bug or added code that should be tested, add tests!
+ - [ ] Documentation in `docs` is updated
+ - [ ] `CHANGELOG.md` is updated
+ - [ ] `README.md` is updated
+ - [ ] `LICENSE.md` is updated with new contributors
diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
index 0000000000000000000000000000000000000000..b2f7b74489fa971ebe155da47a24f10fadf41cbc
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,40 @@
+# Changelog
+
+All notable changes to this project will be documented in this file.
+
+## [Unreleased]
+- Fix references.md link in citation of README.md
+- Add Nextflow to references.md
+- Fix pool_and_psuedoreplicate.py to run single experiment
+- Add test data for test_pool_and_pseudoreplicate
+- Add PlotProfile Option
+- Add Python version to MultiQC
+- Add and Update tests
+- Use GTF files instead of TxDb and org libraries in Annotate Peaks
+- Make gtf and geneName files as param inputs
+- Fix xcor to increase file size for --random-source
+- Fix skip diff test for paired-end data
+
+## [publish_1.0.6 ] - 2019-05-31
+### Added
+- MIT LICENSE
+- Check for correct genome input
+- Check for correct output for motif serach
+
+### Fixed
+- Path to reference fasta file
+
+## [publish_1.0.5 ] - 2019-05-16
+### Fixed
+- Retitled documentation for skipDiff and skipMotif to be more clear
+- Add missing python module for motif search
+- Updated links for phantompeaktools and references
+- Fix callPeaks for single control experiments
+
+## [publish_1.0.3 ] - 2019-04-23
+### Fixed
+- Variable name for design file for Astrocyte
+- File path for design file to work for Astrocyte
+
+## [publish_1.0.0 ] - 2019-04-23
+Initial release of pipeline
diff --git a/LICENSE.md b/LICENSE.md
new file mode 100644
index 0000000000000000000000000000000000000000..5b4b67c8e0190980e6cdcaf7606476815b647c9b
--- /dev/null
+++ b/LICENSE.md
@@ -0,0 +1,15 @@
+MIT License
+
+Copyright (c) 2019, University of Texas Southwestern Medical Center.
+
+All rights reserved.
+
+Contributors: Spencer D. Barnes, Holly Ruess, Jeremy A. Mathews, Beibei Chen, Venkat S. Malladi
+
+Department: Bioinformatic Core Facility, Department of Bioinformatics
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
diff --git a/README.md b/README.md
index 0cbc9d47de6d7ba639749a8353f6fa121c825022..ec8e8c7908178f3db77dfc753fec25b483d65a50 100644
--- a/README.md
+++ b/README.md
@@ -1,3 +1,7 @@
+# **CHIPseq Manual**
+## Version 1.0.6
+## May 31, 2019
+
# BICF ChIP-seq Pipeline
[](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/commits/master)
@@ -5,13 +9,150 @@
[](https://www.nextflow.io/)
[](https://astrocyte-test.biohpc.swmed.edu/static/docs/index.html)
+[](https://doi.org/10.5281/zenodo.2648845)
## Introduction
-BICF ChIPseq is a bioinformatics best-practice analysis pipeline used for ChIP-seq (chromatin immunoprecipitation sequencing) data analysis at [BICF](http://www.utsouthwestern.edu/labs/bioinformatics/) at [UT Southwestern Dept. of Bioinformatics](http://www.utsouthwestern.edu/departments/bioinformatics/).
+BICF ChIPseq is a bioinformatics best-practice analysis pipeline used for ChIP-seq (chromatin immunoprecipitation sequencing) data analysis at [BICF](http://www.utsouthwestern.edu/labs/bioinformatics/) at [UT Southwestern Department of Bioinformatics](http://www.utsouthwestern.edu/departments/bioinformatics/).
The pipeline uses [Nextflow](https://www.nextflow.io), a bioinformatics workflow tool. It pre-processes raw data from FastQ inputs, aligns the reads and performs extensive quality-control on the results.
-This pipeline is primarily used with a SLURM cluster on the [BioHPC Cluster](https://biohpc.swmed.edu/). However, the pipeline should be able to run on any system that Nextflow supports.
+This pipeline is primarily used with a SLURM cluster on the [BioHPC Cluster](https://portal.biohpc.swmed.edu/content/). However, the pipeline should be able to run on any system that supports Nextflow.
+
+Additionally, the pipeline is designed to work with [Astrocyte Workflow System](https://astrocyte.biohpc.swmed.edu/static/docs/index.html) using a simple web interface.
+
+Current version of the software and issue reports are at
+https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis
+
+To download the current version of the software
+```bash
+$ git clone git@git.biohpc.swmed.edu:BICF/Astrocyte/chipseq_analysis.git
+```
+
+## Input files
+##### 1) Fastq Files
+ + You will need the full path to the files for the Bash Scipt
+
+##### 2) Design File
+ + The Design file is a tab-delimited file with 8 columns for Single-End and 9 columns for Paired-End. Letter, numbers, and underlines can be used in the names. However, the names can only begin with a letter. Columns must be as follows:
+ 1. sample_id a short, unique, and concise name used to label output files; will be used as a control_id if it is the control sample
+ 2. experiment_id biosample_treatment_factor; same name given for all replicates of treatment. Will be used for the consensus header.
+ 3. biosample symbol for tissue type or cell line
+ 4. factor symbol for antibody target
+ 5. treatment symbol of treatment applied
+ 6. replicate a number, usually from 1-3 (i.e. 1)
+ 7. control_id sample_id name that is the control for this sample
+ 8. fastq_read1 name of fastq file 1 for SE or PC data
+ 9. fastq_read2 name of fastq file 2 for PE data
+
+
+ + See [HERE](test_data/design_ENCSR729LGA_PE.txt) for an example design file, paired-end
+ + See [HERE](test_data/design_ENCSR238SGC_SE.txt) for an example design file, single-end
+
+##### 3) Bash Script
+ + You will need to create a bash script to run the CHIPseq pipeline on [BioHPC](https://portal.biohpc.swmed.edu/content/)
+ + This pipeline has been optimized for the correct partition
+ + See [HERE](docs/CHIPseq.sh) for an example bash script
+ + The parameters that must be specified are:
+ - --reads '/path/to/files/name.fastq.gz'
+ - --designFile '/path/to/file/design.txt',
+ - --genome 'GRCm38', 'GRCh38', or 'GRCh37' (if you need to use another genome contact the [BICF](mailto:BICF@UTSouthwestern.edu))
+ - --pairedEnd 'true' or 'false' (where 'true' is PE and 'false' is SE; default 'false')
+ - --outDir (optional) path and folder name of the output data, example: /home2/s000000/Desktop/Chipseq_output (if not specficied will be under workflow/output/)
+
+## Pipeline
+ + There are 11 steps to the pipeline
+ 1. Check input files
+ 2. Trim adaptors TrimGalore!
+ 3. Aligned trimmed reads with bwa, and sorts/converts to bam with samtools
+ 4. Mark duplicates with Sambamba, and filter reads with samtools
+ 5. Quality metrics with deep tools
+ 6. Calculate cross-correlation using PhantomPeakQualTools
+ 7. Call peaks with MACS
+ 8. Calculate consensus peaks
+ 9. Annotate all peaks using ChipSeeker
+ 10. Calculate Differential Binding Activity with DiffBind (If more than 1 rep in more than 1 experiment)
+ 11. Use MEME-ChIP to find motifs in original peaks
+
+See [FLOWCHART](docs/flowchart.pdf)
+
+## Output Files
+Folder | File | Description
+--- | --- | ---
+design | N/A | Inputs used for analysis; can ignore
+trimReads | *_trimming_report.txt | report detailing how many reads were trimmed
+trimReads | *_trimmed.fq.gz | trimmed fastq files used for analysis
+alignReads | *.srt.bam.flagstat.qc | QC metrics from the mapping process
+alignReads | *.srt.bam | sorted bam file
+filterReads | *.dup.qc | QC metrics of find duplicate reads (sambamba)
+filterReads | *.filt.nodup.bam | filtered bam file with duplicate reads removed
+filterReads | *.filt.nodup.bam.bai | indexed filtered bam file
+filterReads | *.filt.nodup.flagstat.qc | QC metrics of filtered bam file (mapping stats, samtools)
+filterReads | *.filt.nodup.pbc.qc | QC metrics of library complexity
+convertReads | *.filt.nodup.bedse.gz | bed alignment in BEDPE format
+convertReads | *.filt.nodup.tagAlign.gz | bed alignent in BEDPE format, same as bedse unless samples are paired-end
+multiqcReport | multiqc_report.html | Quality control report of NRF, PBC1, PBC2, NSC, and RSC. Also contains software versions and references to cite.
+experimentQC | coverage.pdf | plot to assess the sequencing depth of a given sample
+experimentQC | *_fingerprint.pdf | plot to determine if the antibody-treatment enriched sufficiently
+experimentQC | heatmeap_SpearmanCorr.pdf | plot of Spearman correlation between samples
+experimentQC | heatmeap_PearsonCorr.pdf | plot of Pearson correlation between samples
+experimentQC | sample_mbs.npz | array of multiple BAM summaries
+crossReads | *.cc.plot.pdf | Plot of cross-correlation to assess signal-to-noise ratios
+crossReads | *.cc.qc | cross-correlation metrics. File [HEADER](docs/xcor_header.txt)
+callPeaksMACS | pooled/*pooled.fc_signal.bw | bigwig data file; raw fold enrichment of sample/control
+callPeaksMACS | pooled/*pooled_peaks.xls | Excel file of peaks
+callPeaksMACS | pooled/*.pvalue_signal.bw | bigwig data file; sample/control signal adjusted for pvalue significance
+callPeaksMACS | pooled/*_pooled.narrowPeak | peaks file; see [HERE](https://genome.ucsc.edu/FAQ/FAQformat.html#format12) for ENCODE narrowPeak header format
+consensusPeaks | *.rejected.narrowPeak | peaks not supported by multiple testing (replicates and pseudo-replicates)
+consensusPeaks | *.replicated.narrowPeak | peaks supported by multiple testing (replicates and pseudo-replicates)
+peakAnnotation | *.chipseeker_annotation.tsv | annotated narrowPeaks file
+peakAnnotation | *.chipseeker_pie.pdf | pie graph of where narrow annotated peaks occur
+peakAnnotation | *.chipseeker_upsetplot.pdf | upsetplot showing the count of overlaps of the genes with different annotated location
+motifSearch | *_memechip/index.html | interactive HTML link of MEME output
+motifSearch | sorted-*.replicated.narrowPeak | Top 600 peaks sorted by p-value; input for motifSearch
+motifSearch | *_memechip/combined.meme | MEME identified motifs
+diffPeaks | heatmap.pdf | Use only for replicated samples; heatmap of relationship of peak location and peak intensity
+diffPeaks | normcount_peaksets.txt | Use only for replicated samples; peak set values of each sample
+diffPeaks | pca.pdf | Use only for replicated samples; PCA of peak location and peak intensity
+diffPeaks | *_diffbind.bed | Use only for replicated samples; bed file of peak locations between replicates
+diffPeaks | *_diffbind.csv | Use only for replicated samples; CSV file of peaks between replicates
+plotProfile | plotProfile.png | Plot profile of the TSS region
+plotProfile | computeMatrix.gz | Compute Matrix from deeptools to create custom plots other than plotProfile
+
+## Common Quality Control Metrics
+ + These are the list of files that should be reviewed before continuing on with the CHIPseq experiment. If your experiment fails any of these metrics, you should pause and re-evaluate whether the data should remain in the study.
+ 1. multiqcReport/multiqc_report.html: follow the ChiP-seq standards [HERE](https://www.encodeproject.org/chip-seq/);
+ 2. experimentQC/*_fingerprint.pdf: make sure the plots information is correct for your antibody/input. See [HERE](https://deeptools.readthedocs.io/en/develop/content/tools/plotFingerprint.html) for more details.
+ 3. crossReads/*cc.plot.pdf: make sure your sample data has the correct signal intensity and location. See [HERE](hhttps://ccg.epfl.ch//var/sib_april15/cases/landt12/strand_correlation.html) for more details.
+ 4. crossReads/*.cc.qc: Column 9 (NSC) should be > 1.1 for experiment and < 1.1 for input. Column 10 (RSC) should be > 0.8 for experiment and < 0.8 for input. See [HERE](https://genome.ucsc.edu/encode/qualityMetrics.html) for more details.
+ 5. experimentQC/coverage.pdf, experimentQC/heatmeap_SpearmanCorr.pdf, experimentQC/heatmeap_PearsonCorr.pdf: See [HERE](https://deeptools.readthedocs.io/en/develop/content/list_of_tools.html) for more details.
+
+
+## Common Errors
+If you find an error, please let the [BICF](mailto:BICF@UTSouthwestern.edu) know and we will add it here.
+
+## Citation
+Please cite individual programs and versions used [HERE](docs/references.md), and the pipeline doi:[10.5281/zenodo.2648844](https://doi.org/10.5281/zenodo.2648844). Please cite in publications: Pipeline was developed by BICF from funding provided by Cancer Prevention and Research Institute of Texas (RP150596).
+
+## Programs and Versions
+ + python/3.6.1-2-anaconda [website](https://www.anaconda.com/download/#linux) [citation](docs/references.md)
+ + trimgalore/0.4.1 [website](https://github.com/FelixKrueger/TrimGalore) [citation](docs/references.md)
+ + cutadapt/1.9.1 [website](https://cutadapt.readthedocs.io/en/stable/index.html) [citation](docs/references.md)
+ + bwa/intel/0.7.12 [website](http://bio-bwa.sourceforge.net/) [citation](docs/references.md)
+ + samtools/1.6 [website](http://samtools.sourceforge.net/) [citation](docs/references.md)
+ + sambamba/0.6.6 [website](http://lomereiter.github.io/sambamba/) [citation](docs/references.md)
+ + bedtools/2.26.0 [website](https://bedtools.readthedocs.io/en/latest/) [citation](docs/references.md)
+ + deeptools/2.5.0.1 [website](https://deeptools.readthedocs.io/en/develop/) [citation](docs/references.md)
+ + phantompeakqualtools/1.2 [website](https://github.com/kundajelab/phantompeakqualtools) [citation](docs/references.md)
+ + macs/2.1.0-20151222 [website](http://liulab.dfci.harvard.edu/MACS/) [citation](docs/references.md)
+ + UCSC_userApps/v317 [website](https://genome.ucsc.edu/util.html) [citation](docs/references.md)
+ + R/3.4.1 [website](https://www.r-project.org/) [citation](docs/references.md)
+ + SPP/1.14
+ + meme/4.11.1-gcc-openmpi [website](http://meme-suite.org/doc/install.html?man_type=web) [citation](docs/references.md)
+ + ChIPseeker [website](https://bioconductor.org/packages/release/bioc/html/ChIPseeker.html) [citation](docs/references.md)
+ + DiffBind [website](https://bioconductor.org/packages/release/bioc/html/DiffBind.html) [citation](docs/references.md)
+
+
-Additionally, the pipeline is designed to work with [Astrocyte Workflow System](https://astrocyte-test.biohpc.swmed.edu/static/docs/index.html) using a simple web interface.
+## Credits
+This example worklow is derived from original scripts kindly contributed by the Bioinformatic Core Facility ([BICF](https://www.utsouthwestern.edu/labs/bioinformatics/)), in the [Department of Bioinformatics](https://www.utsouthwestern.edu/departments/bioinformatics/).
diff --git a/astrocyte_pkg.yml b/astrocyte_pkg.yml
index e1a1846326e06a833f85e0fa6bd19db836fb43f3..704bc8fa415a0e624163971f225f034115f4c9f3 100644
--- a/astrocyte_pkg.yml
+++ b/astrocyte_pkg.yml
@@ -9,7 +9,7 @@
# A unique identifier for the workflow package, text/underscores only
name: 'chipseq_analysis_bicf'
# Who wrote this?
-author: 'Beibei Chen and Venkat Malladi'
+author: 'Holly Ruess, Spencer D. Barnes, Beibei Chen and Venkat Malladi'
# A contact email address for questions
email: 'bicf@utsouthwestern.edu'
# A more informative title for the workflow package
@@ -51,6 +51,7 @@ workflow_modules:
- 'UCSC_userApps/v317'
- 'R/3.3.2-gccmkl'
- 'meme/4.11.1-gcc-openmpi'
+ - 'pandoc/2.7'
# A list of parameters used by the workflow, defining how to present them,
@@ -93,29 +94,25 @@ workflow_parameters:
description: |
One or more input FASTQ files from a ChIP-seq expereiment and a design
file with the link bewetwen the same file name and sample id
- regex: ".*(fastq|fq)*"
+ regex: ".*(fastq|fq)*gz"
min: 2
- id: pairedEnd
type: select
required: true
choices:
- - [ 'true', 'True']
- - [ 'false', 'False']
+ - [ 'true', 'true']
+ - [ 'false', 'false']
description: |
- In single-end sequencing, the sequencer reads a fragment from only one
- end to the other, generating the sequence of base pairs. In paired-end
- reading it starts at one read, finishes this direction at the specified
- read length, and then starts another round of reading from the opposite
- end of the fragment.
+ If Paired-end: True, if Single-end: False.
- - id: design
+ - id: designFile
type: file
required: true
description: |
A design file listing sample id, fastq files, corresponding control id
and additional information about the sample.
- regex: ".*tsv"
+ regex: ".*txt"
- id: genome
type: select
@@ -127,13 +124,41 @@ workflow_parameters:
description: |
Reference species and genome used for alignment and subsequent analysis.
+ - id: skipDiff
+ type: select
+ required: true
+ choices:
+ - [ 'true', 'true']
+ - [ 'false', 'false']
+ description: |
+ Skip differential peak analysis
+
+ - id: skipMotif
+ type: select
+ required: true
+ choices:
+ - [ 'true', 'true']
+ - [ 'false', 'false']
+ description: |
+ Skip motif calling
+
+ - id: skipPlotProfile
+ type: select
+ required: true
+ choices:
+ - [ 'true', 'true']
+ - [ 'false', 'false']
+ description: |
+ Skip Plot Profile Analysis
+
- id: astrocyte
- type: string
+ type: select
+ choices:
+ - [ 'true', 'true' ]
required: true
default: 'true'
- regex: "true"
description: |
- Ensure configuraton for astrocyte
+ Ensure configuraton for astrocyte.
# -----------------------------------------------------------------------------
@@ -154,8 +179,4 @@ vizapp_cran_packages:
# List of any Bioconductor packages, not provided by the modules,
# that must be made available to the vizapp
-vizapp_bioc_packages:
- - none
-
-vizapp_github_packages:
- - js229/Vennerable
+vizapp_bioc_packages: []
diff --git a/docs/CHIPseq.sh b/docs/CHIPseq.sh
new file mode 100644
index 0000000000000000000000000000000000000000..68aa87c7470ce983c9ede9a4918e0c9e3369ae19
--- /dev/null
+++ b/docs/CHIPseq.sh
@@ -0,0 +1,15 @@
+#!/bin/bash
+
+#SBATCH --job-name=CHIPseq
+#SBATCH --partition=super
+#SBATCH --output=CHIPseq.%j.out
+#SBATCH --error=CHIPseq.%j.err
+
+module load nextflow/0.31.0
+module add python/3.6.1-2-anaconda
+
+nextflow run workflow/main.nf \
+--reads '/path/to/*fastq.gz' \
+--designFile '/path/to/design.txt' \
+--genome 'GRCm38' \
+--pairedEnd 'true'
diff --git a/docs/design_example.csv b/docs/design_example.csv
deleted file mode 100644
index 6fa48f3f3f61a8d70dada47f4cfa6c9c0bf1fbfd..0000000000000000000000000000000000000000
--- a/docs/design_example.csv
+++ /dev/null
@@ -1,7 +0,0 @@
-SampleID,Tissue,Factor,Condition,Replicate,Peaks,bamReads,bamControl,ControlID,PeakCaller
-A_1,A,H3K27AC,A,1,A_1.broadPeak,A_1.bam,A_1_input.bam,A_1_input,bed
-A_2,A,H3K27AC,A,2,A_2.broadPeak,A_2.bam,A_2_input.bam,A_2_input,bed
-B_1,B,H3K27AC,B,1,B_1.broadPeak,B_1.bam,B_1_input.bam,B_1_input,bed
-B_2,B,H3K27AC,B,2,B_2.broadPeak,B_2.bam,B_2_input.bam,B_2_input,bed
-C_1,C,H3K27AC,C,1,C_1.broadPeak,C_1.bam,C_1_input.bam,C_1_input,bed
-C_2,C,H3K27AC,C,2,C_2.broadPeak,C_2.bam,C_2_input.bam,C_2_input,bed
diff --git a/docs/design_example.txt b/docs/design_example.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3df687fb3e6ac323eb6967df717505f2323cfe49
--- /dev/null
+++ b/docs/design_example.txt
@@ -0,0 +1,5 @@
+sample_id experiment_id biosample factor treatment replicate control_id fastq_read1
+A1 A tissueA H3K27AC None 1 B1 A1.fastq.gz
+A2 A tissueA H3K27AC None 2 B2 A2.fastq.gz
+B1 B tissueB Input None 1 B1 B1.fastq.gz
+B2 B tissueB Input None 2 B2 B2.fastq.gz
diff --git a/docs/flowchart.pdf b/docs/flowchart.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..845c5fc6048d3b5eabea384a32cf1c75bb8c7024
Binary files /dev/null and b/docs/flowchart.pdf differ
diff --git a/docs/images/phantompeakqualtools.png b/docs/images/phantompeakqualtools.png
new file mode 100644
index 0000000000000000000000000000000000000000..356bec28bcd52272ffb36c3b23c3fab5ebc82ab5
Binary files /dev/null and b/docs/images/phantompeakqualtools.png differ
diff --git a/docs/index.md b/docs/index.md
index b6c66cec174ba4aeca2216c39f18dda3bfa89010..3b4a685cfc112291a221c04d7efa4fb85db9545a 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -1,65 +1,130 @@
-# Astrocyte ChIPseq analysis Workflow Package
+# BICF ChIP-seq Analysis Workflow
## Introduction
**ChIP-seq Analysis** is a bioinformatics best-practice analysis pipeline used for chromatin immunoprecipitation (ChIP-seq) data analysis.
The pipeline uses [Nextflow](https://www.nextflow.io), a bioinformatics workflow tool. It pre-processes raw data from FastQ inputs, aligns the reads and performs extensive quality-control on the results.
-### Pipeline Steps
+Report issues to the Bioinformatic Core Facility [BICF](mailto:BICF@UTSouthwestern.edu)
-1) Trim adaptors TrimGalore!
-2) Align with BWA
-3) Filter reads with Sambamba S
-4) Quality control with DeepTools
-5) Calculate Cross-correlation using SPP and PhantomPeakQualTools
-6) Signal profiling using MACS2
-7) Call consenus peaks
-8) Annotate all peaks using ChipSeeker
-9) Use MEME-ChIP to find motifs in original peaks
-10) Find differential expressed peaks using DiffBind (If more than 1 experiment)
+### Pipeline Steps
+ + There are 11 steps to the pipeline
+ 1. Check input files
+ 2. Trim adaptors TrimGalore!
+ 3. Aligned trimmed reads with bwa, and sorts/converts to bam with samtools
+ 4. Mark duplicates with Sambamba, and filter reads with samtools
+ 5. Quality metrics with deep tools
+ 6. Calculate cross-correlation using PhantomPeakQualTools
+ 7. Call peaks with MACS
+ 8. Calculate consensus peaks
+ 9. Annotate all peaks using ChipSeeker
+ 10. Calculate Differential Binding Activity with DiffBind (If more than 1 rep in more than 1 experiment)
+ 11. Use MEME-ChIP to find motifs in original peaks
## Workflow Parameters
-
- reads - Choose all ChIP-seq fastq files for analysis.
- pairedEnd - Choose True/False if data is paired-end
- design - Choose the file with the experiment design information. TSV format
+ 1. One or more input FASTQ files from a ChIP-seq expereiment and a design file with the link bewetwen the same file name and sample id (required) - Choose all ChIP-seq fastq files for analysis.
+ 2. In single-end sequencing, the sequencer reads a fragment from only one end to the other, generating the sequence of base pairs. In paired-end reading it starts at one read, finishes this direction at the specified read length, and then starts another round of reading from the opposite end of the fragment. (Paired-end: True, Single-end: False) (required)
+ 3. A design file listing sample id, fastq files, corresponding control id and additional information about the sample.
genome - Choose a genomic reference (genome).
-
+ 4. Reference species and genome used for alignment and subsequent analysis. (required)
+ 5. Run differential peak analysis (required). Must have at least 2 replicates per experiment and at least 2 experiments.
+ 6. Run motif calling (required). Top 600 peaks sorted by p-value.
+ 7. Ensure configuration for astrocyte. (required; always true)
## Design file
- The following columns are necessary, must be named as in template. An design file template can be downloaded [HERE](https://git.biohpc.swmed.edu/bchen4/chipseq_analysis/raw/master/docs/design_example.csv)
-
- SampleID
- The id of the sample. This will be the header in output files, please make sure it is concise
- Tissue
- Tissue of the sample
- Factor
- Factor of the experiment
- Condition
- This is the group that will be used for pairwise differential expression analysis
- Replicate
- Replicate id
- Peaks
- The file name of the peak file for this sample
- bamReads
- The file name of the IP BAM for this sample
- bamControl
- The file name of the control BAM for this sample
- ContorlID
- The id of the control sample
- PeakCaller
- The peak caller used
+ + The Design file is a tab-delimited file with 8 columns for Single-End and 9 columns for Paired-End. Letter, numbers, and underlines can be used in the names. However, the names can only begin with a letter. Columns must be as follows:
+ 1. sample_id a short, unique, and concise name used to label output files; will be used as a control_id if it is the control sample
+ 2. experiment_id biosample_treatment_factor; same name given for all replicates of treatment. Will be used for the consensus header.
+ 3. biosample symbol for tissue type or cell line
+ 4. factor symbol for antibody target
+ 5. treatment symbol of treatment applied
+ 6. replicate a number, usually from 1-3 (i.e. 1)
+ 7. control_id sample_id name that is the control for this sample
+ 8. fastq_read1 name of fastq file 1 for SE or PC data
+ 9. fastq_read2 name of fastq file 2 for PE data
+
+
+ + See [HERE](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/test_data/design_ENCSR729LGA_PE.txt) for an example design file, paired-end
+ + See [HERE](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/test_data/design_ENCSR238SGC_SE.txt) for an example design file, single-end
+
+## Output Files
+
+Folder | File | Description
+--- | --- | ---
+design | N/A | Inputs used for analysis; can ignore
+trimReads | *_trimming_report.txt | report detailing how many reads were trimmed
+trimReads | *_trimmed.fq.gz | trimmed fastq files used for analysis
+alignReads | *.srt.bam.flagstat.qc | QC metrics from the mapping process
+alignReads | *.srt.bam | sorted bam file
+filterReads | *.dup.qc | QC metrics of find duplicate reads (sambamba)
+filterReads | *.filt.nodup.bam | filtered bam file with duplicate reads removed
+filterReads | *.filt.nodup.bam.bai | indexed filtered bam file
+filterReads | *.filt.nodup.flagstat.qc | QC metrics of filtered bam file (mapping stats, samtools)
+filterReads | *.filt.nodup.pbc.qc | QC metrics of library complexity
+convertReads | *.filt.nodup.bedse.gz | bed alignment in BEDPE format
+convertReads | *.filt.nodup.tagAlign.gz | bed alignent in BEDPE format, same as bedse unless samples are paired-end
+multiqcReport | multiqc_report.html | Quality control report of NRF, PBC1, PBC2, NSC, and RSC. Also contains software versions and references to cite.
+experimentQC | coverage.pdf | plot to assess the sequencing depth of a given sample
+experimentQC | *_fingerprint.pdf | plot to determine if the antibody-treatment enriched sufficiently
+experimentQC | heatmeap_SpearmanCorr.pdf | plot of Spearman correlation between samples
+experimentQC | heatmeap_PearsonCorr.pdf | plot of Pearson correlation between samples
+experimentQC | sample_mbs.npz | array of multiple BAM summaries
+crossReads | *.cc.plot.pdf | Plot of cross-correlation to assess signal-to-noise ratios
+crossReads | *.cc.qc | cross-correlation metrics. File [HEADER](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/docs/xcor_header.txt)
+callPeaksMACS | pooled/*pooled.fc_signal.bw | bigwig data file; raw fold enrichment of sample/control
+callPeaksMACS | pooled/*pooled_peaks.xls | Excel file of peaks
+callPeaksMACS | pooled/*.pvalue_signal.bw | bigwig data file; sample/control signal adjusted for pvalue significance
+callPeaksMACS | pooled/*_pooled.narrowPeak | peaks file; see [HERE](https://genome.ucsc.edu/FAQ/FAQformat.html#format12) for ENCODE narrowPeak header format
+consensusPeaks | *.rejected.narrowPeak | peaks not supported by multiple testing (replicates and pseudo-replicates)
+consensusPeaks | *.replicated.narrowPeak | peaks supported by multiple testing (replicates and pseudo-replicates)
+peakAnnotation | *.chipseeker_annotation.tsv | annotated narrowPeaks file
+peakAnnotation | *.chipseeker_pie.pdf | pie graph of where narrow annotated peaks occur
+peakAnnotation | *.chipseeker_upsetplot.pdf | upsetplot showing the count of overlaps of the genes with different annotated location
+motifSearch | *_memechip/index.html | interactive HTML link of MEME output
+motifSearch | sorted-*.replicated.narrowPeak | Top 600 peaks sorted by p-value; input for motifSearch
+motifSearch | *_memechip/combined.meme | MEME identified motifs
+diffPeaks | heatmap.pdf | Use only for replicated samples; heatmap of relationship of peak location and peak intensity
+diffPeaks | normcount_peaksets.txt | Use only for replicated samples; peak set values of each sample
+diffPeaks | pca.pdf | Use only for replicated samples; PCA of peak location and peak intensity
+diffPeaks | *_diffbind.bed | Use only for replicated samples; bed file of peak locations between replicates
+diffPeaks | *_diffbind.csv | Use only for replicated samples; CSV file of peaks between replicates
+plotProfile | plotProfile.png | Plot profile of the TSS region
+plotProfile | computeMatrix.gz | Compute Matrix from deeptools to create custom plots other than plotProfile
+
+## Common Quality Control Metrics
+ + These are the list of files that should be reviewed before continuing on with the CHIPseq experiment. If your experiment fails any of these metrics, you should pause and re-evaluate whether the data should remain in the study.
+ 1. multiqcReport/multiqc_report.html: follow the ChiP-seq standards [HERE](https://www.encodeproject.org/chip-seq/);
+ 2. experimentQC/*_fingerprint.pdf: make sure the plots information is correct for your antibody/input. See [HERE](https://deeptools.readthedocs.io/en/develop/content/tools/plotFingerprint.html) for more details.
+ 3. crossReads/*cc.plot.pdf: make sure your sample data has the correct signal intensity and location. See [HERE](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/docs/phantompeaks.md) for more details.
+ 4. crossReads/*.cc.qc: Column 9 (NSC) should be > 1.1 for experiment and < 1.1 for input. Column 10 (RSC) should be > 0.8 for experiment and < 0.8 for input. See [HERE](https://genome.ucsc.edu/encode/qualityMetrics.html) for more details.
+ 5. experimentQC/coverage.pdf, experimentQC/heatmeap_SpearmanCorr.pdf, experimentQC/heatmeap_PearsonCorr.pdf: See [HERE](https://deeptools.readthedocs.io/en/develop/content/list_of_tools.html) for more details.
### Credits
-This example worklow is derived from original scripts kindly contributed by the Bioinformatic Core Facility (BICF), Department of Bioinformatics
+This example worklow is derived from original scripts kindly contributed by the Bioinformatic Core Facility ([BICF](https://www.utsouthwestern.edu/labs/bioinformatics/)), in the [Department of Bioinformatics](https://www.utsouthwestern.edu/departments/bioinformatics/).
+
+Please cite in publications: Pipeline was developed by BICF from funding provided by Cancer Prevention and Research Institute of Texas (RP150596).
### References
-* ChipSeeker: http://bioconductor.org/packages/release/bioc/html/ChIPseeker.html
-* DiffBind: http://bioconductor.org/packages/release/bioc/html/DiffBind.html
-* Deeptools: https://deeptools.github.io/
-* MEME-ChIP: http://meme-suite.org/doc/meme-chip.html
+ + python/3.6.1-2-anaconda [website](https://www.anaconda.com/download/#linux) [citation](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/docs/references.md)
+ + trimgalore/0.4.1 [website](https://github.com/FelixKrueger/TrimGalore) [citation](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/docs/references.md)
+ + cutadapt/1.9.1 [website](https://cutadapt.readthedocs.io/en/stable/index.html) [citation](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/docs/references.md)
+ + bwa/intel/0.7.12 [website](http://bio-bwa.sourceforge.net/) [citation](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/docs/references.md)
+ + samtools/1.6 [website](http://samtools.sourceforge.net/) [citation](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/docs/references.md)
+ + sambamba/0.6.6 [website](http://lomereiter.github.io/sambamba/) [citation](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/docs/references.md)
+ + bedtools/2.26.0 [website](https://bedtools.readthedocs.io/en/latest/) [citation](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/docs/references.md)
+ + deeptools/2.5.0.1 [website](https://deeptools.readthedocs.io/en/develop/) [citation](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/docs/references.md)
+ + phantompeakqualtools/1.2 [website](https://github.com/kundajelab/phantompeakqualtools) [citation](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/docs/references.md)
+ + macs/2.1.0-20151222 [website](http://liulab.dfci.harvard.edu/MACS/) [citation](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/docs/references.md)
+ + UCSC_userApps/v317 [website](https://genome.ucsc.edu/util.html) [citation](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/docs/references.md)
+ + R/3.3.2-gccmkl [website](https://www.r-project.org/) [citation](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/docs/references.md)
+ + meme/4.11.1-gcc-openmpi [website](http://meme-suite.org/doc/install.html?man_type=web) [citation](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/docs/references.md)
+ + ChIPseeker [website](https://bioconductor.org/packages/release/bioc/html/ChIPseeker.html) [citation](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/docs/references.md)
+ + DiffBind [website](https://bioconductor.org/packages/release/bioc/html/DiffBind.html) [citation](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/docs/references.md)
+ + MultiQC [website](https://multiqc.info/) [citation](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/docs/references.md)
+ + BICFChip-seqAnalysisWorkflow [website](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis) [citation](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/docs/references.md)
+ + Nextflow [website](https://www.nextflow.io/) [citation](https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/blob/master/docs/references.md)
diff --git a/docs/phantompeaks.md b/docs/phantompeaks.md
new file mode 100644
index 0000000000000000000000000000000000000000..c819230d3385d357c76fb5a4342c499d6cc06883
--- /dev/null
+++ b/docs/phantompeaks.md
@@ -0,0 +1,10 @@
+## Phantompeakqualtools
+Phantompeakqualtools plots the strand cross-correlation of aligned reads for each sample. In a strand cross-correlation plot, reads are shifted in the direction of the strand they map to by an increasing number of base pairs and the Pearson correlation between the per-position read count vectors for each strand is calculated. Two cross-correlation peaks are usually observed in a ChIP experiment, one corresponding to the read length ("phantom" peak) and one to the average fragment length of the library. The absolute and relative height of the two peaks are useful determinants of the success of a ChIP-seq experiment. A high-quality IP is characterized by a ChIP peak that is much higher than the "phantom" peak, while often very small or no such peak is seen in failed experiments.
+
+
+
+*Source: Landt SG et al, Genome Research (2012)*
+
+Normalized strand coefficient (NSC) is the normalized ratio between the fragment-length cross-correlation peak and the background cross-correlation. NSC values range from a minimum of 1 to larger positive numbers. 1.1 is the critical threshold. Datasets with NSC values much less than 1.1 (< 1.05) tend to have low signal to noise or few peaks (this could be biological eg. a factor that truly binds only a few sites in a particular tissue type OR it could be due to poor quality). ENCODE cutoff: NSC > 1.05.
+
+Relative strand correlation (RSC) is the ratio between the fragment-length peak and the read-length peak. RSC values range from 0 to larger positive values. 1 is the critical threshold. RSC values significantly lower than 1 (< 0.8) tend to have low signal to noise. The low scores can be due to failed and poor quality ChIP, low read sequence quality and hence lots of mismappings, shallow sequencing depth (significantly below saturation) or a combination of these. Like the NSC, datasets with few binding sites (< 200), which is biologically justifiable, also show low RSC scores. ENCODE cutoff: RSC > 0.8.
diff --git a/docs/references.md b/docs/references.md
new file mode 100644
index 0000000000000000000000000000000000000000..ef37be7bee620f256098d0a656b89167c20498ad
--- /dev/null
+++ b/docs/references.md
@@ -0,0 +1,60 @@
+### References
+
+1. **python**:
+ * Anaconda (Anaconda Software Distribution, [https://anaconda.com](https://anaconda.com))
+
+2. **trimgalore**:
+ * trimgalore [https://github.com/FelixKrueger/TrimGalore](https://github.com/FelixKrueger/TrimGalore)
+
+3. **cutadapt**:
+ * Marcel, M. 2011. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal 17(1):10-12. doi:[10.14806/ej.17.1.200](http://dx.doi.org/10.14806/ej.17.1.200)
+
+4. **bwa**:
+ * Li H., and R. Durbin. 2009. Fast and accurate short read alignment with Burrows-Wheeler Transform. Bioinformatics 25: 1754-60. doi:[10.1093/bioinformatics/btp324](http://dx.doi.org/10.1093/bioinformatics/btp324)
+
+5. **samtools**:
+ * Li H., B. Handsaker, A. Wysoker, T. Fennell, J. Ruan, N. Homer, G. Marth, G. Abecasis, R. Durbin, and 1000 Genome Project Data Processing Subgroup. 2009. The Sequence alignment/map (SAM) format and SAMtools. Bioinformatics 25: 2078-9. doi:[10.1093/bioinformatics/btp352](http://dx.doi.org/10.1093/bioinformatics/btp352)
+
+6. **sambamba**:
+ * Tarasov, A., A. J. Vilella, E. Cuppen, I. J. Nijman, and P. Prins. 2015 Sambamba: fast processing of NGS alignment formats. Bioinformatics 31(12): 2032-2034. doi:[10.1093/bioinformatics/btv098](http://dx.doi.org/10.1093/bioinformatics/btv098)
+
+7. **bedtools**:
+ * Quinlan, A. R., and I. M. Hall. 2010. BEDTools: a flexible suite of utilities for comparing genomic feautures. Bioinformatics 26(6): 841-842. doi:[10.1093/bioinformatics/btq033](http://dx.doi.org/10.1093/bioinformatics/btq033)
+
+8. **deeptools**:
+ * RamÃrez, F., D. P. Ryan, B. Grüning, V. Bhardwaj, F. Kilpert, A. S. Richter, S. Heyne, F. Dündar, and T. Manke. 2016. deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Research 44: W160-165. doi:[10.1093/nar/gkw257](http://dx.doi.org/10.1093/nar/gkw257)
+
+9. **phantompeakqualtools**:
+ * Landt S. G., G. K. Marinov, A. Kundaje, et al. 2012. ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia. Genome Res 9: 1813-31. doi:[10.1101/gr.136184.111](http://dx.doi.org/10.1101/gr.136184.111)
+ * Kharchenko P. K., M. Y. Tolstorukov, and P. J. Park. 2008. Design and analysis of ChIP-seq experiments for DNA-binding proteins. Nat Biotechnol 26(12): 1351-1359. doi:[10.1038/nbt.1508](https://dx.doi.org/10.1038/nbt.1508)
+
+10. **macs**:
+ * Zhang Y., T. Liu, C. A. Meyer, J. Eeckhoute, D. S. Johnson, B. E. Bernstein, C. Nusbaum, R. M. Myers, M. Brown, W. Li, and X. S. Liu. 2008. Model-based Analysis of ChIP-Seq (MACS). Genome Biol 9: R137. doi:[10.1186/gb-2008-9-9-r137](https://dx.doi.org/10.1186/gb-2008-9-9-r137)
+
+11. **UCSC(bedGraphToBigWig)**:
+ * Kent W. J., A. S. Zweig, G. Barber, A. S. Hinrichs, and D. Karolchik. BigWig and BigBed: enabling browsing of large distributed data sets. Bioinformatics 26(17): 2204-2207. doi:[10.1093/bioinformatics/btq351](https://dx.doi.org/10.1093/bioinformatics/btq351)
+
+12. **R**:
+ * R Core Team 2014. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL:[http://www.R-project.org/](http://www.R-project.org/).
+
+13. **meme**:
+ * Bailey T. L., M. Bodén, F. A. Buske, M. Frith, C. E. Grant, L. Clementi, J. Ren, W. W. Li, and W. S. Noble. 2009. MEME SUITE: tools for motif discovery and searching. Nucleic Acids Research 37: W202-W208. doi:[10.1093/nar/gkp335](https://dx.doi.org/10.1093/nar/gkp335)
+ * Machanick P., and T. L. Bailey. 2011. MEME-ChIP: motif analysis of large DNA datasets. Bioinformatics 27(12): 1696-1697. doi:[10.1093/bioinformatics/btr189](https://dx.doi.org/10.1093/bioinformatics/btr189)
+
+14. **ChIPseeker**:
+ * Yu G., L. Wang, and Q. He. 2015. ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics 31(14): 2382-2383. doi:[10.1093/bioinformatics/btv145](https://dx.doi.org/10.1093/bioinformatics/btv145)
+
+15. **DiffBind**:
+ * Stark R., and G. Brown. 2011. DiffBind: differential binding analysis of ChIP-Seq peak data. [http://bioconductor.org/packages/release/bioc/vignettes/DiffBind/inst/doc/DiffBind.pdf](http://bioconductor.org/packages/release/bioc/vignettes/DiffBind/inst/doc/DiffBind.pdf). doi:[10.18129/B9.bioc.DiffBind](https://dx.doi.org/10.18129/B9.bioc.DiffBind)
+ * Ross-Innes C. S., R. Stark, A. E. Teschendorff, K. A. Holmes, H. R. Ali, M. J. Dunning, G. D. Brown, O. Gojis, I. O. Ellis, A. R. Green, S. Ali, S. Chin, C. Palmieri, C. Caldas, and J. S. Carroll. 2012. Differential oestrogen receptor binding is associated with clinical outcome in breast cancer. Nature 481: 389-393. doi:[10.1038/nature10730](https://dx.doi.org/10.1038/nature10730)
+
+16. **MultiQc**:
+ * Ewels P., Magnusson M., Lundin S. and Käller M. 2016. MultiQC: Summarize analysis results for multiple tools and samples in a single report. Bioinformatics 32(19): 3047–3048. doi:[10.1093/bioinformatics/btw354](https://dx.doi.org/10.1093/bioinformatics/btw354)
+
+17. **BICF ChIP-seq Analysis Workflow**:
+ * Spencer D. Barnes, Holly Ruess, Jeremy A. Mathews, Beibei Chen, and Venkat S. Malladi. 2019. BICF ChIP-seq Analysis Workflow (publish_1.0.5). Zenodo. doi:[10.5281/zenodo.2648844](https://doi.org/10.5281/zenodo.2648844)
+
+18. **Nextflow**:
+ * Di Tommaso, P., Chatzou, M., Floden, E. W., Barja, P. P., Palumbo, E., and Notredame, C. 2017. Nextflow enables reproducible computational workflows. Nature biotechnology, 35(4), 316.
+
+Please cite in publications: Pipeline was developed by BICF from funding provided by **Cancer Prevention and Research Institute of Texas (RP150596)**.
diff --git a/docs/xcor_header.txt b/docs/xcor_header.txt
new file mode 100644
index 0000000000000000000000000000000000000000..c4c712ba102e3ff79142ae8536003dba04d5eb25
--- /dev/null
+++ b/docs/xcor_header.txt
@@ -0,0 +1,17 @@
+See https://github.com/crazyhottommy/phantompeakqualtools for more details
+
+COL1: Filename: tagAlign/BAM filename
+COL2: numReads: effective sequencing depth i.e. total number of mapped reads in input file
+COL3: estFragLen: comma separated strand cross-correlation peak(s) in decreasing order of correlation.
+ The top 3 local maxima locations that are within 90% of the maximum cross-correlation value are output.
+ In almost all cases, the top (first) value in the list represents the predominant fragment length.
+ If you want to keep only the top value simply run
+ sed -r 's/,[^\t]+//g' <outFile> > <newOutFile>
+COL4: corr_estFragLen: comma separated strand cross-correlation value(s) in decreasing order (col2 follows the same order)
+COL5: phantomPeak: Read length/phantom peak strand shift
+COL6: corr_phantomPeak: Correlation value at phantom peak
+COL7: argmin_corr: strand shift at which cross-correlation is lowest
+COL8: min_corr: minimum value of cross-correlation
+COL9: Normalized strand cross-correlation coefficient (NSC) = COL4 / COL8
+COL10: Relative strand cross-correlation coefficient (RSC) = (COL4 - COL8) / (COL6 - COL8)
+COL11: QualityTag: Quality tag based on thresholded RSC (codes: -2:veryLow,-1:Low,0:Medium,1:High,2:veryHigh)
diff --git a/test_data/A_1.bedse.gz b/test_data/A_1.bedse.gz
new file mode 100644
index 0000000000000000000000000000000000000000..4bbf615c84c514040c206bde70190b83fac1003e
Binary files /dev/null and b/test_data/A_1.bedse.gz differ
diff --git a/test_data/A_1.tagAlign.gz b/test_data/A_1.tagAlign.gz
new file mode 100644
index 0000000000000000000000000000000000000000..b5cc068550fc3c850c014935e21586a89b5e8118
Binary files /dev/null and b/test_data/A_1.tagAlign.gz differ
diff --git a/test_data/B_1.bedse.gz b/test_data/B_1.bedse.gz
new file mode 100644
index 0000000000000000000000000000000000000000..4bbf615c84c514040c206bde70190b83fac1003e
Binary files /dev/null and b/test_data/B_1.bedse.gz differ
diff --git a/test_data/B_1.tagAlign.gz b/test_data/B_1.tagAlign.gz
new file mode 100644
index 0000000000000000000000000000000000000000..b5cc068550fc3c850c014935e21586a89b5e8118
Binary files /dev/null and b/test_data/B_1.tagAlign.gz differ
diff --git a/workflow/conf/bicf_logo.png b/workflow/conf/bicf_logo.png
new file mode 100644
index 0000000000000000000000000000000000000000..0d8015590c5a94f92c39ec2470bd02baa3d09077
Binary files /dev/null and b/workflow/conf/bicf_logo.png differ
diff --git a/workflow/conf/biohpc.config b/workflow/conf/biohpc.config
index 2c0f1c31a9c46b5f9f45699f3c1eb506d6c62e10..037d1e43e37c7daeb21489a00dff85701c8f8fd0 100644
--- a/workflow/conf/biohpc.config
+++ b/workflow/conf/biohpc.config
@@ -41,9 +41,13 @@ process {
executor = 'local'
}
withName: callPeaksMACS {
- module = ['python/3.6.1-2-anaconda', 'macs/2.1.0-20151222', 'phantompeakqualtools/1.2', 'UCSC_userApps/v317', 'bedtools/2.26.0']
+ module = ['python/3.6.1-2-anaconda', 'macs/2.1.0-20151222', 'UCSC_userApps/v317', 'bedtools/2.26.0', 'phantompeakqualtools/1.2']
queue = '128GB,256GB,256GBv1'
}
+ withName: plotProfile {
+ module = ['deeptools/2.5.0.1']
+ cpus = 32
+ }
withName: consensusPeaks {
module = ['python/3.6.1-2-anaconda', 'bedtools/2.26.0']
executor = 'local'
@@ -60,29 +64,38 @@ process {
module = ['python/3.6.1-2-anaconda', 'meme/4.11.1-gcc-openmpi', 'bedtools/2.26.0']
cpus = 32
}
-
+ withName: multiqcReport {
+ module = ['python/3.6.1-2-anaconda', 'pandoc/2.7', 'singularity/3.0.2']
+ executor = 'local'
+ }
}
params {
// Reference file paths on BioHPC
genomes {
'GRCh38' {
- bwa = '/project/shared/bicf_workflow_ref/GRCh38'
+ bwa = '/project/shared/bicf_workflow_ref/human/GRCh38'
genomesize = 'hs'
- chromsizes = '/project/shared/bicf_workflow_ref/GRCh38/genomefile.txt'
- fasta = '/project/shared/bicf_workflow_ref/GRCh38/genome.fa'
+ chromsizes = '/project/shared/bicf_workflow_ref/human/GRCh38/genomefile.txt'
+ fasta = '/project/shared/bicf_workflow_ref/human/GRCh38/genome.fa'
+ gtf = '/project/shared/bicf_workflow_ref/human/GRCh38/gencode.v25.chr_patch_hapl_scaff.annotation.gtf'
+ geneNames = '/project/shared/bicf_workflow_ref/human/GRCh38/genenames.txt'
}
'GRCh37' {
- bwa = '/project/shared/bicf_workflow_ref/GRCh37'
+ bwa = '/project/shared/bicf_workflow_ref/human/GRCh37'
genomesize = 'hs'
- chromsizes = '/project/shared/bicf_workflow_ref/GRCh37/genomefile.txt'
- fasta = '/project/shared/bicf_workflow_ref/GRCh37/genome.fa'
+ chromsizes = '/project/shared/bicf_workflow_ref/human/GRCh37/genomefile.txt'
+ fasta = '/project/shared/bicf_workflow_ref/human/GRCh37/genome.fa'
+ gtf = '/project/shared/bicf_workflow_ref/human/GRCh37/gencode.v19.chr_patch_hapl_scaff.annotation.gtf'
+ geneNames = '/project/shared/bicf_workflow_ref/human/GRCh37/genenames.txt'
}
'GRCm38' {
- bwa = '/project/shared/bicf_workflow_ref/GRCm38'
+ bwa = '/project/shared/bicf_workflow_ref/mouse/GRCm38'
genomesize = 'mm'
- chromsizes = '/project/shared/bicf_workflow_ref/GRCm38/genomefile.txt'
- fasta = '/project/shared/bicf_workflow_ref/GRCm38/genome.fa'
+ chromsizes = '/project/shared/bicf_workflow_ref/mouse/GRCm38/genomefile.txt'
+ fasta = '/project/shared/bicf_workflow_ref/mouse/GRCm38/genome.fa'
+ gtf = '/project/shared/bicf_workflow_ref/mouse/GRCm38/gencode.vM20.annotation.gtf'
+ geneNames = '/project/shared/bicf_workflow_ref/mouse/GRCm38/genenames.txt'
}
}
}
diff --git a/workflow/conf/multiqc_config.yaml b/workflow/conf/multiqc_config.yaml
new file mode 100644
index 0000000000000000000000000000000000000000..244d1e162cb9e1e948e8aa3d77a104a716b5602c
--- /dev/null
+++ b/workflow/conf/multiqc_config.yaml
@@ -0,0 +1,109 @@
+# Custom Logo
+custom_logo: 'bicf_logo.png'
+custom_logo_url: 'https://www.utsouthwestern.edu/labs/bioinformatics/'
+custom_logo_title: 'Bioinformatics Core Facility'
+
+report_header_info:
+ - Contact E-mail: 'bicf@utsouthwestern.edu'
+ - Application Type: 'ChIP-seq'
+ - Department: 'Bioinformatic Core Facility, Department of Bioinformatics'
+ - Contributors and Licensing: 'See https://doi.org/10.5281/zenodo.2648844'
+
+
+# Title to use for the report.
+title: BICF ChIP-seq Analysis Report
+
+report_comment: >
+ This report has been generated by the <a href="https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/" target="_blank">BICF/chipseq_analysis</a>
+ pipeline.
+
+extra_fn_clean_exts:
+ - 'pbc.qc'
+ - 'cc.qc'
+
+fn_ignore_files:
+ - '*dedup.flagstat.qc'
+
+custom_data:
+ library_complexity:
+ file_format: 'tsv'
+ id: 'library_complexity'
+ contents: 'TotalReadPairs DistinctReadPairs OneReadPair TwoReadPairs NRF PBC1 PBC2'
+ section_name: 'Library complexity'
+ plot_type: 'generalstats'
+ pconfig:
+ TotalReadPairs:
+ decimalPlaces: 0
+ shared_key: read_count
+ DistinctReadPairs:
+ decimalPlaces: 0
+ shared_key: read_count
+ NRF:
+ decimalPlaces: 2
+ PBC1:
+ decimalPlaces: 2
+ PBC2:
+ decimalPlaces: 2
+
+sp:
+ phantompeakqualtools/out:
+ fn: '*cc.qc'
+ library_complexity:
+ fn: '*pbc.qc'
+ macs2:
+ fn: '*_peaks.xls'
+
+
+report_section_order:
+ cutadapt:
+ order: -1000
+ Samtools:
+ order: -1100
+ Software_Versions:
+ order: -1200
+ Software_References:
+ order: -1300
+
+table_columns_placement:
+ library_complexity:
+ TotalReadPairs: 1100
+ DistinctReadPairs: 1200
+ NRF: 1300
+ PBC1: 1400
+ PBC2: 1500
+ phantompeakqualtools:
+ Estimated_Fragment_Length_bp: 1600
+ NSC: 1700
+ RSC: 1800
+
+table_columns_visible:
+ cutadapt:
+ percent_trimmed: False
+ library_complexity:
+ OneReadPair: False
+ TwoReadPairs: False
+
+table_cond_formatting_rules:
+ library_complexity_mqc-generalstats-library_complexity-NRF:
+ pass:
+ - gt: 0.8
+ warn:
+ - lt: 0.8
+ fail:
+ - lt: 0.5
+ library_complexity_mqc-generalstats-library_complexity-PBC1:
+ pass:
+ - gt: 0.8
+ warn:
+ - lt: 0.8
+ fail:
+ - lt: 0.5
+ library_complexity_mqc-generalstats-library_complexity-PBC2:
+ pass:
+ - gt: 3
+ warn:
+ - lt: 3
+ fail:
+ - lt: 1
+
+thousandsSep_format: ''
diff --git a/workflow/main.nf b/workflow/main.nf
index 6051ff7cf3d7409db201634a14814a220bab2cf2..bf87941930e7a3737c8610703d224887c113a4d1 100644
--- a/workflow/main.nf
+++ b/workflow/main.nf
@@ -1,38 +1,71 @@
#!/usr/bin/env nextflow
+/*
+
+BICF ChIP-seq Analysis Workflow
+#### Homepage / Documentation
+https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/
+Licensed under MIT (https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/LICENSE.md)
+
+
+*/
+
// Path to an input file, or a pattern for multiple inputs
// Note - $baseDir is the location of this workflow file main.nf
// Define Input variables
params.reads = "$baseDir/../test_data/*.fastq.gz"
-params.pairedEnd = 'false'
+params.pairedEnd = false
params.designFile = "$baseDir/../test_data/design_ENCSR238SGC_SE.txt"
+params.genome = 'GRCm38'
params.cutoffRatio = 1.2
params.outDir= "$baseDir/output"
params.extendReadsLen = 100
params.topPeakCount = 600
-params.astrocyte = 'false'
+params.astrocyte = false
+params.skipDiff = false
+params.skipMotif = false
+params.skipPlotProfile = false
+params.references = "$baseDir/../docs/references.md"
+params.multiqc = "$baseDir/conf/multiqc_config.yaml"
// Assign variables if astrocyte
-params.genome = 'GRCm38'
-if (params.astrocyte == 'false') {
- params.bwaIndex = params.genome ? params.genomes[ params.genome ].bwa ?: false : false
- params.genomeSize = params.genome ? params.genomes[ params.genome ].genomesize ?: false : false
- params.chromSizes = params.genome ? params.genomes[ params.genome ].chromsizes ?: false : false
- params.fasta = params.genome ? params.genomes[ params.genome ].fasta ?: false : false
-} else if (params.astrocyte == 'true') {
+if (params.astrocyte) {
+ print("Running under astrocyte")
referenceLocation = "/project/shared/bicf_workflow_ref"
- params.bwaIndex = "$referenceLocation/$genome"
- params.chromSizes = "$referenceLocation/$genome/genomefile.txt"
- params.fasta = "$referenceLocation/$genome/genome.fa.txt"
- if (params.genome == 'GRCh37' || params.genome == 'GRCh38') {
- params.chromSizes = 'hs'
- } else if (params.chromSizes == 'GRCm38') {
- params.chromSizes = 'mm'
+ if (params.genome == 'GRCh37') {
+ params.bwaIndex = "$referenceLocation/human/$params.genome"
+ params.chromSizes = "$referenceLocation/human/$params.genome/genomefile.txt"
+ params.fasta = "$referenceLocation/human/$params.genome/genome.fa"
+ params.gtf = "$referenceLocation/human/$params.genome/gencode.v19.chr_patch_hapl_scaff.annotation.gtf"
+ params.geneNames = "$referenceLocation/human/$params.genome/genenames.txt"
+ params.genomeSize = 'hs'
+ } else if (params.genome == 'GRCm38') {
+ params.bwaIndex = "$referenceLocation/mouse/$params.genome"
+ params.chromSizes = "$referenceLocation/mouse/$params.genome/genomefile.txt"
+ params.fasta = "$referenceLocation/mouse/$params.genome/genome.fa"
+ params.gtf = "$referenceLocation/mouse/$params.genome/gencode.vM20.annotation.gtf"
+ params.geneNames = "$referenceLocation/mouse/$params.genome/genenames.txt"
+ params.genomeSize = 'mm'
+ } else if (params.genome == 'GRCh38') {
+ params.bwaIndex = "$referenceLocation/human/$params.genome"
+ params.chromSizes = "$referenceLocation/human/$params.genome/genomefile.txt"
+ params.fasta = "$referenceLocation/human/$params.genome/genome.fa"
+ params.gtf = "$referenceLocation/human/$params.genome/gencode.v25.chr_patch_hapl_scaff.annotation.gtf"
+ params.geneNames = "$referenceLocation/human/$params.genome/genenames.txt"
+ params.genomeSize = 'hs'
}
+} else {
+ params.bwaIndex = params.genome ? params.genomes[ params.genome ].bwa ?: false : false
+ params.genomeSize = params.genome ? params.genomes[ params.genome ].genomesize ?: false : false
+ params.chromSizes = params.genome ? params.genomes[ params.genome ].chromsizes ?: false : false
+ params.fasta = params.genome ? params.genomes[ params.genome ].fasta ?: false : false
+ params.gtf = params.genome ? params.genomes[ params.genome ].gtf ?: false : false
+ params.geneNames = params.genome ? params.genomes[ params.genome ].geneNames ?: false : false
}
+
// Check inputs
if( params.bwaIndex ){
bwaIndex = Channel
@@ -50,7 +83,8 @@ readsList = Channel
.collectFile( name: 'fileList.tsv', newLine: true )
// Define regular variables
-designFile = params.designFile
+pairedEnd = params.pairedEnd
+designFile = Channel.fromPath(params.designFile)
genomeSize = params.genomeSize
genome = params.genome
chromSizes = params.chromSizes
@@ -59,12 +93,14 @@ cutoffRatio = params.cutoffRatio
outDir = params.outDir
extendReadsLen = params.extendReadsLen
topPeakCount = params.topPeakCount
-
-if (params.pairedEnd == 'false'){
- pairedEnd = false
-} else {
- pairedEnd = true
-}
+skipDiff = params.skipDiff
+skipMotif = params.skipMotif
+skipPlotProfile = params.skipPlotProfile
+references = params.references
+multiqc = params.multiqc
+gtfFile_plotProfile = Channel.fromPath(params.gtf)
+gtfFile_annotPeaks = Channel.fromPath(params.gtf)
+geneNames = Channel.fromPath(params.geneNames)
// Check design file for errors
process checkDesignFile {
@@ -73,7 +109,7 @@ process checkDesignFile {
input:
- designFile
+ file designFile
file readsList
output:
@@ -84,11 +120,13 @@ process checkDesignFile {
if (pairedEnd) {
"""
+ module load python/3.6.1-2-anaconda
python3 $baseDir/scripts/check_design.py -d $designFile -f $readsList -p
"""
}
else {
"""
+ module load python/3.6.1-2-anaconda
python $baseDir/scripts/check_design.py -d $designFile -f $readsList
"""
}
@@ -110,7 +148,7 @@ rawReads = designFilePaths
process trimReads {
tag "$sampleId-$replicate"
- publishDir "$outDir/${task.process}", mode: 'copy'
+ publishDir "$outDir/${task.process}/${sampleId}", mode: 'copy'
input:
@@ -119,17 +157,22 @@ process trimReads {
output:
set sampleId, file('*.fq.gz'), experimentId, biosample, factor, treatment, replicate, controlId into trimmedReads
- file('*trimming_report.txt') into trimgalore_results
+ file('*trimming_report.txt') into trimgaloreResults
+ file('version_*.txt') into trimReadsVersions
script:
if (pairedEnd) {
"""
+ module load python/3.6.1-2-anaconda
+ module load trimgalore/0.4.1
python3 $baseDir/scripts/trim_reads.py -f ${reads[0]} ${reads[1]} -s $sampleId -p
"""
}
else {
"""
+ module load python/3.6.1-2-anaconda
+ module load trimgalore/0.4.1
python3 $baseDir/scripts/trim_reads.py -f ${reads[0]} -s $sampleId
"""
}
@@ -139,8 +182,9 @@ process trimReads {
// Align trimmed reads using bwa
process alignReads {
+ queue '128GB,256GB,256GBv1'
tag "$sampleId-$replicate"
- publishDir "$outDir/${task.process}", mode: 'copy'
+ publishDir "$outDir/${task.process}/${sampleId}", mode: 'copy'
input:
@@ -151,16 +195,23 @@ process alignReads {
set sampleId, file('*.bam'), experimentId, biosample, factor, treatment, replicate, controlId into mappedReads
file '*.flagstat.qc' into mappedReadsStats
+ file('version_*.txt') into alignReadsVersions
script:
if (pairedEnd) {
"""
+ module load python/3.6.1-2-anaconda
+ module load bwa/intel/0.7.12
+ module load samtools/1.6
python3 $baseDir/scripts/map_reads.py -f ${reads[0]} ${reads[1]} -r ${index}/genome.fa -s $sampleId -p
"""
}
else {
"""
+ module load python/3.6.1-2-anaconda
+ module load bwa/intel/0.7.12
+ module load samtools/1.6
python3 $baseDir/scripts/map_reads.py -f $reads -r ${index}/genome.fa -s $sampleId
"""
}
@@ -170,8 +221,9 @@ process alignReads {
// Dedup reads using sambamba
process filterReads {
+ queue '128GB,256GB,256GBv1'
tag "$sampleId-$replicate"
- publishDir "$outDir/${task.process}", mode: 'copy'
+ publishDir "$outDir/${task.process}/${sampleId}", mode: 'copy'
input:
@@ -184,16 +236,25 @@ process filterReads {
file '*.flagstat.qc' into dedupReadsStats
file '*.pbc.qc' into dedupReadsComplexity
file '*.dedup.qc' into dupReads
+ file('version_*.txt') into filterReadsVersions
script:
if (pairedEnd) {
"""
+ module load python/3.6.1-2-anaconda
+ module load samtools/1.6
+ module load sambamba/0.6.6
+ module load bedtools/2.26.0
python3 $baseDir/scripts/map_qc.py -b $mapped -p
"""
}
else {
"""
+ module load python/3.6.1-2-anaconda
+ module load samtools/1.6
+ module load sambamba/0.6.6
+ module load bedtools/2.26.0
python3 $baseDir/scripts/map_qc.py -b $mapped
"""
}
@@ -210,6 +271,7 @@ dedupReads
// Quality Metrics using deeptools
process experimentQC {
+ queue '128GB,256GB,256GBv1'
publishDir "$outDir/${task.process}", mode: 'copy'
input:
@@ -218,11 +280,14 @@ process experimentQC {
output:
- file '*.{png,npz}' into deepToolsStats
+ file '*.{pdf,npz}' into experimentQCStats
+ file('version_*.txt') into experimentQCVersions
script:
"""
+ module load python/3.6.1-2-anaconda
+ module load deeptools/2.5.0.1
python3 $baseDir/scripts/experiment_qc.py -d $dedupDesign -e $extendReadsLen
"""
@@ -231,8 +296,9 @@ process experimentQC {
// Convert reads to bam
process convertReads {
+ queue '128GB,256GB,256GBv1'
tag "$sampleId-$replicate"
- publishDir "$outDir/${task.process}", mode: 'copy'
+ publishDir "$outDir/${task.process}/${sampleId}", mode: 'copy'
input:
@@ -241,16 +307,23 @@ process convertReads {
output:
set sampleId, file('*.tagAlign.gz'), file('*.bed{pe,se}.gz'), experimentId, biosample, factor, treatment, replicate, controlId into tagReads
+ file('version_*.txt') into convertReadsVersions
script:
if (pairedEnd) {
"""
+ module load python/3.6.1-2-anaconda
+ module load samtools/1.6
+ module load bedtools/2.26.0
python3 $baseDir/scripts/convert_reads.py -b $deduped -p
"""
}
else {
"""
+ module load python/3.6.1-2-anaconda
+ module load samtools/1.6
+ module load bedtools/2.26.0
python3 $baseDir/scripts/convert_reads.py -b $deduped
"""
}
@@ -261,7 +334,7 @@ process convertReads {
process crossReads {
tag "$sampleId-$replicate"
- publishDir "$outDir/${task.process}", mode: 'copy'
+ publishDir "$outDir/${task.process}/${sampleId}", mode: 'copy'
input:
@@ -270,17 +343,22 @@ process crossReads {
output:
set sampleId, seTagAlign, tagAlign, file('*.cc.qc'), experimentId, biosample, factor, treatment, replicate, controlId into xcorReads
- set file('*.cc.qc'), file('*.cc.plot.pdf') into xcorReadsStats
+ set file('*.cc.qc'), file('*.cc.plot.pdf') into crossReadsStats
+ file('version_*.txt') into crossReadsVersions
script:
if (pairedEnd) {
"""
+ module load python/3.6.1-2-anaconda
+ module load phantompeakqualtools/1.2
python3 $baseDir/scripts/xcor.py -t $seTagAlign -p
"""
}
else {
"""
+ module load python/3.6.1-2-anaconda
+ module load phantompeakqualtools/1.2
python3 $baseDir/scripts/xcor.py -t $seTagAlign
"""
}
@@ -309,6 +387,7 @@ process defineExpDesignFiles {
script:
"""
+ module load python/3.6.1-2-anaconda
python3 $baseDir/scripts/experiment_design.py -d $xcorDesign
"""
@@ -334,11 +413,13 @@ process poolAndPsuedoReads {
if (pairedEnd) {
"""
+ module load python/3.6.1-2-anaconda
python3 $baseDir/scripts/pool_and_psuedoreplicate.py -d $experimentObjs -c $cutoffRatio -p
"""
}
else {
"""
+ module load python/3.6.1-2-anaconda
python3 $baseDir/scripts/pool_and_psuedoreplicate.py -d $experimentObjs -c $cutoffRatio
"""
}
@@ -354,7 +435,7 @@ experimentRows = experimentPoolObjs
process callPeaksMACS {
tag "$sampleId-$replicate"
- publishDir "$outDir/${task.process}", mode: 'copy'
+ publishDir "$outDir/${task.process}/${experimentId}/${replicate}", mode: 'copy'
input:
set sampleId, tagAlign, xcor, experimentId, biosample, factor, treatment, replicate, controlId, controlTagAlign from experimentRows
@@ -362,17 +443,29 @@ process callPeaksMACS {
output:
set sampleId, file('*.narrowPeak'), file('*.fc_signal.bw'), file('*.pvalue_signal.bw'), experimentId, biosample, factor, treatment, replicate, controlId into experimentPeaks
- file '*.xls' into summit
+ file '*.xls' into callPeaksMACSsummit
+ file('version_*.txt') into callPeaksMACSVersions
+ file("*.fc_signal.bw") into bigwigs
script:
if (pairedEnd) {
"""
+ module load python/3.6.1-2-anaconda
+ module load macs/2.1.0-20151222
+ module load UCSC_userApps/v317
+ module load bedtools/2.26.0
+ module load phantompeakqualtools/1.2
python3 $baseDir/scripts/call_peaks_macs.py -t $tagAlign -x $xcor -c $controlTagAlign -s $sampleId -g $genomeSize -z $chromSizes -p
"""
}
else {
"""
+ module load python/3.6.1-2-anaconda
+ module load macs/2.1.0-20151222
+ module load UCSC_userApps/v317
+ module load bedtools/2.26.0
+ module load phantompeakqualtools/1.2
python3 $baseDir/scripts/call_peaks_macs.py -t $tagAlign -x $xcor -c $controlTagAlign -s $sampleId -g $genomeSize -z $chromSizes
"""
}
@@ -385,6 +478,30 @@ peaksDesign = experimentPeaks
"$sampleId\t$peak\t$fcSignal\t$pvalueSignal\t$experimentId\t$biosample\t$factor\t$treatment\t$replicate\t$controlId\n"}
.collectFile(name:'design_peak.tsv', seed:"sample_id\tpeaks\tfc_signal\tpvalue_signal\texperiment_id\tbiosample\tfactor\ttreatment\treplicate\tcontrol_id\n", storeDir:"$outDir/design")
+//plotProfile
+process plotProfile {
+ publishDir "$outDir/experimentQC", mode: 'copy'
+
+ input:
+
+ file ("*.pooled.fc_signal.bw") from bigwigs.collect()
+ file gtf from gtfFile_plotProfile
+
+ output:
+
+ file '*.{png,gz}' into plotProfile
+
+ when:
+
+ !skipPlotProfile
+
+ script:
+ """
+ module load deeptools/2.5.0.1
+ bash $baseDir/scripts/plotProfile.sh
+ """
+}
+
// Calculate Consensus Peaks
process consensusPeaks {
@@ -403,10 +520,13 @@ process consensusPeaks {
file 'design_diffPeaks.csv' into designDiffPeaks
file 'design_annotatePeaks.tsv' into designAnnotatePeaks, designMotifSearch
file 'unique_experiments.csv' into uniqueExperiments
+ file('version_*.txt') into consensusPeaksVersions
script:
"""
+ module load python/3.6.1-2-anaconda
+ module load bedtools/2.26.0
python3 $baseDir/scripts/overlap_peaks.py -d $peaksDesign -f $preDiffDesign
"""
@@ -415,28 +535,32 @@ process consensusPeaks {
// Annotate Peaks
process peakAnnotation {
- publishDir "$baseDir/output/${task.process}", mode: 'copy'
+ publishDir "$outDir/${task.process}", mode: 'copy'
input:
file designAnnotatePeaks
+ file gtf from gtfFile_annotPeaks
+ file geneNames
output:
file "*chipseeker*" into peakAnnotation
+ file('version_*.txt') into peakAnnotationVersions
script:
"""
- Rscript $baseDir/scripts/annotate_peaks.R $designAnnotatePeaks $genome
+ module load R/3.3.2-gccmkl
+ Rscript $baseDir/scripts/annotate_peaks.R $designAnnotatePeaks $gtf $geneNames
"""
}
-// Motif Search Peaks
+// Motif Search Peaks
process motifSearch {
- publishDir "$baseDir/output/${task.process}", mode: 'copy'
+ publishDir "$outDir/${task.process}", mode: 'copy'
input:
@@ -446,10 +570,18 @@ process motifSearch {
file "*memechip" into motifSearch
file "*narrowPeak" into filteredPeaks
+ file('version_*.txt') into motifSearchVersions
+
+ when:
+
+ !skipMotif
script:
"""
+ module load python/3.6.1-2-anaconda
+ module load meme/4.11.1-gcc-openmpi
+ module load bedtools/2.26.0
python3 $baseDir/scripts/motif_search.py -d $designMotifSearch -g $fasta -p $topPeakCount
"""
}
@@ -461,7 +593,7 @@ uniqueExperimentsList = uniqueExperiments
// Calculate Differential Binding Activity
process diffPeaks {
- publishDir "$baseDir/output/${task.process}", mode: 'copy'
+ publishDir "$outDir/${task.process}", mode: 'copy'
input:
@@ -474,12 +606,58 @@ process diffPeaks {
file '*_diffbind.csv' into diffPeaksCounts
file '*.pdf' into diffPeaksStats
file 'normcount_peaksets.txt' into normCountPeaks
+ file('version_*.txt') into diffPeaksVersions
when:
- noUniqueExperiments > 1
+
+ noUniqueExperiments > 1 && !skipDiff
script:
+
"""
+ module load python/3.6.1-2-anaconda
+ module load R/3.3.2-gccmkl
Rscript $baseDir/scripts/diff_peaks.R $designDiffPeaks
"""
}
+
+// Generate Multiqc Report, gerernate Software Versions and references
+process multiqcReport {
+ publishDir "$outDir/${task.process}", mode: 'copy'
+
+ input:
+
+ file ('trimReads_vf/*') from trimReadsVersions.first()
+ file ('alignReads_vf/*') from alignReadsVersions.first()
+ file ('filterReads_vf/*') from filterReadsVersions.first()
+ file ('convertReads_vf/*') from convertReadsVersions.first()
+ file ('crossReads_vf/*') from crossReadsVersions.first()
+ file ('callPeaksMACS_vf/*') from callPeaksMACSVersions.first()
+ file ('consensusPeaks_vf/*') from consensusPeaksVersions.first()
+ file ('peakAnnotation_vf/*') from peakAnnotationVersions.first()
+ file ('motifSearch_vf/*') from motifSearchVersions.first().ifEmpty()
+ file ('diffPeaks_vf/*') from diffPeaksVersions.first().ifEmpty()
+ file ('experimentQC_vf/*') from experimentQCVersions.first()
+ file ('trimReads/*') from trimgaloreResults.collect()
+ file ('alignReads/*') from mappedReadsStats.collect()
+ file ('filterReads/*') from dedupReadsComplexity.collect()
+ file ('crossReads/*') from crossReadsStats.collect()
+
+ output:
+
+ file('software_versions_mqc.yaml') into softwareVersions
+ file('software_references_mqc.yaml') into softwareReferences
+ file "multiqc_report.html" into multiqcReport
+ file "*_data" into multiqcData
+
+ script:
+
+ """
+ echo $workflow.nextflow.version > version_nextflow.txt
+ singularity exec /project/shared/bicf_workflow_ref/singularity_images/multiqc.sif multiqc --version > version_multiqc.txt
+ python --version &> version_python.txt
+ python3 $baseDir/scripts/generate_references.py -r $references -o software_references
+ python3 $baseDir/scripts/generate_versions.py -o software_versions
+ singularity exec /project/shared/bicf_workflow_ref/singularity_images/multiqc.sif multiqc -c $multiqc .
+ """
+}
diff --git a/workflow/nextflow.config b/workflow/nextflow.config
index f9fbe846df7cbce559c0e8ec87097921792640a2..e1dd50f4e24db974e055d648003c901258f2c8c2 100644
--- a/workflow/nextflow.config
+++ b/workflow/nextflow.config
@@ -7,8 +7,8 @@ profiles {
manifest {
name = 'chipseq_analysis'
description = 'BICF ChIP-seq Analysis Workflow.'
- homePage = 'https://github.com/nf-core/rnaseq'
- version = '1.0.0'
+ homePage = 'https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis'
+ version = '1.0.6'
mainScript = 'main.nf'
nextflowVersion = '>=0.31.0'
}
diff --git a/workflow/scripts/annotate_peaks.R b/workflow/scripts/annotate_peaks.R
index 4bc8d79296150f64c1a1eed2ecba0b5e33b597d1..853f7aa677a796b6893762b6679573f2049766d3 100644
--- a/workflow/scripts/annotate_peaks.R
+++ b/workflow/scripts/annotate_peaks.R
@@ -1,39 +1,37 @@
#!/bin/Rscript
-# Load libraries
-library("ChIPseeker")
+#*
+#* --------------------------------------------------------------------------
+#* Licensed under MIT (https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/LICENSE.md)
+#* --------------------------------------------------------------------------
+#*
-# Currently mouse or human
-
-library("TxDb.Hsapiens.UCSC.hg19.knownGene")
-library("TxDb.Mmusculus.UCSC.mm10.knownGene")
-library("TxDb.Hsapiens.UCSC.hg38.knownGene")
-library("org.Hs.eg.db")
-library("org.Mm.eg.db")
+#Currently Human or Mouse
+# Load libraries
+library("ChIPseeker")
+library(GenomicFeatures)
# Create parser object
args <- commandArgs(trailingOnly=TRUE)
# Check input args
-if (length(args) != 2) {
- stop("Usage: annotate_peaks.R annotate_design.tsv genome_assembly", call.=FALSE)
+if (length(args) != 3) {
+ stop("Usage: annotate_peaks.R annotate_design.tsv gtf geneNames", call.=FALSE)
}
design_file <- args[1]
-genome_assembly <- args[2]
+gtf <- args[2]
+geneNames <- args[3]
# Load UCSC Known Genes
-if(genome_assembly=='GRCh37') {
- txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
- annodb <- 'org.Hs.eg.db'
-} else if(genome_assembly=='GRCm38') {
- txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
- annodb <- 'org.Mm.eg.db'
-} else if(genome_assembly=='GRCh38') {
- txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
- annodb <- 'org.Hs.eg.db'
-}
+txdb <- makeTxDbFromGFF(gtf)
+sym <- read.table(geneNames, header=T, sep='\t') [,4:5]
+
+# Output version of ChIPseeker
+chipseeker_version = packageVersion('ChIPseeker')
+write.table(paste("Version", chipseeker_version), file = "version_ChIPseeker.txt", sep = "\t",
+ row.names = FALSE, col.names = FALSE)
# Load design file
design <- read.csv(design_file, sep ='\t')
@@ -43,18 +41,18 @@ names(files) <- design$Condition
# Granges of files
peaks <- lapply(files, readPeakFile, as = "GRanges", header = FALSE)
-peakAnnoList <- lapply(peaks, annotatePeak, TxDb=txdb, annoDb=annodb, tssRegion=c(-3000, 3000), verbose=FALSE)
+peakAnnoList <- lapply(peaks, annotatePeak, TxDb=txdb, tssRegion=c(-3000, 3000), verbose=FALSE)
column_names <- c("chr", "start", "end", "width", "strand_1", "name", "score", "strand", "signalValue",
"pValue", "qValue", "peak", "annotation", "geneChr", "geneStart", "geneEnd",
- "geneLength" ,"geneStrand", "geneId", "transcriptId", "distanceToTSS",
- "ENSEMBL", "symbol", "geneName")
+ "geneLength" ,"geneStrand", "geneId", "transcriptId", "distanceToTSS", "symbol")
for(index in c(1:length(peakAnnoList))) {
filename <- paste(names(peaks)[index], ".chipseeker_annotation.tsv", sep="")
df <- as.data.frame(peakAnnoList[[index]])
- colnames(df) <- column_names
- write.table(df[ , !(names(df) %in% c('strand_1'))], filename, sep="\t" ,quote=F, row.names=F)
+ df_final <- merge(df, sym, by.x="geneId", by.y="ensembl", all.x=T)
+ colnames(df_final) <- column_names
+ write.table(df_final[ , !(names(df_final) %in% c('strand_1'))], filename, sep="\t" ,quote=F, row.names=F)
# Draw individual plots
diff --git a/workflow/scripts/call_peaks_macs.py b/workflow/scripts/call_peaks_macs.py
index a68a59a1117786cb87616f0c66483286a3ff66e7..6981262cb637b5e5f5418261f2f6171ad8f560d4 100644
--- a/workflow/scripts/call_peaks_macs.py
+++ b/workflow/scripts/call_peaks_macs.py
@@ -1,10 +1,17 @@
#!/usr/bin/env python3
+#
+# * --------------------------------------------------------------------------
+# * Licensed under MIT (https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/LICENSE.md)
+# * --------------------------------------------------------------------------
+#
+
'''Generate peaks from data.'''
import os
import argparse
import shutil
+import subprocess
import logging
import utils
from xcor import xcor as calculate_xcor
@@ -69,16 +76,19 @@ def check_tools():
logger.info('Checking for required libraries and components on this system')
- r_path = shutil.which("R")
- if r_path:
- logger.info('Found R: %s', r_path)
- else:
- logger.error('Missing R')
- raise Exception('Missing R')
-
macs_path = shutil.which("macs2")
- if r_path:
+ if macs_path:
logger.info('Found MACS2: %s', macs_path)
+
+ # Get Version
+ macs_version_command = "macs2 --version"
+ macs_version = subprocess.check_output(macs_version_command, shell=True, stderr=subprocess.STDOUT)
+
+ # Write to file
+ macs_file = open("version_macs.txt", "wb")
+ macs_file.write(macs_version)
+ macs_file.close()
+
else:
logger.error('Missing MACS2')
raise Exception('Missing MACS2')
@@ -86,6 +96,18 @@ def check_tools():
bg_bw_path = shutil.which("bedGraphToBigWig")
if bg_bw_path:
logger.info('Found bedGraphToBigWig: %s', bg_bw_path)
+
+ # Get Version
+ bg_bw_version_command = "bedGraphToBigWig"
+ try:
+ subprocess.check_output(bg_bw_version_command, shell=True, stderr=subprocess.STDOUT)
+ except subprocess.CalledProcessError as e:
+ bg_bw_version = e.output
+
+ # Write to file
+ bg_bw_file = open("version_bedGraphToBigWig.txt", "wb")
+ bg_bw_file.write(bg_bw_version)
+ bg_bw_file.close()
else:
logger.error('Missing bedGraphToBigWig')
raise Exception('Missing bedGraphToBigWig')
@@ -93,6 +115,16 @@ def check_tools():
bedtools_path = shutil.which("bedtools")
if bedtools_path:
logger.info('Found bedtools: %s', bedtools_path)
+
+ # Get Version
+ bedtools_version_command = "bedtools --version"
+ bedtools_version = subprocess.check_output(bedtools_version_command, shell=True)
+
+ # Write to file
+ bedtools_file = open("version_bedtools.txt", "wb")
+ bedtools_file.write(bedtools_version)
+ bedtools_file.close()
+
else:
logger.error('Missing bedtools')
raise Exception('Missing bedtools')
@@ -111,7 +143,6 @@ def call_peaks_macs(experiment, xcor, control, prefix, genome_size, chrom_sizes)
if int(fragment_length) < 1:
fragment_length = "200"
-
# Generate narrow peaks and preliminary signal tracks
command = 'macs2 callpeak ' + \
@@ -131,7 +162,6 @@ def call_peaks_macs(experiment, xcor, control, prefix, genome_size, chrom_sizes)
narrowpeak_fn = '%s.narrowPeak' % (prefix)
clipped_narrowpeak_fn = 'clipped-%s' % (narrowpeak_fn)
-
steps = ['slopBed -i %s -g %s -b 0' % (int_narrowpeak_fn, chrom_sizes),
'bedClip stdin %s %s' % (chrom_sizes, clipped_narrowpeak_fn)]
diff --git a/workflow/scripts/check_design.py b/workflow/scripts/check_design.py
index 35af1f469826933578ba0df5694f45644836b8b2..1b7d0f9868d6792e79a5c6ae8eeaa8b8b117f452 100644
--- a/workflow/scripts/check_design.py
+++ b/workflow/scripts/check_design.py
@@ -1,5 +1,11 @@
#!/usr/bin/env python3
+#
+# * --------------------------------------------------------------------------
+# * Licensed under MIT (https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/LICENSE.md)
+# * --------------------------------------------------------------------------
+#
+
'''Check if design file is correctly formatted and matches files list.'''
import argparse
diff --git a/workflow/scripts/convert_reads.py b/workflow/scripts/convert_reads.py
index bba7c6f16e49b7ea608d50546ec0b5af21ed4448..63afb364bb02ca859c7ca8c02fe9f9c59a0f6632 100644
--- a/workflow/scripts/convert_reads.py
+++ b/workflow/scripts/convert_reads.py
@@ -1,5 +1,11 @@
#!/usr/bin/env python3
+#
+# * --------------------------------------------------------------------------
+# * Licensed under MIT (https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/LICENSE.md)
+# * --------------------------------------------------------------------------
+#
+
'''Convert tagAlign files for further processing.'''
import os
@@ -54,6 +60,15 @@ def check_tools():
bedtools_path = shutil.which("bedtools")
if bedtools_path:
logger.info('Found bedtools: %s', bedtools_path)
+
+ # Get Version
+ bedtools_version_command = "bedtools --version"
+ bedtools_version = subprocess.check_output(bedtools_version_command, shell=True)
+
+ # Write to file
+ bedtools_file = open("version_bedtools.txt", "wb")
+ bedtools_file.write(bedtools_version)
+ bedtools_file.close()
else:
logger.error('Missing bedtools')
raise Exception('Missing bedtools')
@@ -61,6 +76,15 @@ def check_tools():
samtools_path = shutil.which("samtools")
if samtools_path:
logger.info('Found samtools: %s', samtools_path)
+
+ # Get Version
+ samtools_version_command = "samtools --version"
+ samtools_version = subprocess.check_output(samtools_version_command, shell=True)
+
+ # Write to file
+ samtools_file = open("version_samtools.txt", "wb")
+ samtools_file.write(samtools_version)
+ samtools_file.close()
else:
logger.error('Missing samtools')
raise Exception('Missing samtools')
diff --git a/workflow/scripts/diff_peaks.R b/workflow/scripts/diff_peaks.R
index aae1d24751f24396fff6fe712c28be27e36b0db3..3c5401e20b1c58235767f7c5a1e511c305c0a8db 100644
--- a/workflow/scripts/diff_peaks.R
+++ b/workflow/scripts/diff_peaks.R
@@ -1,5 +1,11 @@
#!/bin/Rscript
+#*
+#* --------------------------------------------------------------------------
+#* Licensed under MIT (https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/LICENSE.md)
+#* --------------------------------------------------------------------------
+#*
+
# Load libraries
library("DiffBind")
@@ -11,6 +17,11 @@ if (length(args) != 1) {
stop("Usage: diff_peaks.R annotate_design.tsv ", call.=FALSE)
}
+# Output version of DiffBind
+diffibind_version = packageVersion('DiffBind')
+write.table(paste("Version", diffibind_version), file = "version_DiffBind.txt", sep = "\t",
+ row.names = FALSE, col.names = FALSE)
+
# Build DBA object from design file
data <- dba(sampleSheet=args[1])
data <- dba.count(data)
diff --git a/workflow/scripts/experiment_design.py b/workflow/scripts/experiment_design.py
index f527b46d32f3dda20be3c84a4a8db822e9cb7310..728aefabdab32bbc80a3aa151ee9a8d8dabc5943 100644
--- a/workflow/scripts/experiment_design.py
+++ b/workflow/scripts/experiment_design.py
@@ -1,5 +1,11 @@
#!/usr/bin/env python3
+#
+# * --------------------------------------------------------------------------
+# * Licensed under MIT (https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/LICENSE.md)
+# * --------------------------------------------------------------------------
+#
+
'''Generate experiment design files for downstream processing.'''
import argparse
diff --git a/workflow/scripts/experiment_qc.py b/workflow/scripts/experiment_qc.py
index d884eaf5c572607dbbf2a4c51e59dcea27fe7087..608afaf038e2504b50695a61fc51a1878fb46b60 100644
--- a/workflow/scripts/experiment_qc.py
+++ b/workflow/scripts/experiment_qc.py
@@ -1,5 +1,11 @@
#!/usr/bin/env python3
+#
+# * --------------------------------------------------------------------------
+# * Licensed under MIT (https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/LICENSE.md)
+# * --------------------------------------------------------------------------
+#
+
'''Experiment correlation and enrichment of reads over genome-wide bins.'''
@@ -52,6 +58,15 @@ def check_tools():
deeptools_path = shutil.which("deeptools")
if deeptools_path:
logger.info('Found deeptools: %s', deeptools_path)
+
+ # Get Version
+ deeptools_version_command = "deeptools --version"
+ deeptools_version = subprocess.check_output(deeptools_version_command, shell=True, stderr=subprocess.STDOUT)
+
+ # Write to file
+ deeptools_file = open("version_deeptools.txt", "wb")
+ deeptools_file.write(deeptools_version)
+ deeptools_file.close()
else:
logger.error('Missing deeptools')
raise Exception('Missing deeptools')
@@ -77,10 +92,10 @@ def generate_read_summary(design, extension):
return mbs_filename
-def check_correlation(mbs):
+def check_spearman_correlation(mbs):
'''Check Spearman pairwise correlation of samples based on read counts.'''
- spearman_filename = 'heatmap_SpearmanCorr.png'
+ spearman_filename = 'heatmap_SpearmanCorr.pdf'
spearman_params = "--corMethod spearman --skipZero" + \
" --plotTitle \"Spearman Correlation of Read Counts\"" + \
" --whatToPlot heatmap --colorMap RdYlBu --plotNumbers"
@@ -96,12 +111,31 @@ def check_correlation(mbs):
out, err = spearman_correlation.communicate()
+def check_pearson_correlation(mbs):
+ '''Check Pearson pairwise correlation of samples based on read counts.'''
+
+ pearson_filename = 'heatmap_PearsonCorr.pdf'
+ pearson_params = "--corMethod pearson --skipZero" + \
+ " --plotTitle \"Pearson Correlation of Read Counts\"" + \
+ " --whatToPlot heatmap --colorMap RdYlBu --plotNumbers"
+
+ pearson_command = (
+ "plotCorrelation -in %s %s -o %s"
+ % (mbs, pearson_params, pearson_filename)
+ )
+
+ logger.info("Running deeptools with %s", pearson_command)
+
+ pearson_correlation = subprocess.Popen(pearson_command, shell=True)
+ out, err = pearson_correlation.communicate()
+
+
def check_coverage(design, extension):
'''Asses the sequencing depth of samples.'''
bam_files = ' '.join(design['bam_reads'])
labels = ' '.join(design['sample_id'])
- coverage_filename = 'coverage.png'
+ coverage_filename = 'coverage.pdf'
coverage_params = "-n 1000000 --plotTitle \"Sample Coverage\"" + \
" --ignoreDuplicates --minMappingQuality 10"
@@ -137,7 +171,7 @@ def update_controls(design):
def check_enrichment(sample_id, control_id, sample_reads, control_reads, extension):
'''Asses the enrichment per sample.'''
- fingerprint_filename = sample_id + '_fingerprint.png'
+ fingerprint_filename = sample_id + '_fingerprint.pdf'
fingerprint_command = (
"plotFingerprint -b %s %s --extendReads %d --labels %s %s --plotFile %s"
@@ -167,7 +201,8 @@ def main():
# Run correlation
mbs_filename = generate_read_summary(design_df, extension)
- check_correlation(mbs_filename)
+ check_spearman_correlation(mbs_filename)
+ check_pearson_correlation(mbs_filename)
# Run coverage
check_coverage(design_df, extension)
diff --git a/workflow/scripts/generate_references.py b/workflow/scripts/generate_references.py
new file mode 100644
index 0000000000000000000000000000000000000000..e20c7e81f2e01a8e694072e1cf2787709d2768e2
--- /dev/null
+++ b/workflow/scripts/generate_references.py
@@ -0,0 +1,73 @@
+#!/usr/bin/env python
+
+#
+# * --------------------------------------------------------------------------
+# * Licensed under MIT (https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/LICENSE.md)
+# * --------------------------------------------------------------------------
+#
+
+'''Make header for HTML of references.'''
+
+import argparse
+import subprocess
+import shlex
+import logging
+
+EPILOG = '''
+For more details:
+ %(prog)s --help
+'''
+
+# SETTINGS
+
+logger = logging.getLogger(__name__)
+logger.addHandler(logging.NullHandler())
+logger.propagate = False
+logger.setLevel(logging.INFO)
+
+
+def get_args():
+ '''Define arguments.'''
+
+ parser = argparse.ArgumentParser(
+ description=__doc__, epilog=EPILOG,
+ formatter_class=argparse.RawDescriptionHelpFormatter)
+
+ parser.add_argument('-r', '--reference',
+ help="The reference file (markdown format).",
+ required=True)
+
+ parser.add_argument('-o', '--output',
+ help="The out file name.",
+ default='references')
+
+ args = parser.parse_args()
+ return args
+
+
+def main():
+ args = get_args()
+ reference = args.reference
+ output = args.output
+
+ out_filename = output + '_mqc.yaml'
+
+ # Header for HTML
+ print('''
+ id: 'Software References'
+ section_name: 'Software References'
+ description: 'This section describes references for the tools used.'
+ plot_type: 'html'
+ data: |
+ '''
+ , file = open(out_filename, "w")
+ )
+
+ # Turn Markdown into HTML
+ references_html = 'bash -c "pandoc -p {} | sed \'s/^/ /\' >> {}"'
+ references_html = references_html.format(reference, out_filename)
+ subprocess.check_call(shlex.split(references_html))
+
+
+if __name__ == '__main__':
+ main()
diff --git a/workflow/scripts/generate_versions.py b/workflow/scripts/generate_versions.py
new file mode 100644
index 0000000000000000000000000000000000000000..4f0d8b143b36efbe207a303b2eb164d0de1a0e8a
--- /dev/null
+++ b/workflow/scripts/generate_versions.py
@@ -0,0 +1,148 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+#
+# * --------------------------------------------------------------------------
+# * Licensed under MIT (https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/LICENSE.md)
+# * --------------------------------------------------------------------------
+#
+
+'''Make YAML of software versions.'''
+
+from __future__ import print_function
+from collections import OrderedDict
+import re
+import os
+import logging
+import glob
+import argparse
+import numpy as np
+
+EPILOG = '''
+For more details:
+ %(prog)s --help
+'''
+
+# SETTINGS
+
+logger = logging.getLogger(__name__)
+logger.addHandler(logging.NullHandler())
+logger.propagate = False
+logger.setLevel(logging.INFO)
+
+SOFTWARE_REGEX = {
+ 'Nextflow': ['version_nextflow.txt', r"(\S+)"],
+ 'Trim Galore!': ['trimReads_vf/version_trimgalore.txt', r"version (\S+)"],
+ 'Cutadapt': ['trimReads_vf/version_cutadapt.txt', r"Version (\S+)"],
+ 'BWA': ['alignReads_vf/version_bwa.txt', r"Version: (\S+)"],
+ 'Samtools': ['alignReads_vf/version_samtools.txt', r"samtools (\S+)"],
+ 'Sambamba': ['filterReads_vf/version_sambamba.txt', r"sambamba (\S+)"],
+ 'BEDTools': ['convertReads_vf/version_bedtools.txt', r"bedtools v(\S+)"],
+ 'R': ['crossReads_vf/version_r.txt', r"R version (\S+)"],
+ 'SPP': ['crossReads_vf/version_spp.txt', r"\[1\] ‘(1.14)’"],
+ 'MACS2': ['callPeaksMACS_vf/version_macs.txt', r"macs2 (\S+)"],
+ 'bedGraphToBigWig': ['callPeaksMACS_vf/version_bedGraphToBigWig.txt', r"bedGraphToBigWig v (\S+)"],
+ 'ChIPseeker': ['peakAnnotation_vf/version_ChIPseeker.txt', r"Version (\S+)\""],
+ 'MEME-ChIP': ['motifSearch_vf/version_memechip.txt', r"Version (\S+)"],
+ 'DiffBind': ['diffPeaks_vf/version_DiffBind.txt', r"Version (\S+)\""],
+ 'deepTools': ['experimentQC_vf/version_deeptools.txt', r"deeptools (\S+)"],
+ 'Python': ['version_python.txt', r"python, version (\S+)"],
+ 'MultiQC': ['version_multiqc.txt', r"multiqc, version (\S+)"],
+}
+
+
+def get_args():
+ '''Define arguments.'''
+
+ parser = argparse.ArgumentParser(
+ description=__doc__, epilog=EPILOG,
+ formatter_class=argparse.RawDescriptionHelpFormatter)
+
+ parser.add_argument('-o', '--output',
+ help="The out file name.",
+ required=True)
+
+ parser.add_argument('-t', '--test',
+ help='Used for testing purposes',
+ default=False,
+ action='store_true')
+
+ args = parser.parse_args()
+ return args
+
+
+def check_files(files, test):
+ '''Check if version files are found.'''
+
+ logger.info("Running file check.")
+
+ software_files = np.array(list(SOFTWARE_REGEX.values()))[:,0]
+
+ extra_files = set(files) - set(software_files)
+
+ if len(extra_files) > 0 and test:
+ logger.error('Missing regex: %s', list(extra_files))
+ raise Exception("Missing regex: %s" % list(extra_files))
+
+
+def main():
+ args = get_args()
+ output = args.output
+ test = args.test
+
+ out_filename = output + '_mqc.yaml'
+
+ results = OrderedDict()
+ results['Nextflow'] = '<span style="color:#999999;\">Not Run</span>'
+ results['Trim Galore!'] = '<span style="color:#999999;\">Not Run</span>'
+ results['Cutadapt'] = '<span style="color:#999999;\">Not Run</span>'
+ results['BWA'] = '<span style="color:#999999;\">Not Run</span>'
+ results['Samtools'] = '<span style="color:#999999;\">Not Run</span>'
+ results['Sambamba'] = '<span style="color:#999999;\">Not Run</span>'
+ results['BEDTools'] = '<span style="color:#999999;\">Not Run</span>'
+ results['R'] = '<span style="color:#999999;\">Not Run</span>'
+ results['SPP'] = '<span style="color:#999999;\">Not Run</span>'
+ results['MACS2'] = '<span style="color:#999999;\">Not Run</span>'
+ results['bedGraphToBigWig'] = '<span style="color:#999999;\">Not Run</span>'
+ results['ChIPseeker'] = '<span style="color:#999999;\">Not Run</span>'
+ results['MEME-ChIP'] = '<span style="color:#999999;\">Not Run</span>'
+ results['DiffBind'] = '<span style="color:#999999;\">Not Run</span>'
+ results['deepTools'] = '<span style="color:#999999;\">Not Run</span>'
+ results['MultiQC'] = '<span style="color:#999999;\">Not Run</span>'
+ results['Python'] = '<span style="color:#999999;\">Not Run</span>'
+
+ # list all files
+ files = glob.glob('**/*.txt', recursive=True)
+
+ # Check for version files:
+ check_files(files, test)
+
+ # Search each file using its regex
+ for k, v in SOFTWARE_REGEX.items():
+ if os.path.isfile(v[0]):
+ with open(v[0]) as x:
+ versions = x.read()
+ match = re.search(v[1], versions)
+ if match:
+ results[k] = "v{}".format(match.group(1))
+
+ # Dump to YAML
+ print(
+ '''
+ id: 'Software Versions'
+ section_name: 'Software Versions'
+ section_href: 'https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/'
+ plot_type: 'html'
+ description: 'are collected at run time from the software output.'
+ data: |
+ <dl class="dl-horizontal">
+ '''
+ , file = open(out_filename, "w"))
+
+ for k, v in results.items():
+ print(" <dt>{}</dt><dd>{}</dd>".format(k, v), file = open(out_filename, "a"))
+ print(" </dl>", file = open(out_filename, "a"))
+
+
+if __name__ == '__main__':
+ main()
diff --git a/workflow/scripts/map_qc.py b/workflow/scripts/map_qc.py
index 920e0090a3fe69d50dd014e8541b5b9916cabecb..70eb84c908991669fea49af14fa3ce137b28dfe0 100644
--- a/workflow/scripts/map_qc.py
+++ b/workflow/scripts/map_qc.py
@@ -1,5 +1,11 @@
#!/usr/bin/env python3
+#
+# * --------------------------------------------------------------------------
+# * Licensed under MIT (https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/LICENSE.md)
+# * --------------------------------------------------------------------------
+#
+
'''Remove duplicates and filter unmapped reads.'''
import os
@@ -62,6 +68,15 @@ def check_tools():
samtools_path = shutil.which("samtools")
if samtools_path:
logger.info('Found samtools: %s', samtools_path)
+
+ # Get Version
+ samtools_version_command = "samtools --version"
+ samtools_version = subprocess.check_output(samtools_version_command, shell=True)
+
+ # Write to file
+ samtools_file = open("version_samtools.txt", "wb")
+ samtools_file.write(samtools_version)
+ samtools_file.close()
else:
logger.error('Missing samtools')
raise Exception('Missing samtools')
@@ -69,6 +84,18 @@ def check_tools():
sambamba_path = shutil.which("sambamba")
if sambamba_path:
logger.info('Found sambamba: %s', sambamba_path)
+
+ # Get Version
+ sambamba_version_command = "sambamba"
+ try:
+ subprocess.check_output(sambamba_version_command, shell=True, stderr=subprocess.STDOUT)
+ except subprocess.CalledProcessError as e:
+ sambamba_version = e.output
+
+ # Write to file
+ sambamba_file = open("version_sambamba.txt", "wb")
+ sambamba_file.write(sambamba_version)
+ sambamba_file.close()
else:
logger.error('Missing sambamba')
raise Exception('Missing sambamba')
@@ -76,6 +103,15 @@ def check_tools():
bedtools_path = shutil.which("bedtools")
if bedtools_path:
logger.info('Found bedtools: %s', bedtools_path)
+
+ # Get Version
+ bedtools_version_command = "bedtools --version"
+ bedtools_version = subprocess.check_output(bedtools_version_command, shell=True)
+
+ # Write to file
+ bedtools_file = open("version_bedtools.txt", "wb")
+ bedtools_file.write(bedtools_version)
+ bedtools_file.close()
else:
logger.error('Missing bedtools')
raise Exception('Missing bedtools')
@@ -167,7 +203,6 @@ def dedup_mapped(bam, bam_basename, paired):
shlex.split(sambamba_markdup_command),
stderr=temp_file)
-
# Remove duplicates
final_bam_prefix = bam_basename + ".dedup"
final_bam_filename = final_bam_prefix + ".bam"
@@ -215,6 +250,7 @@ def compute_complexity(bam, paired, bam_basename):
# Obtain unique count statistics
# PBC File output
+ # Sample Name[tab]
# TotalReadPairs [tab]
# DistinctReadPairs [tab]
# OneReadPair [tab]
@@ -249,9 +285,13 @@ def compute_complexity(bam, paired, bam_basename):
if err:
logger.error("PBC file error: %s", err)
- # Add headers
+ # Add Sample Name and headers
pbc_file = pd.read_csv(tmp_pbc_file_qc_filename, sep='\t', header=None,
names=pbc_headers)
+ pbc_file['Sample'] = bam_basename
+ pbc_headers_new = list(pbc_file)
+ pbc_headers_new.insert(0, pbc_headers_new.pop(pbc_headers_new.index('Sample')))
+ pbc_file = pbc_file[pbc_headers_new]
pbc_file.to_csv(pbc_file_qc_filename, header=True, sep='\t', index=False)
os.remove(bam)
os.remove(bam + '.bai')
diff --git a/workflow/scripts/map_reads.py b/workflow/scripts/map_reads.py
index a1f8161cb4ee71885719d954e9ddf25428499d11..9c09b0a9254ad0749cd36356212dd494b3162bf7 100644
--- a/workflow/scripts/map_reads.py
+++ b/workflow/scripts/map_reads.py
@@ -1,5 +1,11 @@
#!/usr/bin/env python3
+#
+# * --------------------------------------------------------------------------
+# * Licensed under MIT (https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/LICENSE.md)
+# * --------------------------------------------------------------------------
+#
+
'''Align reads to reference genome.'''
import os
@@ -70,6 +76,18 @@ def check_tools():
bwa_path = shutil.which("bwa")
if bwa_path:
logger.info('Found bwa: %s', bwa_path)
+
+ # Get Version
+ bwa_version_command = "bwa"
+ try:
+ subprocess.check_output(bwa_version_command, shell=True, stderr=subprocess.STDOUT)
+ except subprocess.CalledProcessError as e:
+ bwa_version = e.output
+
+ # Write to file
+ bwa_file = open("version_bwa.txt", "wb")
+ bwa_file.write(bwa_version)
+ bwa_file.close()
else:
logger.error('Missing bwa')
raise Exception('Missing bwa')
@@ -77,6 +95,15 @@ def check_tools():
samtools_path = shutil.which("samtools")
if samtools_path:
logger.info('Found samtools: %s', samtools_path)
+
+ # Get Version
+ samtools_version_command = "samtools --version"
+ samtools_version = subprocess.check_output(samtools_version_command, shell=True)
+
+ # Write to file
+ samtools_file = open("version_samtools.txt", "wb")
+ samtools_file.write(samtools_version)
+ samtools_file.close()
else:
logger.error('Missing samtools')
raise Exception('Missing samtools')
@@ -191,6 +218,13 @@ def main():
shlex.split("samtools flagstat %s" % (bam_filename)),
stdout=temp_file)
+ #Genome/Bad fastq File Check
+ file_check = open(bam_mapstats_filename).readlines()
+ percent = file_check[4].split('(')[1]
+ percent = percent.split('%')[0]
+ if float(percent) < 10:
+ raise Exception ('Mapped Genes too low: Check for correct Genotype')
+
# Remove sai files
for sai_file in sai:
os.remove(sai_file)
diff --git a/workflow/scripts/motif_search.py b/workflow/scripts/motif_search.py
index 02e316f67a7bdaf26541ca864ceebfcac50e120a..bc84e68687347162f2c65f065cc547ec617ece26 100644
--- a/workflow/scripts/motif_search.py
+++ b/workflow/scripts/motif_search.py
@@ -1,11 +1,18 @@
#!/usr/bin/env python3
+#
+# * --------------------------------------------------------------------------
+# * Licensed under MIT (https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/LICENSE.md)
+# * --------------------------------------------------------------------------
+#
+
'''Call Motifs on called peaks.'''
import os
import argparse
import logging
import shutil
+import subprocess
from multiprocessing import Pool
import pandas as pd
import utils
@@ -55,6 +62,7 @@ def get_args():
# Functions
+
def check_tools():
'''Checks for required componenets on user system'''
@@ -63,6 +71,15 @@ def check_tools():
meme_path = shutil.which("meme")
if meme_path:
logger.info('Found meme: %s', meme_path)
+
+ # Get Version
+ memechip_version_command = "meme-chip --version"
+ memechip_version = subprocess.check_output(memechip_version_command, shell=True)
+
+ # Write to file
+ meme_file = open("version_memechip.txt", "wb")
+ meme_file.write(b"Version %s" % (memechip_version))
+ meme_file.close()
else:
logger.error('Missing meme')
raise Exception('Missing meme')
@@ -70,6 +87,15 @@ def check_tools():
bedtools_path = shutil.which("bedtools")
if bedtools_path:
logger.info('Found bedtools: %s', bedtools_path)
+
+ # Get Version
+ bedtools_version_command = "bedtools --version"
+ bedtools_version = subprocess.check_output(bedtools_version_command, shell=True)
+
+ # Write to file
+ bedtools_file = open("version_bedtools.txt", "wb")
+ bedtools_file.write(bedtools_version)
+ bedtools_file.close()
else:
logger.error('Missing bedtools')
raise Exception('Missing bedtools')
@@ -95,7 +121,7 @@ def motif_search(filename, genome, experiment, peak):
else:
peak_no = peak
- sorted_fn = '%s.%s.narrowPeak' % (file_basename, peak)
+ sorted_fn = '%s.%s.narrowPeak' % (file_basename, peak_no)
out, err = utils.run_pipe([
'sort -k %dgr,%dgr %s' % (5, 5, filename),
@@ -108,8 +134,7 @@ def motif_search(filename, genome, experiment, peak):
if err:
logger.error("bedtools error: %s", err)
-
- #Call memechip
+ # Call memechip
out, err = utils.run_pipe([
'meme-chip -oc %s -meme-minw 5 -meme-maxw 15 -meme-nmotifs 10 %s -norand' % (out_motif, out_fa)])
if err:
diff --git a/workflow/scripts/overlap_peaks.py b/workflow/scripts/overlap_peaks.py
index 61f0c2e6d11553e54bd9cf5a1cc5c629fa9acb5c..3f1e79d70598f0861c6029939c620ddffdf8cbea 100644
--- a/workflow/scripts/overlap_peaks.py
+++ b/workflow/scripts/overlap_peaks.py
@@ -1,11 +1,18 @@
#!/usr/bin/env python3
+#
+# * --------------------------------------------------------------------------
+# * Licensed under MIT (https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/LICENSE.md)
+# * --------------------------------------------------------------------------
+#
+
'''Generate naive overlap peak files and design file for downstream processing.'''
import os
import argparse
import logging
import shutil
+import subprocess
import pandas as pd
import utils
@@ -49,6 +56,15 @@ def check_tools():
bedtools_path = shutil.which("bedtools")
if bedtools_path:
logger.info('Found bedtools: %s', bedtools_path)
+
+ # Get Version
+ bedtools_version_command = "bedtools --version"
+ bedtools_version = subprocess.check_output(bedtools_version_command, shell=True)
+
+ # Write to file
+ bedtools_file = open("version_bedtools.txt", "wb")
+ bedtools_file.write(bedtools_version)
+ bedtools_file.close()
else:
logger.error('Missing bedtools')
raise Exception('Missing bedtools')
@@ -178,6 +194,9 @@ def main():
handler = logging.FileHandler('consensus_peaks.log')
logger.addHandler(handler)
+ # Check if tools are present
+ check_tools()
+
# Read files as dataframes
design_peaks_df = pd.read_csv(design, sep='\t')
design_files_df = pd.read_csv(files, sep='\t')
diff --git a/workflow/scripts/plotProfile.sh b/workflow/scripts/plotProfile.sh
new file mode 100644
index 0000000000000000000000000000000000000000..0f50501c410d2ab0a7549190ba830aa0a90d7602
--- /dev/null
+++ b/workflow/scripts/plotProfile.sh
@@ -0,0 +1,16 @@
+#!/bin/bash
+#plotProfile.sh
+
+bws=$(ls *.bw)
+gtf=$(ls *.gtf *.bed)
+
+computeMatrix reference-point \
+ --referencePoint TSS \
+ -S $bws \
+ -R $gtf \
+ --skipZeros \
+ -o computeMatrix.gz
+ -p max/2
+
+plotProfile -m computeMatrix.gz \
+ -out plotProfile.png \
diff --git a/workflow/scripts/pool_and_psuedoreplicate.py b/workflow/scripts/pool_and_psuedoreplicate.py
index 7d27b971f038c33e81268728bd9aaec3450c0b33..6c37eed940ca89f71ebd3c49a2dae042f94d5b05 100644
--- a/workflow/scripts/pool_and_psuedoreplicate.py
+++ b/workflow/scripts/pool_and_psuedoreplicate.py
@@ -1,10 +1,19 @@
#!/usr/bin/env python3
+#
+# * --------------------------------------------------------------------------
+# * Licensed under MIT (https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/LICENSE.md)
+# * --------------------------------------------------------------------------
+#
+
'''Generate pooled and pseudoreplicate from data.'''
import argparse
import logging
import os
+import subprocess
+import shutil
+import shlex
import pandas as pd
import numpy as np
import utils
@@ -140,45 +149,30 @@ def self_psuedoreplication(tag_file, prefix, paired):
splits_prefix = 'temp_split'
- out, err = utils.run_pipe([
- 'gzip -dc %s' % (tag_file),
- 'shuf --random-source=%s' % (tag_file),
- 'split -d -l %d - %s' % (lines_per_rep, splits_prefix)])
+ psuedo_command = 'bash -c "zcat {} | shuf --random-source=<(openssl enc -aes-256-ctr -pass pass:$(zcat -f {} | wc -c) -nosalt </dev/zero 2>/dev/null) | '
+ psuedo_command += 'split -d -l {} - {}."'
+ psuedo_command = psuedo_command.format(
+ tag_file,
+ tag_file,
+ int(lines_per_rep),
+ splits_prefix)
+ logger.info("Running psuedo with %s", psuedo_command)
+ subprocess.check_call(shlex.split(psuedo_command))
# Convert read pairs to reads into standard tagAlign file
for i, index in enumerate([0, 1]):
- string_index = '0' + str(index)
+ string_index = '.0' + str(index)
steps = ['cat %s' % (splits_prefix + string_index)]
if paired:
steps.extend([r"""awk 'BEGIN{OFS="\t"}{printf "%s\t%s\t%s\tN\t1000\t%s\n%s\t%s\t%s\tN\t1000\t%s\n",$1,$2,$3,$9,$4,$5,$6,$10}'"""])
steps.extend(['gzip -cn'])
out, err = utils.run_pipe(steps, outfile=pseudoreplicate_dict[i])
- os.remove(splits_prefix + string_index)
return pseudoreplicate_dict
-def main():
- args = get_args()
- paired = args.paired
- design = args.design
- cutoff_ratio = args.cutoff
-
- # Create a file handler
- handler = logging.FileHandler('experiment_generation.log')
- logger.addHandler(handler)
-
- # Read files as dataframes
- design_df = pd.read_csv(design, sep='\t')
-
- # Get current directory to build paths
- cwd = os.getcwd()
-
- # Check Number of replicates and replicates
- no_reps = check_replicates(design_df)
- no_unique_controls = check_controls(design_df)
-
+def generate_design(paired, cutoff_ratio, design_df, cwd, no_reps, no_unique_controls):
if no_reps == 1:
logger.info("No other replicate specified "
"so processing as an unreplicated experiment.")
@@ -210,134 +204,119 @@ def main():
pool_control_tmp = bedpe_to_tagalign(pool_control, "pool_control")
pool_control = pool_control_tmp
- # Psuedoreplicate and update design accordingly
- if not replicated:
-
- # Duplicate rows and update for pool and psuedoreplicates and update tagAlign with single end data
- experiment_id = design_df.at[0, 'experiment_id']
- replicate = design_df.at[0, 'replicate']
- design_new_df = design_df.loc[np.repeat(design_df.index, 4)].reset_index()
-
- # Update tagAlign with single end data
- if paired:
- design_new_df['tag_align'] = design_new_df['se_tag_align']
- design_new_df.drop(labels='se_tag_align', axis=1, inplace=True)
-
- design_new_df['replicate'] = design_new_df['replicate'].astype(str)
- design_new_df.at[1, 'sample_id'] = experiment_id + '_pr1'
- design_new_df.at[1, 'replicate'] = '1_pr'
- design_new_df.at[1, 'xcor'] = 'Calculate'
- design_new_df.at[2, 'sample_id'] = experiment_id + '_pr2'
- design_new_df.at[2, 'replicate'] = '2_pr'
- design_new_df.at[2, 'xcor'] = 'Calculate'
- design_new_df.at[3, 'sample_id'] = experiment_id + '_pooled'
- design_new_df.at[3, 'replicate'] = 'pooled'
- design_new_df.at[3, 'xcor'] = 'Calculate'
- design_new_df.at[3, 'tag_align'] = design_new_df.at[0, 'tag_align']
-
- # Make 2 self psuedoreplicates
- self_pseudoreplicates_dict = {}
- for rep, tag_file in zip(design_df['replicate'], design_df['tag_align']):
- replicate_prefix = experiment_id + '_' + str(rep)
- self_pseudoreplicates_dict = \
- self_psuedoreplication(tag_file, replicate_prefix, paired)
-
- # Update design to include new self pseudo replicates
- for rep, pseudorep_file in self_pseudoreplicates_dict.items():
- path_to_file = cwd + '/' + pseudorep_file
- replicate = rep + 1
- design_new_df.loc[replicate, 'tag_align'] = path_to_file
-
- # Drop index column
- design_new_df.drop(labels='index', axis=1, inplace=True)
+ # Duplicate rows and update for pool and psuedoreplicates and update tagAlign with single end data
+ experiment_id = design_df.at[0, 'experiment_id']
+ replicate_files = design_df.tag_align.unique()
+ pool_experiment = pool(replicate_files, experiment_id + "_pooled", paired)
+
+ # Make 2 self psuedoreplicates
+ pseudoreplicates_dict = {}
+ for rep, tag_file in zip(design_df['replicate'], design_df['tag_align']):
+ replicate_prefix = experiment_id + '_' + str(rep)
+ pr_dict = self_psuedoreplication(tag_file, replicate_prefix, paired)
+ pseudoreplicates_dict[rep] = pr_dict
+
+ # Update design to include new self pseudo replicates
+ pseudoreplicates_df = pd.DataFrame.from_dict(pseudoreplicates_dict)
+ pool_pseudoreplicates_dict = {}
+ for index, row in pseudoreplicates_df.iterrows():
+ replicate_id = index + 1
+ pr_filename = experiment_id + ".pr" + str(replicate_id) + '.tagAlign.gz'
+ pool_replicate = pool(row, pr_filename, False)
+ pool_pseudoreplicates_dict[replicate_id] = pool_replicate
+
+ design_new_df = design_df #.loc[np.repeat(design_df.index, 4)].reset_index()
+ # Update tagAlign with single end data
+ if paired:
+ design_new_df['tag_align'] = design_new_df['se_tag_align']
+ design_new_df.drop(labels='se_tag_align', axis=1, inplace=True)
+ # If paired change to single End
+ if paired:
+ pool_experiment_se = bedpe_to_tagalign(pool_experiment, experiment_id + "_pooled")
else:
- # Make pool of replicates
- replicate_files = design_df.tag_align.unique()
- experiment_id = design_df.at[0, 'experiment_id']
- pool_experiment = pool(replicate_files, experiment_id + "_pooled", paired)
+ pool_experiment_se = pool_experiment
- # If paired change to single End
- if paired:
- pool_experiment_se = bedpe_to_tagalign(pool_experiment, experiment_id + "_pooled")
- else:
- pool_experiment_se = pool_experiment
-
- # Make self psuedoreplicates equivalent to number of replicates
- pseudoreplicates_dict = {}
- for rep, tag_file in zip(design_df['replicate'], design_df['tag_align']):
- replicate_prefix = experiment_id + '_' + str(rep)
- pr_dict = self_psuedoreplication(tag_file, replicate_prefix, paired)
- pseudoreplicates_dict[rep] = pr_dict
-
- # Merge self psuedoreplication for each true replicate
- pseudoreplicates_df = pd.DataFrame.from_dict(pseudoreplicates_dict)
- pool_pseudoreplicates_dict = {}
- for index, row in pseudoreplicates_df.iterrows():
- replicate_id = index + 1
- pr_filename = experiment_id + ".pr" + str(replicate_id) + '.tagAlign.gz'
- pool_replicate = pool(row, pr_filename, False)
- pool_pseudoreplicates_dict[replicate_id] = pool_replicate
-
- design_new_df = design_df
- # Update tagAlign with single end data
- if paired:
- design_new_df['tag_align'] = design_new_df['se_tag_align']
- design_new_df.drop(labels='se_tag_align', axis=1, inplace=True)
# Check controls against cutoff_ratio
# if so replace with pool_control
# unless single control was used
-
- if not single_control:
- path_to_pool_control = cwd + '/' + pool_control
- if control_df.values.max() > cutoff_ratio:
- logger.info("Number of reads in controls differ by " +
- " > factor of %f. Using pooled controls." % (cutoff_ratio))
- design_new_df['control_tag_align'] = path_to_pool_control
- else:
- for index, row in design_new_df.iterrows():
- exp_no_reads = utils.count_lines(row['tag_align'])
- con_no_reads = utils.count_lines(row['control_tag_align'])
- if con_no_reads < exp_no_reads:
- logger.info("Fewer reads in control than experiment." +
- "Using pooled controls for replicate %s."
- % row['replicate'])
+ if not single_control:
+ path_to_pool_control = cwd + '/' + pool_control
+ if control_df.values.max() > cutoff_ratio:
+ logger.info("Number of reads in controls differ by " +
+ " > factor of %f. Using pooled controls." % (cutoff_ratio))
+ design_new_df['control_tag_align'] = path_to_pool_control
+ else:
+ for index, row in design_new_df.iterrows():
+ exp_no_reads = utils.count_lines(row['tag_align'])
+ con_no_reads = utils.count_lines(row['control_tag_align'])
+ if con_no_reads < exp_no_reads:
+ logger.info("Fewer reads in control than experiment." +
+ "Using pooled controls for replicate %s."
+ % row['replicate'])
+ design_new_df.loc[index, 'control_tag_align'] = \
+ path_to_pool_control
+ else:
+ if paired:
+ control = row['control_tag_align']
+ control_basename = os.path.basename(
+ utils.strip_extensions(control, STRIP_EXTENSIONS))
+ control_tmp = bedpe_to_tagalign(control, control_basename)
+ path_to_control = cwd + '/' + control_tmp
design_new_df.loc[index, 'control_tag_align'] = \
- path_to_pool_control
- else:
- if paired:
- control = row['control_tag_align']
- control_basename = os.path.basename(
- utils.strip_extensions(control, STRIP_EXTENSIONS))
- control_tmp = bedpe_to_tagalign(control, control_basename)
- path_to_control = cwd + '/' + control_tmp
- design_new_df.loc[index, 'control_tag_align'] = \
- path_to_control
+ path_to_control
- else:
- path_to_pool_control = pool_control
-
- # Add in pseudo replicates
- tmp_metadata = design_new_df.loc[0].copy()
- tmp_metadata['control_tag_align'] = path_to_pool_control
- for rep, pseudorep_file in pool_pseudoreplicates_dict.items():
- tmp_metadata['sample_id'] = experiment_id + '_pr' + str(rep)
- tmp_metadata['replicate'] = str(rep) + '_pr'
- tmp_metadata['xcor'] = 'Calculate'
- path_to_file = cwd + '/' + pseudorep_file
- tmp_metadata['tag_align'] = path_to_file
- design_new_df = design_new_df.append(tmp_metadata)
-
- # Add in pool experiment
- tmp_metadata['sample_id'] = experiment_id + '_pooled'
- tmp_metadata['replicate'] = 'pooled'
+ else:
+ path_to_pool_control = cwd + '/' + pool_control
+ design_new_df['control_tag_align'] = path_to_pool_control
+
+ # Add in pseudo replicates
+ tmp_metadata = design_new_df.loc[0].copy()
+ tmp_metadata['control_tag_align'] = path_to_pool_control
+ for rep, pseudorep_file in pool_pseudoreplicates_dict.items():
+ tmp_metadata['sample_id'] = experiment_id + '_pr' + str(rep)
+ tmp_metadata['replicate'] = str(rep) + '_pr'
tmp_metadata['xcor'] = 'Calculate'
- path_to_file = cwd + '/' + pool_experiment_se
+ path_to_file = cwd + '/' + pseudorep_file
tmp_metadata['tag_align'] = path_to_file
design_new_df = design_new_df.append(tmp_metadata)
+ # Add in pool experiment
+ tmp_metadata['sample_id'] = experiment_id + '_pooled'
+ tmp_metadata['replicate'] = 'pooled'
+ tmp_metadata['xcor'] = 'Calculate'
+ path_to_file = cwd + '/' + pool_experiment_se
+ tmp_metadata['tag_align'] = path_to_file
+ design_new_df = design_new_df.append(tmp_metadata)
+
+ return design_new_df
+
+
+def main():
+ args = get_args()
+ paired = args.paired
+ design = args.design
+ cutoff_ratio = args.cutoff
+
+ # Create a file handler
+ handler = logging.FileHandler('experiment_generation.log')
+ logger.addHandler(handler)
+
+ # Read files as dataframes
+ design_df = pd.read_csv(design, sep='\t')
+
+ # Get current directory to build paths
+ cwd = os.getcwd()
+
+ # Check Number of replicates and replicates
+ no_reps = check_replicates(design_df)
+ no_unique_controls = check_controls(design_df)
+
+ # Generate new design file
+ design_new_df = generate_design(paired, cutoff_ratio, design_df, cwd, no_reps, no_unique_controls)
# Write out new dataframe
+ experiment_id = design_df.at[0, 'experiment_id']
design_new_df.to_csv(experiment_id + '_ppr.tsv',
header=True, sep='\t', index=False)
diff --git a/workflow/scripts/trim_reads.py b/workflow/scripts/trim_reads.py
index e31ec938b6d7e0ef5570afc394c98a365041f863..717df8e77dffcf4598933fd42ae6a66d2fc277d1 100644
--- a/workflow/scripts/trim_reads.py
+++ b/workflow/scripts/trim_reads.py
@@ -1,5 +1,11 @@
#!/usr/bin/env python3
+#
+# * --------------------------------------------------------------------------
+# * Licensed under MIT (https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/LICENSE.md)
+# * --------------------------------------------------------------------------
+#
+
'''Trim low quality reads and remove sequences less than 35 base pairs.'''
import subprocess
@@ -54,6 +60,15 @@ def check_tools():
trimgalore_path = shutil.which("trim_galore")
if trimgalore_path:
logger.info('Found trimgalore: %s', trimgalore_path)
+
+ # Get Version
+ trim_version_command = "trim_galore --version"
+ trimgalore_version = subprocess.check_output(trim_version_command, shell=True)
+
+ # Write to file
+ trimgalore_file = open("version_trimgalore.txt", "wb")
+ trimgalore_file.write(trimgalore_version)
+ trimgalore_file.close()
else:
logger.error('Missing trimgalore')
raise Exception('Missing trimgalore')
@@ -61,6 +76,15 @@ def check_tools():
cutadapt_path = shutil.which("cutadapt")
if cutadapt_path:
logger.info('Found cutadapt: %s', cutadapt_path)
+
+ # Get Version
+ cutadapt_version_command = "cutadapt --version"
+ cutadapt_version = subprocess.check_output(cutadapt_version_command, shell=True)
+
+ # Write to file
+ cutadapt_file = open("version_cutadapt.txt", "wb")
+ cutadapt_file.write(b"Version %s" % (cutadapt_version))
+ cutadapt_file.close()
else:
logger.error('Missing cutadapt')
raise Exception('Missing cutadapt')
diff --git a/workflow/scripts/utils.py b/workflow/scripts/utils.py
index 6c2da13b056b386e454a309552f5a3623b1ef254..5e4478e3b1cef96e9b4c05033f75122f87aad06e 100644
--- a/workflow/scripts/utils.py
+++ b/workflow/scripts/utils.py
@@ -1,5 +1,11 @@
#!/usr/bin/env python3
+#
+# * --------------------------------------------------------------------------
+# * Licensed under MIT (https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/LICENSE.md)
+# * --------------------------------------------------------------------------
+#
+
'''General utilities.'''
diff --git a/workflow/scripts/xcor.py b/workflow/scripts/xcor.py
index ac2ff65f7f06707e892e1033106dc2219ecb96cd..66fa9aec4d7030288b3bd57d7e156247d82c97ec 100644
--- a/workflow/scripts/xcor.py
+++ b/workflow/scripts/xcor.py
@@ -1,10 +1,17 @@
#!/usr/bin/env python3
+#
+# * --------------------------------------------------------------------------
+# * Licensed under MIT (https://git.biohpc.swmed.edu/BICF/Astrocyte/chipseq_analysis/LICENSE.md)
+# * --------------------------------------------------------------------------
+#
+
'''Compute cross-correlation analysis.'''
import os
import argparse
import shutil
+import subprocess
import logging
from multiprocessing import cpu_count
import utils
@@ -59,10 +66,35 @@ def check_tools():
r_path = shutil.which("R")
if r_path:
logger.info('Found R: %s', r_path)
+
+ # Get Version
+ r_version_command = "R --version"
+ r_version = subprocess.check_output(r_version_command, shell=True)
+
+ # Write to file
+ r_file = open("version_r.txt", "wb")
+ r_file.write(r_version)
+ r_file.close()
else:
logger.error('Missing R')
raise Exception('Missing R')
+ phantompeak_path = shutil.which("run_spp.R")
+ if phantompeak_path:
+ logger.info('Found phantompeak: %s', phantompeak_path)
+
+ # Get Version
+ spp_version_command = "R -e \"packageVersion('spp')\""
+ spp_version = subprocess.check_output(spp_version_command, shell=True)
+
+ # Write to file
+ spp_file = open("version_spp.txt", "wb")
+ spp_file.write(spp_version)
+ spp_file.close()
+ else:
+ logger.error('Missing phantompeak')
+ raise Exception('Missing phantompeak')
+
def xcor(tag, paired):
'''Use spp to calculate cross-correlation stats.'''
@@ -70,17 +102,21 @@ def xcor(tag, paired):
tag_basename = os.path.basename(utils.strip_extensions(tag, STRIP_EXTENSIONS))
uncompressed_tag_filename = tag_basename
-
# Subsample tagAlign file
- number_reads = 15000000
+ number_reads = 20000000
subsampled_tag_filename = \
tag_basename + ".%d.tagAlign.gz" % (number_reads/1000000)
+ tag_extended = 'cat.tagAlign.gz'
+ out, err = utils.run_pipe([
+ "zcat %s %s %s" %
+ (tag, tag, tag)
+ ], outfile=tag_extended)
steps = [
'zcat %s' % (tag),
'grep -v "chrM"',
- 'shuf -n %d --random-source=%s' % (number_reads, tag)]
+ 'shuf -n %d --random-source=%s' % (number_reads, tag_extended)]
if paired:
steps.extend([r"""awk 'BEGIN{OFS="\t"}{$4="N";$5="1000";print $0}'"""])
@@ -116,7 +152,6 @@ def xcor(tag, paired):
return cc_scores_filename
-
def main():
args = get_args()
paired = args.paired
diff --git a/workflow/tests/test_annotate_peaks.py b/workflow/tests/test_annotate_peaks.py
index ef42a9adc332c505715ad70be3b343de901898a8..73656357dc176169246c8196f20c60555e87dfac 100644
--- a/workflow/tests/test_annotate_peaks.py
+++ b/workflow/tests/test_annotate_peaks.py
@@ -19,11 +19,12 @@ def test_pie_singleend():
def test_upsetplot_singleend():
assert os.path.exists(os.path.join(test_output_path, 'ENCSR238SGC.chipseeker_upsetplot.pdf'))
+
@pytest.mark.singleend
def test_annotation_singleend():
annotation_file = test_output_path + 'ENCSR238SGC.chipseeker_annotation.tsv'
assert os.path.exists(annotation_file)
- assert utils.count_lines(annotation_file) == 152840
+ assert utils.count_lines(annotation_file) >= 149284
@pytest.mark.pairedend
@@ -40,4 +41,4 @@ def test_upsetplot_pairedend():
def test_annotation_pairedend():
annotation_file = test_output_path + 'ENCSR729LGA.chipseeker_annotation.tsv'
assert os.path.exists(annotation_file)
- assert utils.count_lines(annotation_file) == 25391
+ assert utils.count_lines(annotation_file) >= 25367
diff --git a/workflow/tests/test_call_peaks_macs.py b/workflow/tests/test_call_peaks_macs.py
index 91c1d0e43460ed663d1c59bf6caa8ad800e8fd61..cd94e1783a25ffa953b3666f7cb3172c33f4f0bd 100644
--- a/workflow/tests/test_call_peaks_macs.py
+++ b/workflow/tests/test_call_peaks_macs.py
@@ -10,41 +10,41 @@ test_output_path = os.path.dirname(os.path.abspath(__file__)) + \
@pytest.mark.singleend
def test_fc_signal_singleend():
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB144FDT.fc_signal.bw'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCSR238SGC/1/', 'ENCLB144FDT.fc_signal.bw'))
@pytest.mark.singleend
def test_pvalue_signal_singleend():
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB144FDT.pvalue_signal.bw'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCSR238SGC/1/', 'ENCLB144FDT.pvalue_signal.bw'))
@pytest.mark.singleend
def test_peaks_xls_singleend():
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB144FDT_peaks.xls'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCSR238SGC/1/', 'ENCLB144FDT_peaks.xls'))
@pytest.mark.singleend
def test_peaks_bed_singleend():
- peak_file = test_output_path + 'ENCLB144FDT.narrowPeak'
- assert utils.count_lines(peak_file) == 227389
+ peak_file = test_output_path + 'ENCSR238SGC/1/' + 'ENCLB144FDT.narrowPeak'
+ assert utils.count_lines(peak_file) == 226738
@pytest.mark.pairedend
def test_fc_signal_pairedend():
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB568IYX.fc_signal.bw'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCSR729LGA/2/', 'ENCLB568IYX.fc_signal.bw'))
@pytest.mark.pairedend
def test_pvalue_signal_pairedend():
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB568IYX.pvalue_signal.bw'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCSR729LGA/2/', 'ENCLB568IYX.pvalue_signal.bw'))
@pytest.mark.pairedend
def test_peaks_xls_pairedend():
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB568IYX_peak.xls'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCSR729LGA/2/', 'ENCLB568IYX_peaks.xls'))
@pytest.mark.pairedend
def test_peaks_bed_pairedend():
- peak_file = test_output_path + 'ENCLB568IYX.narrowPeak'
- assert utils.count_lines(peak_file) == 113821
+ peak_file = test_output_path + 'ENCSR729LGA/2/' + 'ENCLB568IYX.narrowPeak'
+ assert utils.count_lines(peak_file) == 112631
diff --git a/workflow/tests/test_convert_reads.py b/workflow/tests/test_convert_reads.py
index 8dda74678d42c0049eb74caa205ad8b662772bb5..bf4528234807779c38123bf02f66af57ce44bc21 100644
--- a/workflow/tests/test_convert_reads.py
+++ b/workflow/tests/test_convert_reads.py
@@ -9,19 +9,19 @@ test_output_path = os.path.dirname(os.path.abspath(__file__)) + \
@pytest.mark.singleend
def test_tag_reads_singleend():
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB831RUI.tagAlign.gz'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCLB831RUI/ENCLB831RUI.tagAlign.gz'))
@pytest.mark.singleend
def test_bed_reads_singleend():
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB831RUI.bedse.gz'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCLB831RUI/ENCLB831RUI.bedse.gz'))
@pytest.mark.pairedend
def test_tag_reads_pairedend():
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB568IYX.tagAlign.gz'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCLB568IYX/ENCLB568IYX.tagAlign.gz'))
@pytest.mark.pairedend
def test_bed_reads_pairedend():
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB568IYX.bedpe.gz'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCLB568IYX/ENCLB568IYX.bedpe.gz'))
diff --git a/workflow/tests/test_diff_peaks.py b/workflow/tests/test_diff_peaks.py
index 84c5179f021abd676a818ac24e698b54bd7b648d..47940c0b2c103fda1445043915a70b3c4e1464d0 100644
--- a/workflow/tests/test_diff_peaks.py
+++ b/workflow/tests/test_diff_peaks.py
@@ -10,14 +10,16 @@ test_output_path = os.path.dirname(os.path.abspath(__file__)) + \
'/../output/diffPeaks/'
-@pytest.mark.singleend
+@pytest.mark.singleskip_true
def test_diff_peaks_singleend_single_rep():
assert os.path.isdir(test_output_path) == False
-@pytest.mark.pairedend
+
+@pytest.mark.pairedendskip_true
def test_diff_peaks_pairedend_single_rep():
assert os.path.isdir(test_output_path) == False
+
@pytest.mark.singlediff
def test_heatmap_singleend_multiple_rep():
assert os.path.exists(os.path.join(test_output_path, 'heatmap.pdf'))
@@ -42,7 +44,7 @@ def test_diffbind_singleend_multiple_rep():
assert os.path.exists(os.path.join(test_output_path, 'ENCSR238SGC_vs_ENCSR272GNQ_diffbind.bed'))
diffbind_file = test_output_path + 'ENCSR238SGC_vs_ENCSR272GNQ_diffbind.csv'
assert os.path.exists(diffbind_file)
- assert utils.count_lines(diffbind_file) == 201039
+ assert utils.count_lines(diffbind_file) >= 197471
@pytest.mark.paireddiff
@@ -69,4 +71,4 @@ def test_diffbind_pairedend_single_rep():
assert os.path.exists(os.path.join(test_output_path, 'ENCSR729LGA_vs_ENCSR757EMK_diffbind.bed'))
diffbind_file = test_output_path + 'ENCSR729LGA_vs_ENCSR757EMK_diffbind.csv'
assert os.path.exists(diffbind_file)
- assert utils.count_lines(diffbind_file) == 66201
+ assert utils.count_lines(diffbind_file) >= 65182
diff --git a/workflow/tests/test_experiment_qc.py b/workflow/tests/test_experiment_qc.py
index 98853b9d6c2eabc4e4b9887cfb51c07a6af08fa4..c0cd5d7d340b27b3217620ef4b12a1391841820b 100644
--- a/workflow/tests/test_experiment_qc.py
+++ b/workflow/tests/test_experiment_qc.py
@@ -33,32 +33,42 @@ def test_check_update_controls(design_bam):
@pytest.mark.singleend
def test_coverage_singleend():
assert os.path.exists(os.path.join(test_output_path, 'sample_mbs.npz'))
- assert os.path.exists(os.path.join(test_output_path, 'coverage.png'))
+ assert os.path.exists(os.path.join(test_output_path, 'coverage.pdf'))
@pytest.mark.singleend
def test_spearman_singleend():
- assert os.path.exists(os.path.join(test_output_path, 'heatmap_SpearmanCorr.png'))
+ assert os.path.exists(os.path.join(test_output_path, 'heatmap_SpearmanCorr.pdf'))
+
+
+@pytest.mark.singleend
+def test_pearson_singleend():
+ assert os.path.exists(os.path.join(test_output_path, 'heatmap_PearsonCorr.pdf'))
@pytest.mark.singleend
def test_fingerprint_singleend():
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB144FDT_fingerprint.png'))
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB831RUI_fingerprint.png'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCLB144FDT_fingerprint.pdf'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCLB831RUI_fingerprint.pdf'))
@pytest.mark.pairdend
def test_coverage_pairedend():
assert os.path.exists(os.path.join(test_output_path, 'sample_mbs.npz'))
- assert os.path.exists(os.path.join(test_output_path, 'coverage.png'))
+ assert os.path.exists(os.path.join(test_output_path, 'coverage.pdf'))
@pytest.mark.pairdend
def test_spearman_pairedend():
- assert os.path.exists(os.path.join(test_output_path, 'heatmap_SpearmanCorr.png'))
+ assert os.path.exists(os.path.join(test_output_path, 'heatmap_SpearmanCorr.pdf'))
+
+
+@pytest.mark.pairdend
+def test_pearson_pairedend():
+ assert os.path.exists(os.path.join(test_output_path, 'heatmap_PearsonCorr.pdf'))
@pytest.mark.pairdend
def test_fingerprint_pairedend():
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB568IYX_fingerprint.png'))
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB637LZP_fingerprint.png'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCLB568IYX_fingerprint.pdf'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCLB637LZP_fingerprint.pdf'))
diff --git a/workflow/tests/test_generate_software_references.py b/workflow/tests/test_generate_software_references.py
new file mode 100644
index 0000000000000000000000000000000000000000..4d0025dfc3a925d2ec0d7bcd2bae86529d75e574
--- /dev/null
+++ b/workflow/tests/test_generate_software_references.py
@@ -0,0 +1,23 @@
+#!/usr/bin/env python3
+
+import pytest
+import os
+import utils
+import yaml
+
+test_output_path = os.path.dirname(os.path.abspath(__file__)) + \
+ '/../output/multiqcReport/'
+
+
+@pytest.mark.singleend
+def test_software_references():
+ assert os.path.exists(os.path.join(test_output_path, 'software_references_mqc.yaml'))
+
+
+@pytest.mark.singleend
+def test_software_references_output():
+ software_references = os.path.join(test_output_path, 'software_references_mqc.yaml')
+ with open(software_references, 'r') as stream:
+ data_loaded = yaml.load(stream)
+
+ assert len(data_loaded['data'].split('<ul>')) == 19
diff --git a/workflow/tests/test_generate_software_versions.py b/workflow/tests/test_generate_software_versions.py
new file mode 100644
index 0000000000000000000000000000000000000000..04ecdf3f57edee3507abaec27652023c8c5f5716
--- /dev/null
+++ b/workflow/tests/test_generate_software_versions.py
@@ -0,0 +1,23 @@
+#!/usr/bin/env python3
+
+import pytest
+import os
+import utils
+import yaml
+
+test_output_path = os.path.dirname(os.path.abspath(__file__)) + \
+ '/../output/multiqcReport/'
+
+
+@pytest.mark.singleend
+def test_software_versions():
+ assert os.path.exists(os.path.join(test_output_path, 'software_versions_mqc.yaml'))
+
+
+@pytest.mark.singleend
+def test_software_versions_output():
+ software_versions = os.path.join(test_output_path, 'software_versions_mqc.yaml')
+ with open(software_versions, 'r') as stream:
+ data_loaded = yaml.load(stream)
+
+ assert len(data_loaded['data'].split('<dt>')) == 18
diff --git a/workflow/tests/test_map_qc.py b/workflow/tests/test_map_qc.py
index 7bff0601d679dd5426e03c08c1670acbecf0ad22..ab94969b9699c204b30df31680fae088cf021c52 100644
--- a/workflow/tests/test_map_qc.py
+++ b/workflow/tests/test_map_qc.py
@@ -10,14 +10,14 @@ test_output_path = os.path.dirname(os.path.abspath(__file__)) + \
@pytest.mark.singleend
def test_dedup_files_singleend():
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB831RUI.dedup.bam'))
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB831RUI.dedup.bam.bai'))
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB831RUI.dedup.qc'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCLB831RUI/ENCLB831RUI.dedup.bam'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCLB831RUI/ENCLB831RUI.dedup.bam.bai'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCLB831RUI/ENCLB831RUI.dedup.qc'))
@pytest.mark.singleend
def test_map_qc_singleend():
- filtered_reads_report = test_output_path + 'ENCLB831RUI.dedup.flagstat.qc'
+ filtered_reads_report = test_output_path + 'ENCLB831RUI/ENCLB831RUI.dedup.flagstat.qc'
samtools_report = open(filtered_reads_report).readlines()
assert '64962570 + 0 in total' in samtools_report[0]
assert '64962570 + 0 mapped (100.00%:N/A)' in samtools_report[4]
@@ -25,7 +25,7 @@ def test_map_qc_singleend():
@pytest.mark.singleend
def test_library_complexity_singleend():
- library_complexity = test_output_path + 'ENCLB831RUI.pbc.qc'
+ library_complexity = test_output_path + 'ENCLB831RUI/ENCLB831RUI.pbc.qc'
df_library_complexity = pd.read_csv(library_complexity, sep='\t')
assert df_library_complexity["NRF"].iloc[0] == 0.926192
assert df_library_complexity["PBC1"].iloc[0] == 0.926775
@@ -34,14 +34,14 @@ def test_library_complexity_singleend():
@pytest.mark.pairedend
def test_dedup_files_pairedend():
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB568IYX.dedup.bam'))
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB568IYX.dedup.bam.bai'))
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB568IYX.dedup.qc'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCLB568IYX/ENCLB568IYX.dedup.bam'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCLB568IYX/ENCLB568IYX.dedup.bam.bai'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCLB568IYX/ENCLB568IYX.dedup.qc'))
@pytest.mark.pairedend
def test_map_qc_pairedend():
- filtered_reads_report = test_output_path + 'ENCLB568IYX.dedup.flagstat.qc'
+ filtered_reads_report = test_output_path + 'ENCLB568IYX/ENCLB568IYX.dedup.flagstat.qc'
samtools_report = open(filtered_reads_report).readlines()
assert '47388510 + 0 in total' in samtools_report[0]
assert '47388510 + 0 mapped (100.00%:N/A)' in samtools_report[4]
@@ -49,8 +49,8 @@ def test_map_qc_pairedend():
@pytest.mark.pairedend
def test_library_complexity_pairedend():
- library_complexity = test_output_path + 'ENCLB568IYX.pbc.qc'
+ library_complexity = test_output_path + 'ENCLB568IYX/ENCLB568IYX.pbc.qc'
df_library_complexity = pd.read_csv(library_complexity, sep='\t')
assert df_library_complexity["NRF"].iloc[0] == 0.947064
assert round(df_library_complexity["PBC1"].iloc[0],6) == 0.946723
- assert df_library_complexity["PBC2"].iloc[0] == 18.643039
+ assert round(df_library_complexity["PBC2"].iloc[0],6) == 18.642645
diff --git a/workflow/tests/test_map_reads.py b/workflow/tests/test_map_reads.py
index 60a2e39009f64f2a554d200ba41f1f4d28fe92a8..9ff817f452c24420788d6cbe75ce14a6989da167 100644
--- a/workflow/tests/test_map_reads.py
+++ b/workflow/tests/test_map_reads.py
@@ -9,8 +9,8 @@ test_output_path = os.path.dirname(os.path.abspath(__file__)) + \
@pytest.mark.singleend
def test_map_reads_singleend():
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB831RUI.bam'))
- aligned_reads_report = test_output_path + 'ENCLB831RUI.flagstat.qc'
+ assert os.path.exists(os.path.join(test_output_path, 'ENCLB831RUI/ENCLB831RUI.bam'))
+ aligned_reads_report = test_output_path + 'ENCLB831RUI/ENCLB831RUI.flagstat.qc'
samtools_report = open(aligned_reads_report).readlines()
assert '80795025 + 0 in total' in samtools_report[0]
assert '80050072 + 0 mapped (99.08% : N/A)' in samtools_report[4]
@@ -18,8 +18,8 @@ def test_map_reads_singleend():
@pytest.mark.pairedend
def test_map_reads_pairedend():
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB678IDC.bam'))
- aligned_reads_report = test_output_path + 'ENCLB678IDC.flagstat.qc'
+ assert os.path.exists(os.path.join(test_output_path, 'ENCLB678IDC/ENCLB678IDC.bam'))
+ aligned_reads_report = test_output_path + 'ENCLB678IDC/ENCLB678IDC.flagstat.qc'
samtools_report = open(aligned_reads_report).readlines()
assert '72660890 + 0 in total' in samtools_report[0]
assert '72053925 + 0 mapped (99.16% : N/A)' in samtools_report[4]
diff --git a/workflow/tests/test_motif_search.py b/workflow/tests/test_motif_search.py
index 8c1265211b87da7d006b9020087ce04994889fd0..345c77ef0fe2c3f0009d83082b106e8eb039d3c2 100644
--- a/workflow/tests/test_motif_search.py
+++ b/workflow/tests/test_motif_search.py
@@ -21,6 +21,11 @@ def test_limited_peaks_singleend():
def test_motif_search_singleend():
assert os.path.exists(os.path.join(test_output_path, 'ENCSR238SGC_memechip', 'ENCSR238SGC.fa'))
assert os.path.exists(os.path.join(test_output_path, 'ENCSR238SGC_memechip', 'index.html'))
+ assert os.path.getsize(os.path.join(test_output_path, 'ENCSR238SGC_memechip', 'ENCSR238SGC.fa')) > 0
+
+@pytest.mark.singleskip_true
+def test_motif_search_singleend():
+ assert os.path.isdir(test_output_path) == False
@pytest.mark.pairedend
@@ -29,8 +34,9 @@ def test_limited_peaks_pairedend():
assert os.path.exists(peak_file_ENCSR729LGA)
assert utils.count_lines(peak_file_ENCSR729LGA) == 600
-
+
@pytest.mark.pairedend
def test_motif_search_pairedend():
assert os.path.exists(os.path.join(test_output_path, 'ENCSR729LGA_memechip', 'ENCSR729LGA.fa'))
assert os.path.exists(os.path.join(test_output_path, 'ENCSR729LGA_memechip', 'index.html'))
+ assert os.path.getsize(os.path.join(test_output_path, 'ENCSR729LGA_memechip', 'ENCSR729LGA.fa')) > 0
diff --git a/workflow/tests/test_overlap_peaks.py b/workflow/tests/test_overlap_peaks.py
index 9289d165ee8a7cf3bb613616a165e852ad80d462..c450551ab5ac582977163c5d59220f5e807615c6 100644
--- a/workflow/tests/test_overlap_peaks.py
+++ b/workflow/tests/test_overlap_peaks.py
@@ -37,11 +37,11 @@ def test_check_update_design(design_diff):
def test_overlap_peaks_singleend():
assert os.path.exists(os.path.join(test_output_path, 'ENCSR238SGC.rejected.narrowPeak'))
peak_file = test_output_path + 'ENCSR238SGC.replicated.narrowPeak'
- assert utils.count_lines(peak_file) == 152848
+ assert utils.count_lines(peak_file) >= 149291
@pytest.mark.pairedend
def test_overlap_peaks_pairedend():
assert os.path.exists(os.path.join(test_output_path, 'ENCSR729LGA.rejected.narrowPeak'))
peak_file = test_output_path + 'ENCSR729LGA.replicated.narrowPeak'
- assert utils.count_lines(peak_file) == 26281
+ assert utils.count_lines(peak_file) >= 25657
diff --git a/workflow/tests/test_plot_profile.py b/workflow/tests/test_plot_profile.py
new file mode 100644
index 0000000000000000000000000000000000000000..6c9605d6d654f66adec865372223e13ddeaf6b19
--- /dev/null
+++ b/workflow/tests/test_plot_profile.py
@@ -0,0 +1,18 @@
+#!/usr/bin/env python3
+
+import pytest
+import os
+import utils
+
+test_output_path = os.path.dirname(os.path.abspath(__file__)) + \
+ '/../output/experimentQC/'
+
+
+@pytest.mark.singleend
+def test_plot_singleend():
+ assert os.path.exists(os.path.join(test_output_path, 'plotProfile.png'))
+
+
+@pytest.mark.pairedend
+def test_plot_pairedend():
+ assert os.path.exists(os.path.join(test_output_path, 'computeMatrix.gz'))
diff --git a/workflow/tests/test_pool_and_psuedoreplicate.py b/workflow/tests/test_pool_and_psuedoreplicate.py
index 31fffc57e7eaa598a933c97c92b1c901ba731ba7..f251e888f906de1c7e360502cadf45c8c9b85a3c 100644
--- a/workflow/tests/test_pool_and_psuedoreplicate.py
+++ b/workflow/tests/test_pool_and_psuedoreplicate.py
@@ -5,16 +5,18 @@ import pandas as pd
from io import StringIO
import os
import pool_and_psuedoreplicate
+import shutil
+test_design_path = os.path.dirname(os.path.abspath(__file__)) + \
+ '/../../test_data/'
test_output_path = os.path.dirname(os.path.abspath(__file__)) + \
'/../output/design/'
-DESIGN_STRING = """sample_id\ttag_align\txcor\tbiosample\tfactor\ttreatment\treplicate\tcontrol_tag_align
-A_1\tA_1.bedse.gz\tA_1.cc.qc\tLiver\tH3K27ac\tNone\t1\tB_1.bedse.gz
-A_2\tA_2.bedse.gz\tA_2.cc.qc\tLiver\tH3K27ac\tNone\t2\tB_2.bedse.gz
+DESIGN_STRING = """sample_id\tse_tag_align\ttag_align\txcor\texperiment_id\tbiosample\tfactor\ttreatment\treplicate\tcontrol_id\tcontrol_tag_align
+A_1\tA_1.tagAlign.gz\tA_1.bedse.gz\tA_1.cc.qc\tA\tLiver\tH3K27ac\tNone\t1\tB_1\tB_1.bedse.gz
+A_2\tA_2.tagAlign.gz\tA_2.bedse.gz\tA_2.cc.qc\tA\tLiver\tH3K27ac\tNone\t2\tB_2\tB_2.bedse.gz
"""
-
@pytest.fixture
def design_experiment():
design_file = StringIO(DESIGN_STRING)
@@ -60,6 +62,16 @@ def test_check_controls_single(design_experiment_3):
assert no_controls == 1
+@pytest.mark.unit
+def test_single_rep(design_experiment_2):
+ cwd = os.getcwd()
+ shutil.copy(test_design_path + 'A_1.bedse.gz', cwd)
+ shutil.copy(test_design_path + 'B_1.bedse.gz', cwd)
+ shutil.copy(test_design_path + 'A_1.tagAlign.gz', cwd)
+ shutil.copy(test_design_path + 'B_1.tagAlign.gz', cwd)
+ single_rep = pool_and_psuedoreplicate.generate_design('false', 1.2, design_experiment_2, cwd, 1, 1)
+ assert single_rep.shape[0] == 4
+
@pytest.mark.singleend
def test_pool_and_psuedoreplicate_singleend():
design_file = os.path.join(test_output_path, 'ENCSR238SGC_ppr.tsv')
diff --git a/workflow/tests/test_trim_reads.py b/workflow/tests/test_trim_reads.py
index f929d9c20cb95c3c89fdda22096870120195972d..e01707687b244fc3e7bf1b1b7fbd53797a43afbf 100644
--- a/workflow/tests/test_trim_reads.py
+++ b/workflow/tests/test_trim_reads.py
@@ -13,7 +13,7 @@ test_output_path = os.path.dirname(os.path.abspath(__file__)) + \
@pytest.mark.singleend
def test_trim_reads_singleend():
raw_fastq = test_data_path + 'ENCFF833BLU.fastq.gz'
- trimmed_fastq = test_output_path + 'ENCLB144FDT_R1_trimmed.fq.gz'
+ trimmed_fastq = test_output_path + 'ENCLB144FDT/ENCLB144FDT_R1_trimmed.fq.gz'
assert os.path.getsize(raw_fastq) != os.path.getsize(trimmed_fastq)
assert os.path.getsize(trimmed_fastq) == 2512853101
@@ -21,14 +21,14 @@ def test_trim_reads_singleend():
@pytest.mark.singleend
def test_trim_report_singleend():
trimmed_fastq_report = test_output_path + \
- 'ENCLB144FDT_R1_trimmed.fq.gz_trimming_report.txt'
+ 'ENCLB144FDT/ENCLB144FDT_R1.fastq.gz_trimming_report.txt'
assert 'Trimming mode: single-end' in open(trimmed_fastq_report).readlines()[4]
@pytest.mark.pairedend
def test_trim_reads_pairedend():
raw_fastq = test_data_path + 'ENCFF582IOZ.fastq.gz'
- trimmed_fastq = test_output_path + 'ENCLB637LZP_R2_val_2.fq.gz'
+ trimmed_fastq = test_output_path + 'ENCLB637LZP/ENCLB637LZP_R2_val_2.fq.gz'
assert os.path.getsize(raw_fastq) != os.path.getsize(trimmed_fastq)
assert os.path.getsize(trimmed_fastq) == 2229312710
@@ -36,5 +36,5 @@ def test_trim_reads_pairedend():
@pytest.mark.pairedend
def test_trim_report_pairedend():
trimmed_fastq_report = test_output_path + \
- 'ENCLB637LZP_R2.fastq.gz_trimming_report.txt'
+ 'ENCLB637LZP/ENCLB637LZP_R2.fastq.gz_trimming_report.txt'
assert 'Trimming mode: paired-end' in open(trimmed_fastq_report).readlines()[4]
diff --git a/workflow/tests/test_xcor.py b/workflow/tests/test_xcor.py
index 259ecfb1f82ce80e79ee14539c4bb603bb5ddfb4..19777a0be7f045cfb91ca3f1d96fdba7bd39322c 100644
--- a/workflow/tests/test_xcor.py
+++ b/workflow/tests/test_xcor.py
@@ -10,27 +10,27 @@ test_output_path = os.path.dirname(os.path.abspath(__file__)) + \
@pytest.mark.singleend
def test_cross_plot_singleend():
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB144FDT.cc.plot.pdf'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCLB144FDT/ENCLB144FDT.cc.plot.pdf'))
@pytest.mark.singleend
def test_cross_qc_singleend():
- qc_file = os.path.join(test_output_path,"ENCLB144FDT.cc.qc")
+ qc_file = os.path.join(test_output_path,"ENCLB144FDT/ENCLB144FDT.cc.qc")
df_xcor = pd.read_csv(qc_file, sep="\t", header=None)
- assert df_xcor[2].iloc[0] == '190,200,210'
- assert df_xcor[8].iloc[0] == 1.025906
- assert round(df_xcor[9].iloc[0], 6) == 1.139671
+ assert df_xcor[2].iloc[0] == '185,195,205'
+ assert df_xcor[8].iloc[0] == 1.02454
+ assert df_xcor[9].iloc[0] == 0.8098014
@pytest.mark.pairedend
def test_cross_qc_pairedend():
- assert os.path.exists(os.path.join(test_output_path, 'ENCLB568IYX.cc.plot.pdf'))
+ assert os.path.exists(os.path.join(test_output_path, 'ENCLB568IYX/ENCLB568IYX.cc.plot.pdf'))
@pytest.mark.pairedend
def test_cross_plot_pairedend():
- qc_file = os.path.join(test_output_path,"ENCLB568IYX.cc.qc")
+ qc_file = os.path.join(test_output_path,"ENCLB568IYX/ENCLB568IYX.cc.qc")
df_xcor = pd.read_csv(qc_file, sep="\t", header=None)
- assert df_xcor[2].iloc[0] == '220,430,475'
- assert round(df_xcor[8].iloc[0],6) == 1.060018
- assert df_xcor[9].iloc[0] == 4.099357
+ assert df_xcor[2].iloc[0] == '215,225,455'
+ assert round(df_xcor[8].iloc[0],6) == 1.056201
+ assert df_xcor[9].iloc[0] == 3.599357