# QC and Trimming for CUT&RUN Data

CUT&RUN produces paired-end reads with a characteristic nucleosomal fragment
size distribution. Reads are typically shorter than ChIP-seq due to the
MNase cleavage mechanism.

## Pre-Trimming QC with FastQC

```bash
fastqc --threads 4 --outdir fastqc_raw/ sample_R1.fastq.gz sample_R2.fastq.gz
```

Key checks:
- Per-base quality (expect Phred >28)
- Adapter content (Nextera or Illumina universal adapter)
- Insert size (CUT&RUN fragments are often <150 bp, causing adapter read-through)
- Sequence duplication (moderate levels expected with low input)

**Note**: CUT&RUN libraries from low cell input often show higher duplication
than ChIP-seq. This is expected and not necessarily a quality issue.

## Adapter Trimming with Trim Galore

CUT&RUN fragments are often shorter than read length, causing adapter
read-through. Aggressive adapter trimming is important.

```bash
trim_galore \
    --paired \
    --quality 20 \
    --phred33 \
    --length 20 \
    --cores 4 \
    --fastqc \
    --nextera \
    --output_dir trim_galore/ \
    sample_R1.fastq.gz \
    sample_R2.fastq.gz
```

### Parameter Rationale

| Parameter | Value | Reason |
|-----------|-------|--------|
| `--quality 20` | Phred 20 | Standard quality cutoff |
| `--length 20` | 20 bp | Keep short sub-nucleosomal fragments |
| `--nextera` | Flag | CUT&RUN often uses Nextera adapters (check protocol) |
| `--cores 4` | 4 | Parallel processing |

### Adapter Type Selection

Check which adapter was used in library preparation:
- **Nextera**: Most CUT&Tag protocols (use `--nextera` flag)
- **Illumina TruSeq**: Some CUT&RUN protocols (default Trim Galore detection)
- If unsure, let Trim Galore auto-detect (omit `--nextera`)

## Post-Trimming Verification

After trimming, verify:
- >90% of reads pass quality filter
- Adapter contamination removed (often 20-50% for CUT&RUN)
- Read length distribution: 20-150 bp (many short reads are expected)

High adapter contamination rate (>30%) is NORMAL for CUT&RUN because
many fragments are shorter than the read length.

## Fragment Size Check (Post-Alignment)

After alignment, verify the fragment size distribution:

```bash
samtools view -f 2 -F 1804 sample.bam | \
    awk '{if($9 > 0) print $9}' | \
    sort -n | uniq -c | \
    awk '{print $2, $1}' > fragment_sizes.txt
```

Expected CUT&RUN fragment distribution:
- **TF targets**: Peak at <120 bp (sub-nucleosomal)
- **Histone marks**: Strong peak at ~150 bp (mononucleosomal)
- **Both**: Should show nucleosomal ladder pattern
- Broad smear with no peaks suggests protocol failure

## CUT&Tag vs CUT&RUN Trimming

CUT&Tag uses Tn5 transposase which adds 19 bp Mosaic End (ME) sequences.
For CUT&Tag data specifically:

```bash
trim_galore \
    --paired \
    --quality 20 \
    --length 20 \
    --cores 4 \
    --nextera \
    --fastqc \
    sample_R1.fastq.gz \
    sample_R2.fastq.gz
```

The `--nextera` flag is especially important for CUT&Tag since Tn5 inserts
Nextera-compatible adapters.

# Bowtie2 Alignment for CUT&RUN (Genome + Spike-in)

CUT&RUN uses Bowtie2 for alignment with specific settings optimized for the
short, paired-end fragments produced by MNase or Tn5 cleavage. Two separate
alignments are required: one to the target genome and one to the E. coli
spike-in genome.

## Genome Index Preparation

```bash
# Target genome (one-time)
bowtie2-build --threads 8 genome.fa genome_index

# E. coli spike-in genome (one-time)
# Use E. coli K12 MG1655 (GenBank: U00096.3)
bowtie2-build --threads 4 ecoli_K12.fa ecoli_index
```

## Target Genome Alignment

```bash
bowtie2 \
    --very-sensitive \
    --no-mixed \
    --no-discordant \
    --dovetail \
    --phred33 \
    -I 10 -X 700 \
    --threads 8 \
    -x genome_index \
    -1 sample_R1_val_1.fq.gz \
    -2 sample_R2_val_2.fq.gz \
    2> sample_bowtie2.log \
    | samtools view -@ 4 -bS - \
    | samtools sort -@ 4 -o sample_sorted.bam

samtools index sample_sorted.bam
```

### Bowtie2 Parameters for CUT&RUN

| Parameter | Value | Reason |
|-----------|-------|--------|
| `--very-sensitive` | Preset | Maximum sensitivity for short fragments |
| `--no-mixed` | Flag | Both mates must align |
| `--no-discordant` | Flag | Mates must be properly paired |
| `--dovetail` | Flag | Allow dovetail alignments (overlapping PE reads) |
| `-I 10` | 10 bp | Minimum insert size (very short fragments exist) |
| `-X 700` | 700 bp | Maximum insert size |

**Critical**: The `--dovetail` flag is essential because CUT&RUN fragments
are often shorter than read length, causing R1 and R2 to extend past each
other (dovetailing).

## Spike-in Alignment

Align reads that did NOT map to the target genome to the E. coli spike-in:

```bash
# Extract unmapped reads from genome alignment
samtools view -b -f 12 -F 256 sample_sorted.bam \
    | samtools sort -@ 4 -n -o unmapped_sorted.bam

bedtools bamtofastq -i unmapped_sorted.bam \
    -fq unmapped_R1.fq -fq2 unmapped_R2.fq

# Align to E. coli
bowtie2 \
    --very-sensitive \
    --no-mixed \
    --no-discordant \
    --dovetail \
    --phred33 \
    -I 10 -X 700 \
    --threads 4 \
    -x ecoli_index \
    -1 unmapped_R1.fq \
    -2 unmapped_R2.fq \
    2> sample_spikein.log \
    | samtools view -@ 2 -bS -q 10 -F 1804 -f 2 - \
    | samtools sort -@ 2 -o sample_spikein.bam

samtools index sample_spikein.bam
```

### Alternative: Direct Spike-in Alignment

Some workflows align ALL reads to both genomes simultaneously using a
concatenated index. This is simpler but less precise:

```bash
# Concatenate genomes (one-time)
cat genome.fa ecoli_K12.fa > combined.fa
bowtie2-build --threads 8 combined.fa combined_index

# Align to combined genome
bowtie2 --very-sensitive --no-mixed --no-discordant --dovetail \
    -I 10 -X 700 --threads 8 \
    -x combined_index \
    -1 R1.fq.gz -2 R2.fq.gz \
    | samtools view -bS - > combined.bam

# Separate by genome
samtools view -b combined.bam chr1 chr2 ... chrX chrY > genome.bam
samtools view -b combined.bam ecoli_chr > spikein.bam
```

## Spike-in Read Counts

```bash
# Count spike-in reads
spikein_count=$(samtools view -c -F 1804 -f 2 sample_spikein.bam)
genome_count=$(samtools view -c -F 1804 -f 2 sample_sorted.bam)
echo "Spike-in reads: ${spikein_count}"
echo "Genome reads: ${genome_count}"
echo "Spike-in fraction: $(echo "scale=4; $spikein_count / ($spikein_count + $genome_count)" | bc)"
```

Expected spike-in fraction:
- **1-10%**: Optimal range for normalization
- **<0.1%**: Too few spike-in reads for reliable normalization
- **>30%**: Excessive spike-in; may indicate poor target enrichment

## Alignment QC

```bash
# Parse Bowtie2 log for mapping rate
grep "overall alignment rate" sample_bowtie2.log
```

Expected mapping rates:
- **Target genome**: >80% for standard samples
- **Spike-in**: spike-in reads should be 1-10% of all reads (spike-in / (spike-in + genome), as computed above)

Low genome mapping rate may indicate:
- Contamination (check FastQ Screen)
- Wrong genome build
- Very poor library quality

# Filtering, Deduplication, and Spike-in Normalization

CUT&RUN data processing requires standard quality filtering, duplicate removal,
blacklist filtering, AND spike-in normalization for quantitative analysis.

## Quality Filtering

```bash
samtools view -b -h \
    -q 10 \
    -F 1804 \
    -f 2 \
    sample_sorted.bam \
    | samtools sort -@ 4 -o sample_filtered.bam
```

### Filter Parameters

| Flag | Meaning |
|------|---------|
| `-q 10` | MAPQ >= 10 (CUT&RUN uses lower threshold than ChIP-seq) |
| `-F 4` | Remove unmapped |
| `-F 256` | Remove secondary |
| `-F 512` | Remove QC-fail |
| `-F 1024` | Remove duplicates -- a no-op here; duplicates are removed by Picard in the next step |
| `-f 2` | Keep properly paired only |

**Note**: MAPQ 10 instead of 30 for CUT&RUN. The lower threshold retains
more signal because CUT&RUN targets can be in repetitive regions.

## Duplicate Marking

```bash
picard MarkDuplicates \
    INPUT=sample_filtered.bam \
    OUTPUT=sample_dedup.bam \
    METRICS_FILE=sample_dup_metrics.txt \
    REMOVE_DUPLICATES=true \
    VALIDATION_STRINGENCY=LENIENT \
    ASSUME_SORTED=true

samtools index sample_dedup.bam
```

CUT&RUN from low cell numbers may have higher duplication. Accept up to 40%.

## Blacklist Filtering

The workflow filters the BAM against the single file given as `--blacklist`
and stops there: peak files are never filtered, and there is no separate
suspect-list parameter.

```bash
# What the workflow runs, with --blacklist as -b
bedtools intersect \
    -a sample_dedup.bam \
    -b hg38-blacklist.v2.bed \
    -v \
    > sample_final.bam

samtools index sample_final.bam
```

To also exclude the CUT&RUN-specific suspect list (Nordin 2023) -- ~400
regions enriched in CUT&RUN controls that produce false positive peaks,
independent of the ENCODE blacklist -- merge the two files once and pass the
result as `--blacklist`:

```bash
cat hg38-blacklist.v2.bed CUTandRUN.suspectlist.hg38.bed \
    | sort -k1,1 -k2,2n | bedtools merge > combined_blacklist.bed
```

## Spike-in Normalization

### Calculate Scale Factors

The workflow writes one count file per sample
(`spikein/{sample}.spikein_counts.txt`, columns sample and count) and then one
combined `spikein/scale_factors.txt` for the whole run:

```bash
# Per-sample counts, concatenated
cat *.spikein_counts.txt > all_counts.txt

# Scale every sample to the smallest non-zero count (factor = min / count).
# A sample with no spike-in reads is left unscaled (factor 1).
min_count=$(awk -F'\t' '$2 > 0 {print $2}' all_counts.txt | sort -n | head -1)

awk -F'\t' -v min="${min_count:-0}" 'BEGIN {OFS="\t"} {
    factor = ($2 > 0 && min > 0) ? min / $2 : 1
    print $1, $2, factor
}' all_counts.txt > scale_factors.txt
```

### Apply Spike-in Scaling to Signal

The workflow scales the bigWig only, with deepTools:

```bash
# Read scale factor for one sample
scale=$(awk -F'\t' -v s="sample1" '$1==s {print $3}' scale_factors.txt)

bamCoverage \
    --bam sample_final.bam \
    --outFileName sample_normalized.bw \
    --scaleFactor ${scale} \
    --binSize 10 \
    --normalizeUsing None \
    --extendReads \
    --numberOfProcessors 4
```

Without a spike-in index the workflow drops `--scaleFactor` and uses
`--normalizeUsing RPKM` instead.

### Alternative: bedGraph route (not used by the workflow)

```bash
bedtools genomecov \
    -ibam sample_final.bam \
    -bg \
    -pc \
    -scale ${scale} \
    -g hg38.chrom.sizes \
    | sort -k1,1 -k2,2n > sample_normalized.bedGraph

bedGraphToBigWig sample_normalized.bedGraph hg38.chrom.sizes sample_normalized.bw
```

## Generate Fragment BED File

SEACR takes a fragment bedGraph, built from a fragment BED. `bedtools
bamtobed -bedpe` needs mates on adjacent lines, so name-sort the BAM first --
on a coordinate-sorted BAM it warns per read that the mate does not occur next
to it and emits a near-empty BED:

```bash
# Name-sort first
samtools sort -n -@ 2 -o namesorted.bam sample_final.bam

# Keep properly paired fragments on one chromosome and shorter than 1 kb
bedtools bamtobed -bedpe -i namesorted.bam \
    | awk 'BEGIN {OFS="\t"} $1 == $4 && $6 - $2 < 1000 {print $1, $2, $6}' \
    | sort -k1,1 -k2,2n -k3,3n \
    > sample.fragments.bed

# Fragment bedGraph for SEACR -- deliberately unscaled
bedtools genomecov \
    -i sample.fragments.bed \
    -g hg38.chrom.sizes \
    -bg \
    > sample_fragments.bedGraph
```

The workflow does not apply the spike-in factor here: SEACR sees unscaled
fragment coverage, and only the bigWig is calibrated. The BED is published as
`signal/{sample}.fragments.bed`; the bedGraph is an intermediate.

## QC Statistics

```bash
# Final read count
echo "Final reads: $(samtools view -c sample_final.bam)"

# Flagstat (the workflow publishes this as alignment/{sample}.flagstat.txt)
samtools flagstat sample_final.bam > sample_flagstat.txt

# Spike-in fraction; the spike-in BAM is an intermediate, so read the
# published count file instead
genome=$(samtools view -c -F 1804 -f 2 sample_final.bam)
spikein=$(cut -f2 spikein/sample.spikein_counts.txt)
echo "Spike-in fraction: $(echo "scale=4; $spikein / ($genome + $spikein)" | bc)"
```

# SEACR Peak Calling for CUT&RUN

SEACR (Sparse Enrichment Analysis for CUT&RUN) is specifically designed for
the sparse, low-background signal profile of CUT&RUN data. It outperforms
MACS2 on CUT&RUN data because it does not assume a Poisson background model.

## SEACR with IgG Control

When an IgG control is available:

```bash
SEACR_1.3.sh \
    sample_fragments.bedGraph \
    control_IgG_fragments.bedGraph \
    norm \
    stringent \
    sample_seacr
```

### Parameters

| Argument | Value | Description |
|----------|-------|-------------|
| Arg 1 | sample bedGraph | Treatment signal (fragment coverage) |
| Arg 2 | control bedGraph | IgG/no-antibody control signal |
| Arg 3 | `norm` | Normalize control to treatment depth |
| Arg 4 | `stringent` | Stringent mode (conservative peaks) |
| Arg 5 | output prefix | Output file prefix |

## SEACR without Control (Top %)

When no control is available, SEACR uses a numeric threshold:

```bash
SEACR_1.3.sh \
    sample_fragments.bedGraph \
    0.01 \
    non \
    stringent \
    sample_seacr_noctrl
```

`non` is required here: the SEACR v1.3 README states that a numeric threshold
must be paired with `non`. The v1.3 script does not error on `norm` -- it
silently skips normalization -- so the mistake is invisible. The workflow
always passes `non` when there is no `--control`.

The `0.01` value means the top 1% of signal is used as the enrichment threshold.
Adjust based on expected peak count:
- `0.01` (1%): Conservative, fewer peaks
- `0.05` (5%): Moderate
- `0.10` (10%): Permissive, more peaks

## SEACR Stringent vs Relaxed

Run both modes to compare:

```bash
# Stringent: Only peaks that pass both global and local enrichment
SEACR_1.3.sh sample.bedGraph control.bedGraph norm stringent sample_stringent

# Relaxed: Peaks passing global enrichment threshold only
SEACR_1.3.sh sample.bedGraph control.bedGraph norm relaxed sample_relaxed
```

| Mode | Description | Typical Peak Count |
|------|-------------|-------------------|
| Stringent | Global AND local enrichment | 5,000-20,000 |
| Relaxed | Global enrichment only | 10,000-50,000 |

## SEACR Output Format

SEACR v1.3 outputs a 6-column BED-like file:

```
chr1  1000  2000  500.5  8.2  chr1:1450-1480
```

Columns:
1. Chromosome
2. Peak start
3. Peak end
4. Total signal in the peak (AUC)
5. Max signal in the peak
6. Region of maximum signal, as `chr:start-end`

The workflow names these files `{sample}.seacr.stringent.bed` and, with
`--seacr_mode relaxed` or `both`, `{sample}.seacr.relaxed.bed`, published to
`peaks/`.

## Alternative: MACS2 Peak Calling

MACS2 can be used as an alternative or validation:

```bash
macs2 callpeak \
    -t sample_final.bam \
    -c control_IgG.bam \
    -f BAMPE \
    -g hs \
    -n sample_macs2 \
    --nomodel \
    --keep-dup all \
    -q 0.05 \
    --outdir macs2_peaks/
```

### MACS2 Parameters for CUT&RUN

| Parameter | Value | Reason |
|-----------|-------|--------|
| `-f BAMPE` | Paired-end | Use fragment information |
| `--nomodel` | Skip model | CUT&RUN fragments don't follow ChIP model |
| `--keep-dup all` | Keep all | Duplicates already removed |
| `-q 0.05` | FDR 0.05 | Standard threshold |

**Caution**: MACS2 may overcall peaks on CUT&RUN data due to the low
background. Compare with SEACR results and use the intersection for
high-confidence peaks.

## Blacklist + Suspect List Filtering (manual, not run by this workflow)

The workflow applies `--blacklist` to the BAM only and never filters peak
files. To filter the published peaks against both the blacklist and the
CUT&RUN suspect list, run this yourself:

```bash
bedtools intersect \
    -a results/peaks/sample.seacr.stringent.bed \
    -b hg38-blacklist.v2.bed CUTandRUN.suspectlist.hg38.bed \
    -v \
    > sample_peaks_filtered.bed
```

## Peak Overlap Between Callers

When using both SEACR and MACS2, assess concordance:

```bash
# Convert SEACR to 3-column BED
cut -f1-3 results/peaks/sample.seacr.stringent.bed > seacr_peaks.bed
cut -f1-3 results/peaks/sample.macs2_peaks.narrowPeak > macs2_peaks.bed

# Overlap
bedtools intersect -a seacr_peaks.bed -b macs2_peaks.bed -u | wc -l
total_seacr=$(wc -l < seacr_peaks.bed)
total_macs2=$(wc -l < macs2_peaks.bed)
echo "SEACR peaks: $total_seacr"
echo "MACS2 peaks: $total_macs2"
```

Typical overlap: 60-80% of SEACR peaks overlap MACS2 peaks.
High-confidence set: intersection of both callers.

## FRiP Calculation

The FRIP process runs this for every peak set called for a sample and writes
the results to `qc/{sample}.frip_mqc.tsv`:

```bash
total_reads=$(samtools view -c results/alignment/sample.filtered.bam)
reads_in_peaks=$(bedtools intersect \
    -u \
    -a results/alignment/sample.filtered.bam \
    -b results/peaks/sample.seacr.stringent.bed \
    | samtools view -c -)
awk -v a="$reads_in_peaks" -v b="$total_reads" 'BEGIN { printf "FRiP: %.4f\n", a / b }'
```

Both counts are alignments, with no flag filter, so mates count separately.
Judge the result against the QC table in `SKILL.md`; see
`05-qc-metrics.md` for the published file's columns.

# CUT&RUN QC Metrics

Quality assessment for CUT&RUN data includes standard alignment metrics,
CUT&RUN-specific fragment analysis, spike-in validation, and peak quality.

Pass/warn/fail thresholds live in one place: the QC table in `SKILL.md`. The
commands below show how to derive each value, including the ones the workflow
does not compute (spike-in fraction, peak statistics).

## Fragment Size Distribution

The fragment size distribution is the most informative CUT&RUN QC metric:

```bash
# Extract fragment sizes from properly paired reads
samtools view -f 2 -F 1804 sample_final.bam | \
    awk '{if($9 > 0 && $9 < 1000) print $9}' | \
    sort -n | uniq -c | \
    awk '{print $2, $1}' > fragment_sizes.txt
```

The workflow runs exactly this and publishes
`qc/{sample}.fragment_sizes.txt`.

### Expected Patterns by Target

| Target Type | Fragment Pattern | Example |
|-------------|-----------------|---------|
| TF (CTCF, etc.) | Peak <120 bp, some at 150 bp | Sharp sub-nucleosomal |
| Active histone (H3K4me3) | Strong 150 bp peak | Mononucleosomal |
| Repressive histone (H3K27me3) | 150 bp + 300 bp | Mono + dinucleosomal |
| IgG control | Flat distribution | No enrichment pattern |

### Red Flags in Fragment Distribution

- No nucleosomal periodicity: Protocol may have failed
- Only large fragments (>300 bp): Over-digestion or poor tagmentation
- Spike at exact read length: Adapter trimming incomplete
- Identical to IgG: No target enrichment

## Spike-in QC

### Spike-in Fraction

The spike-in BAM is an intermediate; the published per-sample count file is
`spikein/{sample}.spikein_counts.txt` (columns sample and count):

```bash
genome=$(samtools view -c -F 1804 -f 2 results/alignment/sample.filtered.bam)
spikein=$(cut -f2 results/spikein/sample.spikein_counts.txt)
fraction=$(echo "scale=4; $spikein / ($genome + $spikein)" | bc)
echo "Spike-in fraction: $fraction"
```

Interpretation of the fraction (thresholds: see the QC table in `SKILL.md`):
too little spike-in makes the normalization imprecise, too much indicates
poor target enrichment.

### Spike-in Consistency Across Samples

For reliable normalization, spike-in counts should vary across samples
(reflecting different amounts of target material), but not be zero. The
workflow already collects them:

```bash
sort -k2 -n results/spikein/scale_factors.txt
```

Columns: sample, spike-in count, scale factor.

## Alignment Statistics

The workflow publishes `alignment/{sample}.flagstat.txt` and
`alignment/{sample}.dup_metrics.txt`, and the Bowtie2 logs reach
`multiqc/multiqc_report.html`. To recompute:

```bash
samtools flagstat results/alignment/sample.filtered.bam > flagstat.txt
```

Mapping rate, properly paired fraction and duplication rate are judged against
the QC table in `SKILL.md`.

## FRiP (Fraction of Reads in Peaks)

The FRIP process computes this and publishes `qc/{sample}.frip_mqc.tsv`, one
row per peak set called for the sample (SEACR stringent and/or relaxed, and/or
MACS2). Columns: `Peak set` (the peak file), `FRiP`, `reads_in_peaks`,
`total_reads`. MultiQC shows the same table as "Fraction of reads in peaks".

This is what it runs for each peak set:

```bash
total=$(samtools view -c results/alignment/sample.filtered.bam)
in_peaks=$(bedtools intersect \
    -u \
    -a results/alignment/sample.filtered.bam \
    -b results/peaks/sample.seacr.stringent.bed \
    | samtools view -c -)
awk -v a="$in_peaks" -v b="$total" 'BEGIN { printf "FRiP: %.4f\n", a / b }'
```

Both counts are alignments of the filtered BAM with no further flag filter, so
the two mates of a fragment count separately. Counting differently (for
example with `-F 1804 -f 2`) gives a number that does not match the published
one.

CUT&RUN typically has higher FRiP than ChIP-seq because of lower background.
Use the FRiP row of the QC table in `SKILL.md` for the thresholds.

## Peak Count and Size

```bash
# Peak statistics
total_peaks=$(wc -l < results/peaks/sample.seacr.stringent.bed)
echo "Total peaks: $total_peaks"

# Peak size distribution
awk '{print $3-$2}' results/peaks/sample.seacr.stringent.bed | \
    awk '{sum+=$1; n++; a[n]=$1} END {
        asort(a);
        print "Median peak size:", a[int(n/2)];
        print "Mean peak size:", sum/n;
        print "Min:", a[1];
        print "Max:", a[n]
    }'
```

## MultiQC Aggregation

```bash
multiqc \
    --title "CUT&RUN Pipeline QC" \
    --filename multiqc_report \
    --outdir multiqc/ \
    fastqc/ trim_galore/ alignment/ qc/
```

## Summary QC Table

Generate a per-sample summary:

```bash
echo -e "Sample\tReads\tMap_Rate\tDedup_Rate\tSpikein_Frac\tPeaks\tFRiP\tFrag_Peak"
echo -e "${SAMPLE}\t${TOTAL}\t${MAP_RATE}\t${DUP_RATE}\t${SPIKEIN}\t${PEAKS}\t${FRIP}\t${FRAG}"
```

# CUT&RUN Pipeline — Literature References

**Last updated:** 2026-03-07
**Purpose:** Reference catalog for the pipeline-cutandrun skill — papers defining CUT&RUN and CUT&Tag methods, SEACR peak calling, spike-in normalization, and the CUT&RUN suspect list for artifact filtering.

---

## CUT&RUN / CUT&Tag Method Development

---

### Skene & Henikoff 2017 — CUT&RUN: targeted nuclease strategy

- **Citation:** Skene PJ, Henikoff S. An efficient targeted nuclease strategy for high-resolution mapping of DNA binding sites. eLife, 6:e21856, 2017.
- **DOI:** [10.7554/eLife.21856](https://doi.org/10.7554/eLife.21856)
- **PMID:** 28079019 | **PMC:** PMC5310842
- **Citations:** ~1,500
- **Key findings:** Introduced Cleavage Under Targets and Release Using Nuclease (CUT&RUN), which tethers protein A-Micrococcal Nuclease (pA-MNase) to an antibody bound to a target protein in situ, then activates cleavage by calcium addition. Released chromatin fragments diffuse out of the nucleus and are collected from the supernatant, eliminating the need for chromatin solubilization or immunoprecipitation. Produces extremely low background compared to ChIP-seq because only targeted fragments are released, while bulk chromatin remains in the nucleus. Requires as few as 100 cells and generates high-quality profiles for both histone marks and transcription factors with minimal sequencing depth (~5M reads vs 20-40M for ChIP-seq). The E. coli DNA carried over from pA-MNase production serves as an internal spike-in for normalization.

---

### Kaya-Okur et al. 2019 — CUT&Tag: efficient epigenomic profiling

- **Citation:** Kaya-Okur HS, Wu SJ, Codomo CA, Pledger ES, Bryson TD, Henikoff JG, Ahmad K, Henikoff S. CUT&Tag for efficient epigenomic profiling of small samples and single cells. Nature Communications, 10(1):1930, 2019.
- **DOI:** [10.1038/s41467-019-09982-5](https://doi.org/10.1038/s41467-019-09982-5)
- **PMID:** 31036827 | **PMC:** PMC6488672
- **Citations:** ~1,200
- **Key findings:** Introduced Cleavage Under Targets and Tagmentation (CUT&Tag), which replaces pA-MNase with protein A-Tn5 transposase (pA-Tn5). After antibody tethering, Tn5 simultaneously fragments and tags target chromatin with sequencing adapters. Key advantages over CUT&RUN: streamlined single-tube workflow, compatible with single-cell applications, and no separate library preparation step needed. However, CUT&Tag has higher background than CUT&RUN due to Tn5 insertion preferences in accessible chromatin. Both methods are supported by this pipeline with the same alignment and peak calling workflow, differing only in spike-in source (E. coli from Tn5 vs MNase production).

---

### Skene & Henikoff 2018 — CUT&RUN protocol optimization

- **Citation:** Skene PJ, Henikoff JG, Henikoff S. Targeted in situ genome-wide profiling with high efficiency for low cell numbers. Nature Protocols, 13(5):1006-1019, 2018.
- **DOI:** [10.1038/nprot.2018.015](https://doi.org/10.1038/nprot.2018.015)
- **PMID:** 29651053
- **Citations:** ~500
- **Key findings:** Detailed step-by-step CUT&RUN protocol with optimization guidelines for different targets (histone marks vs TFs). Established the standard experimental parameters: ConA bead binding for cell immobilization, 0°C calcium-activated cleavage for 30 minutes, fragment release at 37°C, and library preparation with NEBNext reagents. Provided troubleshooting guidance for common issues including low yield, high background, and incomplete cleavage. This protocol forms the basis for most published CUT&RUN experiments and informs the QC expectations in this pipeline.

---

## Peak Calling

---

### Meers et al. 2019 — SEACR: peak calling for CUT&RUN

- **Citation:** Meers MP, Tenenbaum D, Henikoff S. Peak calling by Sparse Enrichment Analysis for CUT&RUN chromatin profiling. Epigenetics & Chromatin, 12(1):42, 2019.
- **DOI:** [10.1186/s13072-019-0287-4](https://doi.org/10.1186/s13072-019-0287-4)
- **PMID:** 31300027 | **PMC:** PMC6626385
- **Citations:** ~800
- **Key findings:** Introduced SEACR (Sparse Enrichment Analysis for CUT&RUN), a peak caller specifically designed for the sparse, low-background signal profile of CUT&RUN data. Unlike MACS2 (which models background as a Poisson distribution), SEACR uses the empirical distribution of signal in control (IgG) or target data to identify enriched regions without parametric assumptions. Offers two modes: stringent (peaks must exceed both the global threshold and a local enrichment test) and relaxed (global threshold only). When no IgG control is available, SEACR uses the top 1% of target signal as a threshold. Benchmarking showed SEACR produces fewer false positives than MACS2 on CUT&RUN data, particularly in regions with low but genuine enrichment.

---

## Quality Control

---

### Nordin et al. 2023 — CUT&RUN suspect list

- **Citation:** Nordin A, Zambanini G, Pagella P, Bhatt DK, Bjork P, Nilsson J, Mead P, Boyle AP. The CUT&RUN suspect list of problematic regions of the genome. Genome Biology, 24:185, 2023.
- **DOI:** [10.1186/s13059-023-02960-3](https://doi.org/10.1186/s13059-023-02960-3)
- **PMID:** 37580722 | **PMC:** PMC10424377
- **Citations:** ~50
- **Key findings:** Identified genomic regions that produce artifactual signal in CUT&RUN and CUT&Tag experiments, distinct from the standard ENCODE blacklist. These "suspect" regions include areas with high MNase/Tn5 accessibility, specific repeat families, and regions prone to antibody-independent cleavage. The suspect list should be applied in addition to the ENCODE blacklist v2 (Amemiya et al. 2019) when filtering CUT&RUN/CUT&Tag peaks. Without suspect list filtering, up to 20% of called peaks may be artifacts, particularly for targets with moderate enrichment. This workflow applies only the single BED given as `--blacklist`, to the BAM; pass a merged blacklist + suspect list there, or filter the peaks manually.

---

### Meers et al. 2019 — Spike-in normalization for CUT&RUN/CUT&Tag

- **Citation:** Meers MP, Bryson TD, Henikoff JG, Henikoff S. Improved CUT&RUN chromatin profiling tools. eLife, 8:e46314, 2019.
- **DOI:** [10.7554/eLife.46314](https://doi.org/10.7554/eLife.46314)
- **PMID:** 31232687 | **PMC:** PMC6632061
- **Citations:** ~400
- **Key findings:** Established the spike-in calibration framework for quantitative CUT&RUN. E. coli DNA carried over from pA-MNase production provides an internal standard: the ratio of spike-in reads between samples inversely correlates with target enrichment efficiency. Scale factors computed from spike-in counts enable quantitative comparison of signal intensity across samples and conditions. Demonstrated that spike-in normalization is essential for detecting global changes in histone modifications (e.g., drug treatments that globally increase or decrease a mark), which would be invisible with standard library-size normalization. This calibration approach is implemented in the pipeline's spike-in normalization step.

---

## Supplementary Tools (Non-CUT&RUN-Specific)

See pipeline-chipseq/references/literature.md for detailed descriptions of shared tools, and pipeline-atacseq/references/literature.md for Bowtie2.

---

### Langmead & Salzberg 2012 — Bowtie 2

- **DOI:** [10.1038/nmeth.1923](https://doi.org/10.1038/nmeth.1923) | **PMID:** 22388286 | **Citations:** ~47,300
- **CUT&RUN role:** Alignment of CUT&RUN/CUT&Tag reads to both the target genome and the E. coli spike-in genome. Bowtie2 is preferred over BWA-MEM for CUT&RUN because the short fragments typical of CUT&RUN (~150 bp mononucleosomal, <120 bp sub-nucleosomal) are well-suited to Bowtie2's alignment algorithm. Run with --very-sensitive and --dovetail flags.

---

### Zhang et al. 2008 — MACS2

- **DOI:** [10.1186/gb-2008-9-9-r137](https://doi.org/10.1186/gb-2008-9-9-r137) | **PMID:** 18798982 | **Citations:** ~7,000
- **CUT&RUN role:** Alternative peak caller for CUT&RUN data. MACS2 can be used instead of or alongside SEACR, but tends to overcall peaks due to CUT&RUN's low background. When using MACS2 on CUT&RUN data, use --nomodel and set appropriate --shift/--extsize parameters. Useful for comparison with ChIP-seq results processed with the same peak caller.

---

### Li et al. 2009 — SAMtools

- **DOI:** [10.1093/bioinformatics/btp352](https://doi.org/10.1093/bioinformatics/btp352) | **PMID:** 19505943 | **Citations:** ~53,700
- **CUT&RUN role:** BAM sorting, filtering, indexing, and alignment statistics.

---

### Broad Institute — Picard MarkDuplicates

- **URL:** [https://broadinstitute.github.io/picard/](https://broadinstitute.github.io/picard/)
- **CUT&RUN role:** PCR duplicate marking. CUT&RUN typically has lower duplication rates than ChIP-seq because it requires less input material and fewer PCR cycles. Duplication rates >20% suggest overamplification.

---

### Quinlan & Hall 2010 — BEDTools

- **DOI:** [10.1093/bioinformatics/btq033](https://doi.org/10.1093/bioinformatics/btq033) | **PMID:** 20110278 | **Citations:** ~12,000
- **CUT&RUN role:** Used by the workflow to filter the BAM against `--blacklist`, to build the fragment BED and bedGraph, and to count reads in peaks for FRiP (`bedtools intersect -u`); used manually for peak filtering, which the workflow does not run.

---

### Ramírez et al. 2016 — deepTools2

- **DOI:** [10.1093/nar/gkw257](https://doi.org/10.1093/nar/gkw257) | **PMID:** 27079975 | **Citations:** ~6,100
- **CUT&RUN role:** Generates spike-in normalized signal tracks (bigWig) using bamCoverage with the --scaleFactor parameter computed from spike-in read counts. Also used for fingerprint plots and heatmaps.

---

### Amemiya et al. 2019 — ENCODE Blacklist

- **DOI:** [10.1038/s41598-019-45839-z](https://doi.org/10.1038/s41598-019-45839-z) | **PMID:** 31249361 | **Citations:** ~1,372
- **CUT&RUN role:** The file normally passed as `--blacklist`, which the workflow applies to the BAM. Both filters are worth applying: the blacklist addresses general sequencing/alignment artifacts while the CUT&RUN suspect list (Nordin 2023) addresses CUT&RUN-specific enzyme cleavage artifacts. Merge the two into one BED and pass it as `--blacklist`, or filter the peaks manually.

---

### Ewels et al. 2016 — MultiQC

- **DOI:** [10.1093/bioinformatics/btw354](https://doi.org/10.1093/bioinformatics/btw354) | **PMID:** 27312411 | **Citations:** ~6,800
- **CUT&RUN role:** Aggregates QC metrics from FastQC, Bowtie2, Picard, SEACR, and spike-in statistics into a unified HTML report.

---

### Andrews 2010 — FastQC

- **URL:** [https://www.bioinformatics.babraham.ac.uk/projects/fastqc/](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)
- **CUT&RUN role:** Raw read quality assessment. CUT&RUN libraries typically show high-quality bases and characteristic nucleosomal fragment sizes.

---

### Martin 2011 — Cutadapt (basis for Trim Galore)

- **DOI:** [10.14806/ej.17.1.200](https://doi.org/10.14806/ej.17.1.200) | **Citations:** ~13,000
- **CUT&RUN role:** Adapter trimming. CUT&RUN sub-nucleosomal fragments frequently read through into adapters, making trimming important for accurate alignment.