Workflow: Trim Galore ATAC-Seq pipeline single-read

Fetched 2023-07-20 16:35:57 GMT

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. The pipeline was adapted for ATAC-Seq single-read data analysis by updating genome coverage step. _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. To adapt the pipeline for ATAC-Seq data analysis we calculate genome coverage using only the first 9 bp from every read. 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.

children parents
Workflow as SVG
  • 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 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

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

exclude_chromosome String (Optional) Exclude chromosomes

Space separated list of chromosomes to be excluded

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/python-get-stat-chipseq.cwl (CommandLineTool)
python-get-stat-rnaseq

Processing log files generated by BioWardrobe's ChIP-Seq workflows to produce expiriment quality table.

filter_bam
../tools/samtools-filter.cwl (CommandLineTool)

Excludes chromosomes from the input BAM file

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.

`default_log_name` function returns names for generated log files (for both paired-end and single-end cases). `trim_galore` itself doesn't support setting custom names for output files.

For paired-end data processing both `input_file_pair` and `paired` should be set. If either of them is not set, the other one becomes unset automatically.

bam_to_bigwig

Workflow converts input BAM file into bigWig and bedGraph files

extract_fastq
../tools/extract-fastq.cwl (CommandLineTool)

Tool to decompress input FASTQ file 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 \"\" - 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.

`default_output_filename` function returns default name for SAM output and log files. In case when `sam` and `output_filename` inputs are not set, default filename will have `.sam` extension but format may not correspond SAM specification. To set output filename manually use `output_filename` input. Default output filename is based on `output_filename` or basename of `upstream_filelist`, `downstream_filelist` or `crossbow_filelist` file (if array, the first file in array is taken). If function is called without argenments and `output_filename` input is set, it will be returned from the function.

For single-end input data any of the `upstream_filelist` or `downstream_filelist` inputs can be used.

Log filename (`log_file` output) is generated by `default_output_filename` function with ex='.bw'

`indices_folder` defines folder to contain Bowtie indices. Based on the first found file with `rev.1.ebwt` or `rev.1.ebwtl` extension, bowtie index prefix is returned from input's `valueFrom` field.

macs2_callpeak
../tools/macs2-callpeak-biowardrobe-only.cwl (CommandLineTool)

Tool runs peak calling using MACS2 within the custom script `run.py`

`default_output_filename` function returns default output prefix for all generated output files. Function is called when `output_prefix` input is not set. Optional argument `sufix` is equal to `_macs2` by default. `sufix` is hardcoded and should not be changed (if necessary, update `_macs2` in all tool outputs and `run.py` script too). `sufix` is included in the filenames of all output files except those, who ends with `_fragment_stat.tsv`. For the rest of the output files `_macs2` is already included as a part of the `sufix` argument when `default_output_filename` function is called. This is mostly done because of the naming conventions for output files in supplemental project `https://github.com/Barski-lab/biowardrobe-analysis` (if `_fragment_stat.tsv` is updated to `_macs2_fragment_stat.tsv` in `biowardrobe-analysis`, the way we use `default_output_filename` function can be simplified). The basename for output filename is taken from `treatment_file` input (if array - the first item is used).

`baseCommand` runs Python script `run.py` staged into the output directory. Within this script MACS2 is called one ortwo times as a subprocess with all provided input parameters.

`extsize` is mandatory parameter. It will be ignored by MASC2, when `--nomodel` is not set.

The aim of `run.py` script is to return `calculated_fragment_size`, `expected_fragment_size` and `island_count` as a text file.

Scenarios: 1. Run MACS2 with all input parameters. If MACS2 first run failed and `--nomodel` wasn't set (it means that `--extsize` wasn't used), set `expected_fragment_size` to `--extsize` value and try to rerun MACS2 with `--nomodel` parameter. If MACS2 second run failed too - exit workflow, if it didn't fail - get values for `calculated_fragment_size` and `island_count` from generated output file. Returned data: `calculated_fragment_size` - the value of fragment size from MACS2 second run `expected_fragment_size` - the value of fragment size equal to `--extsize` (set when MACS2 first run failed) `island_count` - the number of islands from MACS2 second run

2. Run MACS2 with all input parameters. If MACS2 first run failed and `--nomodel` was already set, exit workflow, because there is no need to rerun MACS2, it will fail too.

3. Run MACS2 with all input parameters. If MACS2 first run didn't fail, we get `calculated_fragment_size` and `island_count` from generated output file, and set `expected_fragment_size` to be equal to `calculated_fragment_size`. If `calculated_fragment_size` is bigger or equal to 80, there is no need to rerun MACS2. Returned data: `calculated_fragment_size` - the value of fragment size from MACS2 first run `expected_fragment_size` - the value of fragment size equal to `calculated_fragment_size` `island_count` - the number of islands from MACS2 first run

4. Run MACS2 with all input parameters. If MACS2 first run didn't fail, we get `calculated_fragment_size` and `island_count` from generated output file, and set `expected_fragment_size` to be equal to `calculated_fragment_size`. If `calculated_fragment_size` is less then 80 and MACS2 was run without '--nomodel' (it means we didn't force it to use the fixed fragment size `--extsize`), we rerun MACS2 with `--nomodel` argument. If MACS2 second run failed too - exit workflow, if not - get values for `calculated_fragment_size` and `island_count` from generated output file. Returned data: `calculated_fragment_size` - the value of fragment size from MACS2 second run `expected_fragment_size` - the value of fragment size from MACS2 first run `island_count` - the number of islands from MACS2 second run

5. Run MACS2 with all input parameters. If MACS2 first run didn't fail and was already run with '--nomodel' parameter (forced to use fixed fragment size `--extsize`), we get `calculated_fragment_size` and `island_count` from generated output file, and set `expected_fragment_size` to be equal to `calculated_fragment_size`. No mater which value we got for `calculated_fragment_size`, we don't rerun MACS2. Returned data: `calculated_fragment_size` - the value of fragment size got from MACS2 first run `expected_fragment_size` - the value of fragment size equal to `calculated_fragment_size` `island_count` - the number of islands got from MACS2 first run of MACS2

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.

Before execution `baseCommand` input BAM and index files (if provided in secondaryFiles) are staged into directory set as docker parameter `--workdir` (tool's output directory), using `InitialWorkDirRequirement`. Setting `writable: true` makes cwl-runner to make copies of the BAM and index files (if provided in secondaryFiles) and mount them to docker container with `rw` mode as part of `--workdir` (if set to false, the files staged into output directory will be mounted to docker container separately with `ro` mode). Because `bam_file` is copied into `--workdir` (where `samtools rmdup` should be run) we need to change its name - set `entryname`. This prevents `samtools rmdup` from reading and writing to the same file (which happens only when `output_filename` is not set and we generated output filename automatically on the base of the `bam_file` basename). Index file which is optionally set as `secondaryFiles` for `bam_file` shouln't be renamed, because `samtools rmdup` don't work with this file at all and is completely ignored when `trigger` is true.

Trigger logic is implemented in bash script set by default in input `bash_script`. If first argment $0 (which is `trigger` input) is true, run `samtools rmdup` with the rest of the arguments. If $0 is not true, skip samtools running and rename staged into output directory BAM file. Current filename of staged BAM file and the filename which we are waiting to be the output of tool are always set as two last arguments of the script.

Input `trigger` is Boolean, but returns String, because of `valueFrom` field. The `valueFrom` is used, because if `trigger` is false, cwl-runner doesn't append this argument at all to the the `baseCommand` - new feature of CWL v1.0.2. Alternatively, `prefix` field could be used, but it causes changing logic in bash script saved in `bash_script` input.

`default_output_filename` function is used for generating output filename if input `output_filename` is not set or in case when `trigger` is false and we need to return BAM and index files (if provided in secondaryFiles) staged into output directory.

Output `rmdup_output` returns `secondaryFiles` only in case when `trigger` was set to false (we need to rerun index).

island_intersect
../tools/iaintersect.cwl (CommandLineTool)

Tool assigns each peak obtained from MACS2 to a gene and region (upstream, promoter, exon, intron, intergenic)

`default_output_filename` function returns output filename with sufix set as `ext` argument. Function is called when either `output_filename` or `log_filename` inputs are not provided.

average_tag_density
../tools/atdp.cwl (CommandLineTool)

Tool calculates average tag density profile around all annotated TSS.

`default_output_filename` function returns output filename with sufix set as `ext` argument. Function is called when either `output_filename` or `log_filename` inputs are not provided.

Before running `baseCommand`, annotaion file `annotation_filename` is staged into output directory (Docker's `--workdir`) with `\"writable\": true` (to allow to overwrite it by `refgene-sort`).

To run `atdp` index file should be provided (either in `secondaryFiles` of `input_file` or as separate input `index_file`)

`baseCommand` runs bash script from `script` input. Script runs `refgene-sort` to sort annotaion file and then runs `atdp`. `refgene-sort` - utility to sort refgene annotation files, using MySQL syntax.

Optionally (with cwltool==1.0.20171107133715), script can be simplified to #!/bin/bash set -- \"$0\" \"$@\" refgene-sort -i \"${2:4}\" -o \"${2:4}\" -s \"ORDER BY chrom, strand, CASE strand WHEN '+' THEN txStart WHEN '-' THEN txEnd END\" atdp \"$@\" because `set -- \"$1\" --a=$(basename \"${2:4}\") \"${@:3}\"` is not needed anymore.

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.

Before execution `baseCommand`, `sort_input` and `secondaryFiles` (if provided) are staged into directory set as docker parameter `--workdir` (tool's output directory), using `InitialWorkDirRequirement`. Setting `writable: true` makes cwl-runner to make copies of the `sort_input` and `secondaryFiles` (if provided) and mount them to docker container with `rw` mode as part of `--workdir` (if set to false, the files staged into output directory will be mounted to docker container separately with `ro` mode). Because both `samtools sort` and `samtools index` can overwrite files with the same names (and in case of `samtools sort` even the input file can be overwritten), we don't need to rename any of the staged files.

Trigger logic is implemented in two bash scripts set by default as `bash_script_sort` and `bash_script_index` inputs. For both of then, if the first argument $0 (which is `trigger` input) is true, run `samtools sort/index` with the rest of the arguments. If $0 is not true, skip `samtools sort/index` and return `sort_input` and `secondaryFiles` (if provided) staged into output directory.

Input `trigger` is Boolean, but returns String, because of `valueFrom` field. The `valueFrom` is used, because if `trigger` is false, cwl-runner doesn't append this argument at all to the the `baseCommand` - new feature of CWL v1.0.2. Alternatively, `prefix` field could be used, but it causes changing in script logic.

If using `sort_output_filename`, the output file extension should be `*.bam`, because `samtools sort` defines the output file format on the base of the file extension. If `*.sam` is sed as output filename, it cannot be usefully indexed by `samtools index`.

`default_bam` function is used to generate output filename for `samtools sort` if input `sort_output_filename` is not set or when `trigger` is false and we need to return `sort_input` and `secondaryFiles` (if provided) files staged into output directory. Output filename is generated on the base of `sort_input` basename with `.bam` extension by default.

`ext` function is used to return the index file extension (BAI/CSI) based on `csi` and `bai` inputs according to the following logic `csi` && `bai` => BAI !`csi` && !`bai ` => BAI `csi` && !`bai ` => CSI

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.

Before execution `baseCommand`, `sort_input` and `secondaryFiles` (if provided) are staged into directory set as docker parameter `--workdir` (tool's output directory), using `InitialWorkDirRequirement`. Setting `writable: true` makes cwl-runner to make copies of the `sort_input` and `secondaryFiles` (if provided) and mount them to docker container with `rw` mode as part of `--workdir` (if set to false, the files staged into output directory will be mounted to docker container separately with `ro` mode). Because both `samtools sort` and `samtools index` can overwrite files with the same names (and in case of `samtools sort` even the input file can be overwritten), we don't need to rename any of the staged files.

Trigger logic is implemented in two bash scripts set by default as `bash_script_sort` and `bash_script_index` inputs. For both of then, if the first argument $0 (which is `trigger` input) is true, run `samtools sort/index` with the rest of the arguments. If $0 is not true, skip `samtools sort/index` and return `sort_input` and `secondaryFiles` (if provided) staged into output directory.

Input `trigger` is Boolean, but returns String, because of `valueFrom` field. The `valueFrom` is used, because if `trigger` is false, cwl-runner doesn't append this argument at all to the the `baseCommand` - new feature of CWL v1.0.2. Alternatively, `prefix` field could be used, but it causes changing in script logic.

If using `sort_output_filename`, the output file extension should be `*.bam`, because `samtools sort` defines the output file format on the base of the file extension. If `*.sam` is sed as output filename, it cannot be usefully indexed by `samtools index`.

`default_bam` function is used to generate output filename for `samtools sort` if input `sort_output_filename` is not set or when `trigger` is false and we need to return `sort_input` and `secondaryFiles` (if provided) files staged into output directory. Output filename is generated on the base of `sort_input` basename with `.bam` extension by default.

`ext` function is used to return the index file extension (BAI/CSI) based on `csi` and `bai` inputs according to the following logic `csi` && `bai` => BAI !`csi` && !`bai ` => BAI `csi` && !`bai ` => CSI

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) [Textual format] BioWardrobe compatible Bowtie & Samtools Rmdup combined log

BioWardrobe compatible processed and combined Bowtie aligner and Samtools rmdup 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

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

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

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

get_stat_formatted_log File (Optional) [TSV] Bowtie & Samtools Rmdup combined formatted log

Processed and combined Bowtie aligner and Samtools rmdup formatted log

Permalink: https://w3id.org/cwl/view/git/e238d1756f1db35571e84d72e1699e5d1540f10c/workflows/trim-atacseq-se.cwl