Workflow: Trim Galore ChIP-Seq pipeline single-read
. This ChIP-Seq pipeline is based on the original [BioWardrobe's](https://biowardrobe.com) [PubMed ID:26248465](https://www.ncbi.nlm.nih.gov/pubmed/26248465) **ChIP-Seq** basic analysis workflow for a **single-read** experiment with Trim Galore. ### Data Analysis SciDAP starts from the .fastq files which most DNA cores and commercial NGS companies return. Starting from raw data allows us to ensure that all experiments have been processed in the same way and simplifies the deposition of data to GEO upon publication. The data can be uploaded from users computer, downloaded directly from an ftp server of the core facility by providing a URL or from GEO by providing SRA accession number. Our current pipelines include the following steps: 1. Trimming the adapters with TrimGalore. This step is particularly important when the reads are long and the fragments are short-resulting in sequencing adapters at the end of read. If adapter is not removed the read will not map. TrimGalore can recognize standard adapters, such as Illumina or Nexterra/Tn5 adapters. 2. QC 3. (Optional) trimming adapters on 5' or 3' end by the specified number of bases. 4. Mapping reads with BowTie. Only uniquely mapped reads with less than 3 mismatches are used in the downstream analysis. Results are saved as a .bam file. 5. (Optional) Removal of duplicates (reads/pairs of reads mapping to exactly same location). This step is used to remove reads overamplified in PCR. Unfortunately, it may also remove \"good\" reads. We usually do not remove duplicates unless the library is heavily duplicated. Please note that MACS2 will remove 'excessive' duplicates during peak calling ina smart way (those not supported by other nearby reads). 6. Peakcalling by MACS2. (Optionally), it is possible to specify read extension length for MACS2 to use if the length determined automatically is wrong. 7. Generation of BigWig coverage files for display on the browser. The coverage shows the number of fragments at each base in the genome normalized to the number of millions of mapped reads. In the case of PE sequencing the fragments are real, but in the case of single reads the fragments are estimated by extending reads to the average fragment length found by MACS2 or specified by the user in 6. ### Details _Trim Galore_ is a wrapper around [Cutadapt](https://github.com/marcelm/cutadapt) and [FastQC](http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) to consistently apply adapter and quality trimming to FastQ files, with extra functionality for RRBS data. In outputs it returns coordinate sorted BAM file alongside with index BAI file, quality statistics of the input FASTQ file, reads coverage in a form of BigWig file, peaks calling data in a form of narrowPeak or broadPeak files, islands with the assigned nearest genes and region type, data for average tag density plot (on the base of BAM file). Workflow starts with step *fastx\_quality\_stats* from FASTX-Toolkit to calculate quality statistics for input FASTQ file. At the same time `bowtie` is used to align reads from input FASTQ file to reference genome *bowtie\_aligner*. The output of this step is unsorted SAM file which is being sorted and indexed by `samtools sort` and `samtools index` *samtools\_sort\_index*. Based on workflow’s input parameters indexed and sorted BAM file can be processed by `samtools rmdup` *samtools\_rmdup* to get rid of duplicated reads. If removing duplicates is not required the original input BAM and BAI files return. Otherwise step *samtools\_sort\_index\_after\_rmdup* repeat `samtools sort` and `samtools index` with BAM and BAI files. Right after that `macs2 callpeak` performs peak calling *macs2\_callpeak*. On the base of returned outputs the next step *macs2\_island\_count* calculates the number of islands and estimated fragment size. If the last one is less that 80bp (hardcoded in the workflow) `macs2 callpeak` is rerun again with forced fixed fragment size value (*macs2\_callpeak\_forced*). If at the very beginning it was set in workflow input parameters to force run peak calling with fixed fragment size, this step is skipped and the original peak calling results are saved. In the next step workflow again calculates the number of islands and estimates fragment size (*macs2\_island\_count\_forced*) for the data obtained from *macs2\_callpeak\_forced* step. If the last one was skipped the results from *macs2\_island\_count\_forced* step are equal to the ones obtained from *macs2\_island\_count* step. Next step (*macs2\_stat*) is used to define which of the islands and estimated fragment size should be used in workflow output: either from *macs2\_island\_count* step or from *macs2\_island\_count\_forced* step. If input trigger of this step is set to True it means that *macs2\_callpeak\_forced* step was run and it returned different from *macs2\_callpeak* step results, so *macs2\_stat* step should return [fragments\_new, fragments\_old, islands\_new], if trigger is False the step returns [fragments\_old, fragments\_old, islands\_old], where sufix \"old\" defines results obtained from *macs2\_island\_count* step and sufix \"new\" - from *macs2\_island\_count\_forced* step. The following two steps (*bamtools\_stats* and *bam\_to\_bigwig*) are used to calculate coverage on the base of input BAM file and save it in BigWig format. For that purpose bamtools stats returns the number of mapped reads number which is then used as scaling factor by bedtools genomecov when it performs coverage calculation and saves it in BED format. The last one is then being sorted and converted to BigWig format by bedGraphToBigWig tool from UCSC utilities. Step *get\_stat* is used to return a text file with statistics in a form of [TOTAL, ALIGNED, SUPRESSED, USED] reads count. Step *island\_intersect* assigns genes and regions to the islands obtained from *macs2\_callpeak\_forced*. Step *average\_tag\_density* is used to calculate data for average tag density plot on the base of BAM file.
- Selected
- |
- Default Values
- Nested Workflows
- Tools
- Inputs/Outputs
Inputs
ID | Type | Title | Doc |
---|---|---|---|
threads | Integer (Optional) | Number of threads |
Number of threads for those steps that support multithreading |
broad_peak | Boolean (Optional) | Callpeak broad |
Set to call broad peak for MACS2 |
fastq_file | File [FASTQ] | FASTQ input file |
Reads data in a FASTQ format, received after single end sequencing |
clip_3p_end | Integer (Optional) | Clip from 3p end |
Number of bases to clip from the 3p end |
clip_5p_end | Integer (Optional) | Clip from 5p end |
Number of bases to clip from the 5p end |
genome_size | String | Effective genome size |
MACS2 effective genome size: hs, mm, ce, dm or number, for example 2.7e9 |
chrom_length | File [Textual format] | Chromosomes length file |
Chromosomes length file |
control_file | File (Optional) [BAM] | Use experiment as a control |
Use experiment as a control for MACS2 peak calling |
promoter_dist | Integer (Optional) | Max distance from gene TSS (in both direction) overlapping which the peak will be assigned to the promoter region |
Max distance from gene TSS (in both direction) overlapping which the peak will be assigned to the promoter region |
upstream_dist | Integer (Optional) | Max distance from the promoter (only in upstream direction) overlapping which the peak will be assigned to the upstream region |
Max distance from the promoter (only in upstream direction) overlapping which the peak will be assigned to the upstream region |
indices_folder | Directory | Indexed genome folder (bowtie) |
Path to indexed genome folder by **bowtie** |
annotation_file | File [TSV] | Annotation file |
Tab-separated annotation file |
exp_fragment_size | Integer (Optional) | Expected fragment size |
Expected fragment size for MACS2 |
remove_duplicates | Boolean (Optional) | Remove duplicates |
Calls samtools rmdup to remove duplicates from sortesd BAM file |
force_fragment_size | Boolean (Optional) | Force fragment size |
Force MACS2 to use exp_fragment_size |
Steps
ID | Runs | Label | Doc |
---|---|---|---|
preseq |
../tools/preseq-lc-extrap.cwl
(CommandLineTool)
|
Tool runs preseq lc_extrap. Only BAM input file is supported (-B option is used by default) successCodes: [1] - is used to pass this tool as a step in a workflow in case the BAM file was not correct for Preseq Discarded arguments: -V, -vals input is a text file containing only the observed counts -H, -hist input is a text file containing the observed histogram |
|
rename |
../tools/rename.cwl
(CommandLineTool)
|
Tool renames `source_file` to `target_filename`. Input `target_filename` should be set as string. If it's a full path, only basename will be used. If BAI file is present, it will be renamed too |
|
get_stat |
../tools/collect-statistics-chip-seq.cwl
(CommandLineTool)
|
Tool processes and combines log files generated by Trimgalore, Bowtie, Samtools and MACS2. |
|
trim_fastq |
../tools/trimgalore.cwl
(CommandLineTool)
|
Tool runs Trimgalore - the wrapper around Cutadapt and FastQC to consistently apply adapter and quality trimming
to FastQ files. |
|
bypass_trim |
../tools/bypass-trimgalore-se.cwl
(CommandLineTool)
|
If the number of reads in the trimmed_fastq_file is less then min_reads_count, tool will return original_fastq_file and null as selected_report_file. Otherwise, the trimmed_fastq_file and trimming_report_file will be returned. Might be usefull in case of trimgalore removed all reads from the original_fastq_file |
|
bam_to_bigwig |
../tools/bam-bedgraph-bigwig.cwl
(Workflow)
|
Workflow converts input BAM file into bigWig and bedGraph files. |
|
extract_fastq |
../tools/extract-fastq.cwl
(CommandLineTool)
|
Tool to decompress input FASTQ file(s). If several FASTQ files are provided, they will be concatenated in the order that corresponds to files in input. Bash script's logic: - disable case sensitive glob check - check if root name of input file already include '.fastq' or '.fq' extension. If yes, set DEFAULT_EXT to \"\", otherwise use '.fastq' - check file type, decompress if needed - return 1, if file type is not recognized This script also works of input file doesn't have any extension at all |
|
bowtie_aligner |
../tools/bowtie-alignreads.cwl
(CommandLineTool)
|
Tool maps input raw reads files to reference genome using Bowtie. |
|
macs2_callpeak |
../tools/macs2-callpeak-biowardrobe-only.cwl
(CommandLineTool)
|
Tool runs peak calling using MACS2 within the custom script `run.py` |
|
samtools_rmdup |
../tools/samtools-rmdup.cwl
(CommandLineTool)
|
Tool to remove duplicates from coordinate sorted BAM file set as input `bam_file`.
If input `trigger` is set to `true` or isn't set at all (`true` is used by default), run `samtools rmdup`, return
newly generated BAM file without duplicates and log as outputs `rmdup_output` and `rmdup_log`.
If input `trigger` is set to `false`, return unchanged BAM and index files (if provided in secondaryFiles),
previously staged into output directory, and log. |
|
island_intersect |
../tools/iaintersect.cwl
(CommandLineTool)
|
Tool assigns each peak obtained from MACS2 to a gene and region (upstream, promoter, exon, intron, intergenic) |
|
get_bam_statistics |
../tools/samtools-stats.cwl
(CommandLineTool)
|
Generates statistics for the input BAM file. |
|
average_tag_density |
../tools/atdp.cwl
(CommandLineTool)
|
Tool calculates average tag density profile around all annotated TSS. |
|
fastx_quality_stats |
../tools/fastx-quality-stats.cwl
(CommandLineTool)
|
Tool calculates statistics on the base of FASTQ file quality scores. If `output_filename` is not provided call function `default_output_filename` to return default output file name generated as `input_file` basename + `.fastxstat` extension. |
|
samtools_sort_index |
../tools/samtools-sort-index.cwl
(CommandLineTool)
|
Tool to sort and index input BAM/SAM/CRAM.
If input `trigger` is set to `true` or isn't set at all (`true` is used by default), run `samtools sort` and
`samtools index`, return sorted BAM and BAI/CSI index file.
If input `trigger` is set to `false`, return unchanged `sort_input` (BAM/SAM/CRAM) and index (BAI/CSI, if provided in
`secondaryFiles`) files, previously staged into output directory. |
|
samtools_sort_index_after_rmdup |
../tools/samtools-sort-index.cwl
(CommandLineTool)
|
Tool to sort and index input BAM/SAM/CRAM.
If input `trigger` is set to `true` or isn't set at all (`true` is used by default), run `samtools sort` and
`samtools index`, return sorted BAM and BAI/CSI index file.
If input `trigger` is set to `false`, return unchanged `sort_input` (BAM/SAM/CRAM) and index (BAI/CSI, if provided in
`secondaryFiles`) files, previously staged into output directory. |
|
get_bam_statistics_after_filtering |
../tools/samtools-stats.cwl
(CommandLineTool)
|
Generates statistics for the input BAM file. |
Outputs
ID | Type | Label | Doc |
---|---|---|---|
bigwig | File [bigWig] | BigWig file |
Generated BigWig file |
atdp_log | File [TSV] | ATDP log |
Average Tag Density generated log |
macs2_log | File (Optional) [Textual format] | MACS2 log |
MACS2 output log |
bowtie_log | File [Textual format] | BOWTIE alignment log |
BOWTIE generated alignment log |
atdp_result | File [TSV] | ATDP results |
Average Tag Density generated results |
bambai_pair | File [BAM] | Coordinate sorted BAM alignment file (+index BAI) |
Coordinate sorted BAM file and BAI index file |
trim_report | File | TrimGalore report |
TrimGalore generated log |
get_stat_log | File (Optional) [YAML] | YAML formatted combined log |
YAML formatted combined log |
macs2_moder_r | File (Optional) [Textual format] | MACS2 generated R script |
R script to produce a PDF image about the model based on your data |
iaintersect_log | File [TSV] | Island intersect log |
Iaintersect generated log |
unaligned_fastq | File[] (Optional) [FASTQ] | Unaligned FASTQ file(s) |
Unaligned FASTQ file(s) |
fastx_statistics | File [Textual format] | FASTQ statistics |
fastx_quality_stats generated FASTQ file quality statistics file |
preseq_estimates | File (Optional) [TSV] | Preseq estimates |
Preseq estimated results |
get_stat_markdown | File (Optional) [TIDE TXT] | Markdown formatted combined log |
Markdown formatted combined log |
macs2_broad_peaks | File (Optional) [ENCODE broad peak format] | Broad peaks |
Contains the peak locations together with peak summit, pvalue and qvalue |
macs2_gapped_peak | File (Optional) [bed12] | Gapped peaks |
Contains both the broad region and narrow peaks |
multimapped_fastq | File[] (Optional) [FASTQ] | Multimapped FASTQ file(s) |
Multimapped FASTQ file(s) |
iaintersect_result | File [TSV] | Island intersect results |
Iaintersect generated results |
macs2_called_peaks | File (Optional) [xls] | Called peaks |
XLS file to include information about called peaks |
macs2_narrow_peaks | File (Optional) [ENCODE narrow peak format] | Narrow peaks |
Contains the peak locations together with peak summit, pvalue and qvalue |
macs2_peak_summits | File (Optional) [BED] | Peak summits |
Contains the peak summits locations for every peaks |
samtools_rmdup_log | File [Textual format] | Remove duplicates log |
Samtools rmdup generated log |
macs2_fragment_stat | File (Optional) [Textual format] | FRAGMENT, FRAGMENTE, ISLANDS |
fragment, calculated fragment, islands count from MACS2 results |
mapped_reads_number | Integer | Mapped reads number |
Mapped reads number for downstream analyses |
bam_statistics_report | File [Textual format] | BAM statistics report (original) |
BAM statistics report (right after alignment and sorting) |
get_stat_formatted_log | File (Optional) [TSV] | Bowtie & Samtools Rmdup combined formatted log |
Processed and combined Bowtie aligner and Samtools rmdup formatted log |
estimated_fragment_size | Integer | Estimated fragment size |
Estimated fragment size for downstream analyses |
bam_statistics_report_after_filtering | File [Textual format] | BAM statistics report (after filtering) |
BAM statistics report (after all filters applied) |
https://w3id.org/cwl/view/git/8a92669a566589d80fde9d151054ffc220ed4ddd/workflows/trim-chipseq-se.cwl