Differential Region Analysis with DESeq2

Overview

This skill performs differential region analysis between experimental conditions using DESeq2 in a count-based framework. Main steps include:

• Initialize the project directory.
• Refer to the Inputs & Outputs section to check inputs and build the output architecture. All the output file should located in

  ${proj_dir}

  in Step 0.
• Always prompt user if required files are missing.
• Always prompt user for the threshold of

  qvalues

  and

  log2foldchange

  to define significant regions.
• Merge peaks across replicates or samples to build a consensus peak set.
• Generate read count matrix over peaks using featureCounts or bedtools.
• Prepare sample metadata file describing conditions and replicates.
• Perform differential analysis using DESeq2.
• Visualize and interpret results (PCA, volcano plot).
• Output significantly up and down accessible regions.

When to use this skill

Use the differential-region-analysis pipeline when your goal is to identify genomic regions with condition-dependent changes in signal intensity, provided the signal can be represented as raw read counts per region.

Recommended scenarios include:

• Comparing treated vs. control samples to identify regulatory regions responsive to a drug, signaling molecule, or environmental change.
• Investigating cell differentiation or developmental trajectories to reveal dynamic chromatin remodeling.
• Analyzing disease vs. normal tissues to pinpoint dysregulated enhancer or promoter accessibility.
• Integrating with RNA-seq or ChIP-seq data to connect chromatin accessibility with transcriptional or epigenetic regulation.

The pipeline performs best with datasets containing biological replicates (≥2 per condition) and moderate to high sequencing depth (~20–50 million reads per sample).

Inputs & Outputs

Inputs (choose one)

• If starting from BAM files and BED peak files → Generate consensus peaks and count matrix.
• If starting from existing count matrix → Go directly to DESeq2 analysis.
• If multiple conditions or batches → Include batch/condition in design

Outputs

${sample}_DAR_analysis/ # or ${tf}_${sample}_DB_analysis in differential TF binding detection task
    tables/
      all_peaks.bed
      consensus_peaks.bed # Unified peak set
      atac_counts.txt # Count matrix of reads per peak
      samples.csv # Sample metadata
    DARs/
      DAR_results.csv # DESeq2 results (log2FC, p-values)
      DAR_sig.bed # Significantly diffential accessible regions
      DAR_up.bed
      DAR_down.bed  
    plots/ # visualization outputs
      PCA.pdf
      Volcano.pdf
    logs/ # analysis logs 
    temp/ # other temp files

Decision Tree

Step 0: Initialize Project

1. Make director for this project:

Call:

• mcp__project-init-tools__project_init

with:

• sample

  : sample name (e.g. c1_vs_c2)

• task

  : DAR_analysis

The tool will:

• Create

  ${sample}_DAR_analysis

  (or

  ${tf}_${sample}_DB_analysis

  ) directory.
• Return the full path of the

  ${sample}_DAR_analysis

  (or

  ${tf}_${sample}_DB_analysis

  ) directory, which will be used as

  ${proj_dir}

  .

Step 1: Generate Consensus Peaks

Combine peaks from replicates to define a shared feature space. Call:

• mcp__pydeseq2-tools__generate_consensus_peaks with:

• bed_files

  : List of paths to peak BED files from replicates.

• output_bed

  : Output path for the merged consensus BED file.

• output_saf

  : Output path for the SAF file (needed for featureCounts)

Output:

consensus_peaks.bed

consensus_peaks.saf

Step 2: Generate Count Matrix

Call:

• mcp__pydeseq2-tools__count_reads_featurecounts

with:

• saf_file

  : SAF file output from Step 1.

• bam_files

  : List of paths to BAM files.

• output_counts

  : Path to output count matrix.

• is_paired_end

  : Whether the BAM file is pair end or not.

• threads

Output:

atac_counts.txt

Step 3: Prepare Metadata

Prepare

samples.csv

sample,condition,replicate
sample1.bam,c1,1
sample2.bam,c1,2
sample3.bam,c2,1
sample4.bam,c2,2

Step 4: Differential Accessibility with pyDESeq2

Call:

• mcp__pydeseq2-tools__run_pydeseq2_analysis

with:

• counts_file: Path to featureCounts from Step 2.
• metadata_file: Path to metadata CSV from Step 3.
• design_factors: Design formula columns (e.g. 'condition' or 'batch,condition').
• contrast_column: Column name for contrast (e.g. 'condition').
• contrast_control: Control group name (e.g. 'Control').
• contrast_treatment: Treatment group name (e.g. 'Treated').
• output_csv: Output path for results CSV.

Output:

DAR_results.csv

${tf}_DB_results.csv

Step 5: Visualization and QC

Call:

• mcp__pydeseq2-tools__visualize_results

with:

• results_csv

  : Path to DESeq2 results CSV.

• counts_file

  : Path to original counts file (for PCA).

• metadata_file

  : Path to metadata (for PCA grouping).

• output_dir

  : Directory to save plots.

• condition_col

  : (e.g."condition")

Step 6: Output significantly up and down accessible regions

Call:

• mcp__pydeseq2-tools__filter_and_export_bed

with:

• results_csv

  : Path to DESeq2 results CSV.

• output_prefix

  : Prefix for output BED files.

• padj_cutoff

  : Provided by user

• log2fc_cutoff

  : Provided by user

Output:

DAR_sig.bed

DAR_up.bed

DAR_down.bed

${tf}_DB_sig.bed

${tf}_DB_up.bed

${tf}_DB_down.bed

Advanced Usage

• Batch effects:

  design = ~ batch + condition

• Multi-group comparison:

  contrast=c("condition","A","B")

• Time series:

  DESeq(dds, test="LRT", reduced=~1)

• Filter low counts:

  dds[rowSums(counts(dds)) >= 20, ]

Notes & Troubleshooting

rowSums >= 20

fitType="local"

betaPrior=FALSE