Explore Workflows

# ChIP-Seq basic analysis workflow for single-read data Reads are aligned to the reference genome with [Bowtie](http://bowtie-bio.sourceforge.net/index.shtml). Results are saved as coordinate sorted [BAM](http://samtools.github.io/hts-specs/SAMv1.pdf) alignment and index BAI files. Optionally, PCR duplicates can be removed. To obtain coverage in [bigWig](https://genome.ucsc.edu/goldenpath/help/bigWig.html) format, average fragment length is calculated by [MACS2](https://github.com/taoliu/MACS), and individual reads are extended to this length in the 3’ direction. Areas of enrichment identified by MACS2 are saved in ENCODE [narrow peak](http://genome.ucsc.edu/FAQ/FAQformat.html#format12) or [broad peak](https://genome.ucsc.edu/FAQ/FAQformat.html#format13) formats. Called peaks together with the nearest genes are saved in TSV format. In addition to basic statistics (number of total/mapped/multi-mapped/unmapped/duplicate reads), pipeline generates several quality control measures. Base frequency plots are used to estimate adapter contamination, a frequent occurrence in low-input ChIP-Seq experiments. Expected distinct reads count from [Preseq](http://smithlabresearch.org/software/preseq/) can be used to estimate read redundancy for a given sequencing depth. Average tag density profiles can be used to estimate ChIP enrichment for promoter proximal histone modifications. Use of different parameters for different antibodies (calling broad or narrow peaks) is possible. Additionally, users can elect to use BAM file from another experiment as control for MACS2 peak calling. ## Cite as *Kartashov AV, Barski A. BioWardrobe: an integrated platform for analysis of epigenomics and transcriptomics data. Genome Biol. 2015;16(1):158. Published 2015 Aug 7. [doi:10.1186/s13059-015-0720-3](https://www.ncbi.nlm.nih.gov/pubmed/26248465)* ## Software versions - Bowtie 1.2.0 - Samtools 1.4 - Preseq 2.0 - MACS2 2.1.1.20160309 - Bedtools 2.26.0 - UCSC userApps v358 ## Inputs | ID | Label | Description | Required | Default | Upstream analyses | | ------------------------- | ---------------------------------------------- | ---------------------------------------------------------------------------------------------------------------------------------------------------------------- | :------: | ------- | ------------------------------- | | **fastq\_file** | FASTQ file | Single-read sequencing data in FASTQ format (fastq, fq, bzip2, gzip, zip) | + | | | | **indices\_folder** | Genome indices | Directory with the genome indices generated by Bowtie | + | | genome\_indices/bowtie\_indices | | **annotation\_file** | Genome annotation file | Genome annotation file in TSV format | + | | genome\_indices/annotation | | **genome\_size** | Effective genome size | The length of the mappable genome (hs, mm, ce, dm or number, for example 2.7e9) | + | | genome\_indices/genome\_size | | **chrom\_length** | Chromosome lengths file | Chromosome lengths file in TSV format | + | | genome\_indices/chrom\_length | | **broad\_peak** | Call broad peaks | Make MACS2 call broad peaks by linking nearby highly enriched regions | + | | | | **control\_file** | Control ChIP-Seq single-read experiment | Indexed BAM file from the ChIP-Seq single-read experiment to be used as a control for MACS2 peak calling | | Null | control\_file/bambai\_pair | | **exp\_fragment\_size** | Expected fragment size | Expected fragment size for read extenstion towards 3' end if *force\_fragment\_size* was set to True or if calculated by MACS2 fragment size was less that 80 bp | | 150 | | | **force\_fragment\_size** | Force peak calling with expected fragment size | Make MACS2 don't build the shifting model and use expected fragment size for read extenstion towards 3' end | | False | | | **clip\_3p\_end** | Clip from 3' end | Number of base pairs to clip from 3' end | | 0 | | | **clip\_5p\_end** | Clip from 5' end | Number of base pairs to clip from 5' end | | 0 | | | **remove\_duplicates** | Remove PCR duplicates | Remove PCR duplicates from sorted BAM file | | False | | | **threads** | Number of threads | Number of threads for those steps that support multithreading | | 2 | | ## Outputs | ID | Label | Description | Required | Visualization | | ------------------------ | ---------------------------------- | ------------------------------------------------------------------------------------ | :------: | ------------------------------------------------------------------ | | **fastx\_statistics** | FASTQ quality statistics | FASTQ quality statistics in TSV format | + | *Base Frequency* and *Quality Control* plots in *QC Plots* tab | | **bambai\_pair** | Aligned reads | Coordinate sorted BAM alignment and index BAI files | + | *Nucleotide Sequence Alignments* track in *IGV Genome Browser* tab | | **bigwig** | Genome coverage | Genome coverage in bigWig format | + | *Genome Coverage* track in *IGV Genome Browser* tab | | **iaintersect\_result** | Gene annotated peaks | MACS2 peak file annotated with nearby genes | + | *Peak Coordinates* table in *Peak Calling* tab | | **atdp\_result** | Average Tag Density Plot | Average Tag Density Plot file in TSV format | + | *Average Tag Density Plot* in *QC Plots* tab | | **macs2\_called\_peaks** | Called peaks | Called peaks file with 1-based coordinates in XLS format | + | | | **macs2\_narrow\_peaks** | Narrow peaks | Called peaks file in ENCODE narrow peak format | | *Narrow peaks* track in *IGV Genome Browser* tab | | **macs2\_broad\_peaks** | Broad peaks | Called peaks file in ENCODE broad peak format | | *Broad peaks* track in *IGV Genome Browser* tab | | **preseq\_estimates** | Expected Distinct Reads Count Plot | Expected distinct reads count file from Preseq in TSV format | | *Expected Distinct Reads Count Plot* in *QC Plots* tab | | **workflow\_statistics** | Workflow execution statistics | Overall workflow execution statistics from bowtie\_aligner and samtools\_rmdup steps | + | *Overview* tab and experiment's preview | | **bowtie\_log** | Read alignment log | Read alignment log file from Bowtie | + | |

https://github.com/datirium/workflows.git

Path: workflows/chipseq-se.cwl

Branch/Commit ID: 954bb2f213d97dfef1cddaf9e830169a92ad0c6b

https://github.com/genome/analysis-workflows.git

Path: definitions/pipelines/germline_wgs.cwl

Branch/Commit ID: 40097e1ed094c5b42b68f3db2ff2cbe78c182479

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 **paired-end** experiment. A [FASTQ](http://maq.sourceforge.net/fastq.shtml) input file has to be provided. The pipeline produces a sorted BAM file alongside with index BAI file, quality statistics of the input FASTQ file, coverage by estimated fragments as a 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. 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 an unsorted SAM file which is being sorted and indexed by `samtools sort` and `samtools index` *samtools\_sort\_index*. Depending 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 BAM and BAI files are returned. Otherwise step *samtools\_sort\_index\_after\_rmdup* repeat `samtools sort` and `samtools index` with BAM and BAI files without duplicates. Next `macs2 callpeak` performs peak calling *macs2\_callpeak* and the next step reports *macs2\_island\_count* the number of islands and estimated fragment size. If the latter is less that 80bp (hardcoded in the workflow) `macs2 callpeak` is rerun again with forced fixed fragment size value (*macs2\_callpeak\_forced*). It is also possible to force MACS2 to use pre set fragment size in the first place. 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 from BAM file and save it in BigWig format. For that purpose bamtools stats returns the number of mapped reads which is then used as scaling factor by bedtools genomecov when it performs coverage calculation and saves it as a BEDgraph file whichis then 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 nearest 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 from the BAM file.

https://github.com/datirium/workflows.git

Path: workflows/chipseq-pe.cwl

Branch/Commit ID: 44214a9d02e6d85b03eb708552ed812ae3d4a733

https://github.com/common-workflow-language/cwltool.git

Path: tests/wf/conflict.cwl

Branch/Commit ID: 83038feb2a6fc3bab952e1ecc2a11bfbc8c557b4

Packed ID: main

https://github.com/Duke-GCB/bespin-cwl.git

Path: subworkflows/exomeseq-03-organizedirectories.cwl

Branch/Commit ID: 216ff9bf78130add564f7bcfba6385d5dab4c77d

Cross-Linking ImmunoPrecipitation ================================= `CLIP` (`cross-linking immunoprecipitation`) is a method used in molecular biology that combines UV cross-linking with immunoprecipitation in order to analyse protein interactions with RNA or to precisely locate RNA modifications (e.g. m6A). (Uhl|Houwaart|Corrado|Wright|Backofen|2017)(Ule|Jensen|Ruggiu|Mele|2003)(Sugimoto|König|Hussain|Zupan|2012)(Zhang|Darnell|2011) (Ke| Alemu| Mertens| Gantman|2015) CLIP-based techniques can be used to map RNA binding protein binding sites or RNA modification sites (Ke| Alemu| Mertens| Gantman|2015)(Ke| Pandya-Jones| Saito| Fak|2017) of interest on a genome-wide scale, thereby increasing the understanding of post-transcriptional regulatory networks. The identification of sites where RNA-binding proteins (RNABPs) interact with target RNAs opens the door to understanding the vast complexity of RNA regulation. UV cross-linking and immunoprecipitation (CLIP) is a transformative technology in which RNAs purified from _in vivo_ cross-linked RNA-protein complexes are sequenced to reveal footprints of RNABP:RNA contacts. CLIP combined with high-throughput sequencing (HITS-CLIP) is a generalizable strategy to produce transcriptome-wide maps of RNA binding with higher accuracy and resolution than standard RNA immunoprecipitation (RIP) profiling or purely computational approaches. The application of CLIP to Argonaute proteins has expanded the utility of this approach to mapping binding sites for microRNAs and other small regulatory RNAs. Finally, recent advances in data analysis take advantage of cross-link–induced mutation sites (CIMS) to refine RNA-binding maps to single-nucleotide resolution. Once IP conditions are established, HITS-CLIP takes ~8 d to prepare RNA for sequencing. Established pipelines for data analysis, including those for CIMS, take 3–4 d. Workflow -------- CLIP begins with the in-vivo cross-linking of RNA-protein complexes using ultraviolet light (UV). Upon UV exposure, covalent bonds are formed between proteins and nucleic acids that are in close proximity. (Darnell|2012) The cross-linked cells are then lysed, and the protein of interest is isolated via immunoprecipitation. In order to allow for sequence specific priming of reverse transcription, RNA adapters are ligated to the 3' ends, while radiolabeled phosphates are transferred to the 5' ends of the RNA fragments. The RNA-protein complexes are then separated from free RNA using gel electrophoresis and membrane transfer. Proteinase K digestion is then performed in order to remove protein from the RNA-protein complexes. This step leaves a peptide at the cross-link site, allowing for the identification of the cross-linked nucleotide. (König| McGlincy| Ule|2012) After ligating RNA linkers to the RNA 5' ends, cDNA is synthesized via RT-PCR. High-throughput sequencing is then used to generate reads containing distinct barcodes that identify the last cDNA nucleotide. Interaction sites can be identified by mapping the reads back to the transcriptome.

https://github.com/datirium/workflows.git

Path: workflows/clipseq-se.cwl

Branch/Commit ID: 9ee330737f4603e4e959ffe786fbb2046db70a00

A basic analysis workflow for paired-read CUT&RUN and CUT&TAG sequencing experiments. These sequencing library prep methods are ultra-sensitive chromatin mapping technologies compared to the ChIP-Seq methodology. Its primary benefits include 1) length filtering, 2) a higher signal-to-noise ratio, and 3) built-in normalization for between sample comparisons. This workflow utilizes the tool MACS2 which calls enriched regions in the target sequence data by identifying the top regions by area under a poisson distribution (of the alignment pileup). This workflow is loosely based on the [CUT-RUNTools-2.0 pipeline](https://github.com/fl-yu/CUT-RUNTools-2.0) pipeline, and the ChIP-Seq pipeline from [BioWardrobe](https://biowardrobe.com) [PubMed ID:26248465](https://www.ncbi.nlm.nih.gov/pubmed/26248465) was used as a CWL template. ### __Inputs__ *General Info (required\*):* - Experiment short name/Alias* - a unique name for the sample (e.g. what was used on tubes while processing it) - Cells* - sample cell type or organism name - Conditions* - experimental condition name - Catalog # - catalog number for cells from vender/supplier - Primary [genome index](https://scidap.com/tutorials/basic/genome-indices) for peak calling* - preprocessed genome index of sample organism for primary alignment and peak calling - Secondary [genome index](https://scidap.com/tutorials/basic/genome-indices) for spike-in normalization* - preprocessed genome index of spike-in organism for secondary alignment (of unaligned reads from primary alignment) and spike-in normalization, default should be E. coli K-12 - FASTQ file for R1* - read 1 file of a pair-end library - FASTQ file for R2* - read 2 file of a pair-end library *Advanced:* - - Number of bases to clip from the 3p end - used by bowtie aligner to trim <int> bases from 3' (right) end of reads - Number of bases to clip from the 5p end - used by bowtie aligner to trim <int> bases from 5' (left) end of reads - Call samtools rmdup to remove duplicates from sorted BAM file? - toggle on/off to remove duplicate reads from analysis - Fragment Length Filter will retain fragments between set base pair (bp) ranges for peak analysis - drop down menu - `default_below_1000` retains fragments <1000 bp - `histones_130_to_300` retains fragments between 130-300 bp - `TF_below_130` retains fragments <130 bp - Max distance (bp) from gene TSS (in both directions) overlapping which the peak will be assigned to the promoter region - default set to `1000` - Max distance (bp) from the promoter (only in upstream directions) overlapping which the peak will be assigned to the upstream region - default set to `20000` - Number of threads for steps that support multithreading - default set to `2` ### __Outputs__ Intermediate and final downloadable outputs include: - IGV with gene, BigWig (raw and normalized), and stringent peak tracks - quality statistics and visualizations for both R1/R2 input FASTQ files - coordinate sorted BAM file with associated BAI file for primary alignment - read pileup/coverage in BigWig format (raw and normalized) - cleaned bed files (containing fragment coordinates), and spike-in normalized peak-called BED files (also includes \"narrow\" and \"broad\" peaks). - stringent peak call bed file with nearest gene annotations per peak ### __Data Analysis 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 ends of reads. If adapter is not removed the read will not map. TrimGalore can recognize standard adapters, such as Illumina or Nextera/Tn5 adapters. 2. Generate quality control statistics of trimmed, unmapped sequence data 3. (Optional) Clipping of 5' and/or 3' end by the specified number of bases. 4. Mapping reads to primary genome index with Bowtie. - Only uniquely mapped reads with less than 3 mismatches are used in the downstream analysis. Results are then sorted and indexed. Final outputs are in bam/bai format, which are also used to extrapolate effects of additional sequencing based on library complexity. 5. (Optional) Removal of duplicates (reads/pairs of reads mapping to exactly the same location). - This step is used to remove reads overamplified during amplification of the library. Unfortunately, it may also remove \"good\" reads. We usually do not remove duplicates unless the library is heavily duplicated. 6. Mapping unaligned reads from primary alignment to secondary genome index with Bowtie. - This step is used to obtain the number of reads for normalization, used to scale the pileups from the primary alignment. After normalization, sample pileups/peak may then be appropriately compared to one another assuming an equal use of spike-in material during library preparation. Note the default genome index for this step should be *E. coli* K-12 if no spike-in material was called out in the library protocol. Refer to [Step 16](https://www.protocols.io/view/cut-amp-tag-data-processing-and-analysis-tutorial-e6nvw93x7gmk/v1?step=16#step-4A3D8C70DC3011EABA5FF3676F0827C5) of the \"CUT&Tag Data Processing and Analysis Tutorial\" by Zheng Y et al (2020). Protocol.io. 7. Formatting alignment file to account for fragments based on paired-end BAM. - Generates a filtered and normalized bed file to be used as input for peak calling. 8. Call enriched regions using MACS2. - This step called peaks (broad and narrow) using the MACS2 tool with default parameters and no normalization to a control sample. 9. Generation and formatting of output files. - This step collects read, alignment, and peak statistics, as well asgenerates BigWig coverage/pileup files for display on the browser using IGV. The coverage shows the number of fragments that cover each base in the genome both normalized and unnormalized to the calculated spike-in scaling factor. ### __References__ - Meers MP, Tenenbaum D, Henikoff S. (2019). Peak calling by Sparse Enrichment Analysis for CUT&RUN chromatin profiling. Epigenetics and Chromatin 12(1):42. - Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol 10:R25.

https://github.com/datirium/workflows.git

Path: workflows/cutandrun-macs2-pe.cwl

Branch/Commit ID: 22880e0f41d0420a17d643e8a6e8ee18165bbfbf

Generates tag density heatmap and histogram for the list of features in a headerless regions file. Inputs used are the bigWig file(s) of one or more ChIP/ATAC/C&R samples, and one or more filtered feature file(s) from the filtering and/or set operation workflows. The latter format contains `chrom start end name score strand`, only the first 3 columns are used in deeptools computeMatrix tool. The matrix is then used as input to plotHeatmap to generate the tag density plot and tag enrichment heatmap. computeMatrix paramters: --regionsFileName, -R File name, in BED format, containing the regions to plot. If multiple bed files are given, each one is considered a group that can be plotted separately. Also, adding a “#” symbol in the bed file causes all the regions until the previous “#” to be considered one group. --scoreFileName, -S bigWig file(s) containing the scores to be plotted. BigWig files can be obtained by using the bamCoverage or bamCompare tools. More information about the bigWig file format can be found at http://genome.ucsc.edu/goldenPath/help/bigWig.html --outFileName, -o File name to save the gzipped matrix file needed by the “plotHeatmap” and “plotProfile” tools. --beforeRegionStartLength=0, -b=0, --upstream=0 Distance upstream of the start site of the regions defined in the region file. If the regions are genes, this would be the distance upstream of the transcription start site. --regionBodyLength=1000, -m=1000 Distance in bases to which all regions will be fit. --afterRegionStartLength=0, -a=0, --downstream=0 Distance downstream of the end site of the given regions. If the regions are genes, this would be the distance downstream of the transcription end site. --numberOfProcessors=max/2, -p=max/2 Number of processors to use. Type “max/2” to use half the maximum number of processors or “max” to use all available processors. plotHeatmap parameters: --matrixFile, -m Matrix file from the computeMatrix tool. --outFileName, -out File name to save the image to. The file ending will be used to determine the image format. The available options are: “png”, “eps”, “pdf” and “svg”, e.g., MyHeatmap.png. --sortRegions=descend Whether the heatmap should present the regions sorted. The default is to sort in descending order based on the mean value per region. Possible choices: descend, ascend, no --sortUsing=mean Indicate which method should be used for sorting. For each row the method is computed. Possible choices: mean, median, max, min, sum, region_length --colorMap=RdYlBu Color map to use for the heatmap. Available values can be seen here: http://matplotlib.org/users/colormaps.html The available options are: ‘Spectral’, ‘summer’, ‘coolwarm’, ‘Set1’, ‘Set2’, ‘Set3’, ‘Dark2’, ‘hot’, ‘RdPu’, ‘YlGnBu’, ‘RdYlBu’, ‘gist_stern’, ‘cool’, ‘gray’, ‘GnBu’, ‘gist_ncar’, ‘gist_rainbow’, ‘CMRmap’, ‘bone’, ‘RdYlGn’, ‘spring’, ‘terrain’, ‘PuBu’, ‘spectral’, ‘gist_yarg’, ‘BuGn’, ‘bwr’, ‘cubehelix’, ‘YlOrRd’, ‘Greens’, ‘PRGn’, ‘gist_heat’, ‘Paired’, ‘hsv’, ‘Pastel2’, ‘Pastel1’, ‘BuPu’, ‘copper’, ‘OrRd’, ‘brg’, ‘gnuplot2’, ‘jet’, ‘gist_earth’, ‘Oranges’, ‘PiYG’, ‘YlGn’, ‘Accent’, ‘gist_gray’, ‘flag’, ‘BrBG’, ‘Reds’, ‘RdGy’, ‘PuRd’, ‘Blues’, ‘Greys’, ‘autumn’, ‘pink’, ‘binary’, ‘winter’, ‘gnuplot’, ‘RdBu’, ‘prism’, ‘YlOrBr’, ‘rainbow’, ‘seismic’, ‘Purples’, ‘ocean’, ‘PuOr’, ‘PuBuGn’, ‘nipy_spectral’, ‘afmhot’ --kmeans Number of clusters to compute. When this option is set, the matrix is split into clusters using the k-means algorithm. Only works for data that is not grouped, otherwise only the first group will be clustered. If more specific clustering methods are required, then save the underlying matrix and run the clustering using other software. The plotting of the clustering may fail with an error if a cluster has very few members compared to the total number or regions.

https://github.com/datirium/workflows.git

Path: workflows/deeptools.cwl

Branch/Commit ID: fa4f172486288a1a9d23864f1d6962d85a453e16

https://github.com/ncbi/pgap.git

Path: task_types/tt_kmer_seq_entry_extract_wnode.cwl

Branch/Commit ID: a2d6cd4c53bf3501f6bd79edebb7ca30bba8456f

https://github.com/kids-first/kf-alignment-workflow.git

Path: workflows/kf_alignment_fq_input_wf.cwl

Branch/Commit ID: 55315b6abb488f1f25fe725407814e8d4c23ba81