Explore Workflows

# Differential gene expression analysis This differential gene expression (DGE) analysis takes as input samples from two experimental conditions that have been processed with an RNA-Seq workflow (see list of \"Upstream workflows\" below). DESeq estimates variance-mean dependence in count data from high-throughput sequencing assays, then tests for DGE based on a model which assumes a negative binomial distribution of gene expression (aligned read count per gene). ### Experimental Setup and Results Interpretation The workflow design uses as its fold change (FC) calculation: condition 1 (c1, e.g. treatment) over condition 2 (c2, e.g. control). In other words: `FC == (c1/c2)` Therefore: - if FC<1 the log2(FC) is <0 (negative), meaning expression in condition1<condition2 (gene is downregulated in c1) - if FC>1 the log2(FC) is >0 (positive), meaning expression in condition1>condition2 (gene is upregulated in c1) In other words, if you have input TREATMENT samples as condition 1, and CONTROL samples as condition 2, a positive L2FC for a gene indicates that expression of the gene in TREATMENT is greater (or upregulated) compared to CONTROL. Next, threshold the p-adjusted values with your FDR (false discovery rate) cutoff to determine if the change may be considered significant or not. It is important to note when DESeq1 or DESeq2 is used in our DGE analysis workflow. If a user inputs only a single sample per condition DESeq1 is used for calculating DGE. In this experimental setup, there are no repeated measurements per gene per condition, therefore biological variability in each condition cannot be captured so the output p-values are assumed to be purely \"technical\". On the other hand, if >1 sample(s) are input per condition DESeq2 is used. In this case, biological variability per gene within each condition is available to be incorporated into the model, and resulting p-values are assumed to be \"biological\". Additionally, DESeq2 fold change is \"shrunk\" to account for sample variability, and as Michael Love (DESeq maintainer) puts it, \"it looks at the largest fold changes that are not due to low counts and uses these to inform a prior distribution. So the large fold changes from genes with lots of statistical information are not shrunk, while the imprecise fold changes are shrunk. This allows you to compare all estimated LFC across experiments, for example, which is not really feasible without the use of a prior\". In either case, the null hypothesis (H0) tested is that there are no significantly differentially expressed genes between conditions, therefore a smaller p-value indicates a lower probability of the H0 occurring by random chance and therefore, below a certain threshold (traditionally <0.05), H0 should be rejected. Additionally, due to the many thousands of independent hypotheses being tested (each gene representing an independent test), the p-values attained by the Wald test are adjusted using the Benjamini and Hochberg method by default. These \"padj\" values should be used for determination of significance (a reasonable value here would be <0.10, i.e. below a 10% FDR). Further Analysis: Output from the DESeq workflow may be used as input to the GSEA (Gene Set Enrichment Analysis) workflow for identifying enriched marker gene sets between conditions. ### DESeq1 High-throughput sequencing assays such as RNA-Seq, ChIP-Seq or barcode counting provide quantitative readouts in the form of count data. To infer differential signal in such data correctly and with good statistical power, estimation of data variability throughout the dynamic range and a suitable error model are required. Simon Anders and Wolfgang Huber propose a method based on the negative binomial distribution, with variance and mean linked by local regression and present an implementation, [DESeq](http://www.bioconductor.org/packages/3.8/bioc/html/DESeq.html), as an R/Bioconductor package. ### DESeq2 In comparative high-throughput sequencing assays, a fundamental task is the analysis of count data, such as read counts per gene in RNA-seq, for evidence of systematic changes across experimental conditions. Small replicate numbers, discreteness, large dynamic range and the presence of outliers require a suitable statistical approach. [DESeq2](http://www.bioconductor.org/packages/release/bioc/html/DESeq2.html), a method for differential analysis of count data, using shrinkage estimation for dispersions and fold changes to improve stability and interpretability of estimates. This enables a more quantitative analysis focused on the strength rather than the mere presence of differential expression. ### __References__ - Anders S, Huber W (2010). “Differential expression analysis for sequence count data.” Genome Biology, 11, R106. doi: 10.1186/gb-2010-11-10-r106, http://genomebiology.com/2010/11/10/R106/. - Love MI, Huber W, Anders S (2014). “Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.” Genome Biology, 15, 550. doi: 10.1186/s13059-014-0550-8.

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

Path: workflows/deseq.cwl

Branch/Commit ID: bf80c9339d81a78aefb8de661bff998ed86e836e

Wrapped in a workflow SoupX tool for easy access to Cell Ranger pipeline compressed outputs.

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

Path: tools/soupx-subworkflow.cwl

Branch/Commit ID: ebbf23764ede324cabc064bd50647c1f643726fa

Single-cell ATAC-Seq Dimensionality Reduction Analysis Integrates multiple single-cell ATAC-Seq datasets, reduces dimensionality using LSI.

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

Path: workflows/sc-atac-reduce.cwl

Branch/Commit ID: 22880e0f41d0420a17d643e8a6e8ee18165bbfbf

Sequence reads are first cleaned from adapters and transformed into fully bisulfite-converted forward (C->T) and reverse read (G->A conversion of the forward strand) versions, before they are aligned to similarly converted versions of the genome (also C->T and G->A converted). Sequence reads that produce a unique best alignment from the four alignment processes against the bisulfite genomes (which are running in parallel) are then compared to the normal genomic sequence and the methylation state of all cytosine positions in the read is inferred. A read is considered to align uniquely if an alignment has a unique best alignment score (as reported by the AS:i field). If a read produces several alignments with the same number of mismatches or with the same alignment score (AS:i field), a read (or a read-pair) is discarded altogether. On the next step we extract the methylation call for every single C analysed. The position of every single C will be written out to a new output file, depending on its context (CpG, CHG or CHH), whereby methylated Cs will be labelled as forward reads (+), non-methylated Cs as reverse reads (-). The output of the methylation extractor is then transformed into a bedGraph and coverage file. The bedGraph counts output is then used to generate a genome-wide cytosine report which reports the number on every single CpG (optionally every single cytosine) in the genome, irrespective of whether it was covered by any reads or not. As this type of report is informative for cytosines on both strands the output may be fairly large (~46mn CpG positions or >1.2bn total cytosine positions in the human genome).

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

Path: workflows/bismark-methylation-se.cwl

Branch/Commit ID: 87f213456b3f966b773d396cce1fe5a272dad858

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

Path: task_types/tt_cache_asnb_entries.cwl

Branch/Commit ID: 16e3915d2a357e2a861b30911c832e5ddc0c1784

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: 675a3ff982091faf304931e9261aacdbabf51702

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

Path: cwltool/schemas/v1.0/v1.0/step-valuefrom3-wf.cwl

Branch/Commit ID: 7c7615c44b80f8e76e659433f8c7875603ae0b25

DESeq2 Multi-factor Analysis Runs DeSeq2 multi-factor analysis with manual control over major parameters

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

Path: workflows/deseq-multi-factor.cwl

Branch/Commit ID: 7030da528559c7106d156284e50ff0ecedab0c4e

Sequence reads are first cleaned from adapters and transformed into fully bisulfite-converted forward (C->T) and reverse read (G->A conversion of the forward strand) versions, before they are aligned to similarly converted versions of the genome (also C->T and G->A converted). Sequence reads that produce a unique best alignment from the four alignment processes against the bisulfite genomes (which are running in parallel) are then compared to the normal genomic sequence and the methylation state of all cytosine positions in the read is inferred. A read is considered to align uniquely if an alignment has a unique best alignment score (as reported by the AS:i field). If a read produces several alignments with the same number of mismatches or with the same alignment score (AS:i field), a read (or a read-pair) is discarded altogether. On the next step we extract the methylation call for every single C analysed. The position of every single C will be written out to a new output file, depending on its context (CpG, CHG or CHH), whereby methylated Cs will be labelled as forward reads (+), non-methylated Cs as reverse reads (-). The output of the methylation extractor is then transformed into a bedGraph and coverage file. The bedGraph counts output is then used to generate a genome-wide cytosine report which reports the number on every single CpG (optionally every single cytosine) in the genome, irrespective of whether it was covered by any reads or not. As this type of report is informative for cytosines on both strands the output may be fairly large (~46mn CpG positions or >1.2bn total cytosine positions in the human genome).

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

Path: workflows/bismark-methylation-se.cwl

Branch/Commit ID: 9ee330737f4603e4e959ffe786fbb2046db70a00

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

Path: v1.0/v1.0/search.cwl

Branch/Commit ID: 9a23706ec061c5d2c02ff60238d218aadf0b5db9

Packed ID: main