# Stage 1: QC and Trimming

## Input
- Raw paired-end FASTQ files (ENCODE standard is paired-end stranded RNA-seq)
- Adapter sequences: TruSeq adapters (Illumina TruSeq Stranded mRNA kit)

## Tools
- **FastQC 0.12.1**: Per-base quality, adapter content, sequence duplication estimate,
  GC bias
- **Trim Galore 0.6.10** (wraps cutadapt 4.6): Adapter trimming + quality filtering

## Key Difference from ATAC-seq / ChIP-seq
RNA-seq uses **Illumina TruSeq** adapters, not Nextera. Trim Galore auto-detects
TruSeq adapters by default. RNA-seq reads tend to be longer (75-150 bp) and have
higher quality than ATAC-seq, but poly-A tails in adapter read-through are common
with short inserts.

## RNA-seq Specific QC Checks
- **Poly-A tail contamination**: Adapter content plot may show poly-A sequences if
  library inserts are shorter than read length. Trim Galore handles this automatically.
- **rRNA contamination**: High duplication rates combined with skewed GC content can
  indicate failed rRNA depletion. The workflow does not measure the rRNA rate; see the
  manual count in `05-qc-metrics.md`.
- **GC bias**: RNA-seq GC content should reflect transcriptome composition, not genome.
  A bimodal GC plot may indicate contamination or degraded RNA.

## Parameters

| Parameter | Default | Notes |
|-----------|---------|-------|
| Quality cutoff | 20 | Phred score minimum |
| Min length | 36 | Longer minimum than ATAC-seq due to longer fragments |
| Adapter | TruSeq (auto-detect) | Illumina TruSeq adapters |
| Stringency | 1 | Overlap with adapter sequence required (Trim Galore default) |

## Commands

```bash
# Raw QC
fastqc -t 4 -o qc_raw/ sample_R1.fastq.gz sample_R2.fastq.gz

# Paired-end trimming with TruSeq adapters (auto-detected)
trim_galore --paired --quality 20 --length 36 --fastqc \
  --cores 4 -o trimmed/ sample_R1.fastq.gz sample_R2.fastq.gz
```

## Expected Output

The workflow runs FastQC on the raw reads and Trim Galore with `--fastqc` on top, so both
sets of reports exist:

- `trimmed/*_trimming_report.txt` -- trimming statistics (one per input file)
- `trimmed/*_val_1.fq.gz`, `trimmed/*_val_2.fq.gz` -- trimmed paired-end reads
  (`*_trimmed.fq.gz` with `--single_end`)
- `fastqc/*.{html,zip}` -- FastQC reports for the raw reads and for the trimmed reads,
  both published to the same directory and both fed to MultiQC

## QC Checkpoints

| Check | Threshold | Action if Failed |
|-------|-----------|------------------|
| Per-base quality | >Q20 after trimming | Check sequencing run quality |
| Adapter content | <5% after trimming | Verify adapter detection |
| GC content | Unimodal, ~50% for human transcriptome | Check for contamination or rRNA |
| Duplication rate | <60% (variable in RNA-seq) | High duplication is common for abundant transcripts |
| Read count | >=30M PE reads recommended | May need more sequencing |

## Troubleshooting
- **High adapter content**: Short library inserts cause adapter read-through. Ensure
  `--length 36` retains enough reads after trimming.
- **Bimodal GC content**: Likely rRNA contamination or DNA contamination. Confirm with the
  manual rRNA count in `05-qc-metrics.md` and the intronic/intergenic fractions in
  `qc/rseqc/<sample>.read_distribution.txt`.
- **Very high duplication**: Some duplication is biological (highly expressed genes).
  True PCR duplicates are only a concern if NRF <0.5. Mark but do not remove duplicates
  for RNA-seq quantification.

# Stage 2: STAR Alignment

## Tools
- **STAR 2.7.11b**: Splice-aware aligner (Dobin et al. 2013, ~12,000 citations)
- **Samtools 1.19**: BAM conversion, sorting, indexing, and statistics

## Why STAR Instead of HISAT2 or Bowtie2
STAR is preferred for ENCODE RNA-seq because:
1. Native splice junction awareness with 2-pass mode for novel junction discovery
2. Simultaneous output of genome BAM + transcriptome BAM (for RSEM)
3. Built-in gene count quantification (`--quantMode GeneCounts`)
4. Native bedGraph signal generation (`--outWigType bedGraph`)
5. Chimeric read detection for gene fusion analysis
6. Standard in the ENCODE RNA-seq uniform pipeline

## Reference Genome and Annotation

| Organism | Assembly | STAR Index Source | GTF Annotation |
|----------|----------|-------------------|----------------|
| Human | GRCh38 (hg38) | Build from ENCODE FASTA | GENCODE v38+ (comprehensive) |
| Mouse | mm10 (GRCm38) | Build from ENCODE FASTA | GENCODE vM27+ (comprehensive) |

Use the ENCODE "no alt" analysis set FASTA with GENCODE comprehensive gene annotation.
The GTF is consumed here, when the index is built, and is not a parameter of the
workflow — the directory produced below is what you pass as `--star_index`, and its
`chrNameLength.txt` is the default `--chrom_sizes`.

## STAR Index Generation (one-time prep, outside the workflow)

```bash
# Generate genome index (run once, requires ~32 GB RAM for human)
STAR --runMode genomeGenerate \
  --genomeDir star_index/ \
  --genomeFastaFiles GRCh38.primary_assembly.genome.fa \
  --sjdbGTFfile gencode.v38.primary_assembly.annotation.gtf \
  --sjdbOverhang 100 \
  --runThreadN 8
```

The `--sjdbOverhang` should ideally be `read_length - 1`. The default of 100 works well
for 101 bp reads and is robust across typical read lengths (75-150 bp).

STAR 2.7.11b, the version in the image, reads genome indexes built with STAR 2.7.4a or
later. An index generated by anything older has to be rebuilt.

## STAR 2-Pass Alignment

```bash
STAR --genomeDir star_index/ \
  --readFilesIn R1.fq.gz R2.fq.gz \
  --readFilesCommand zcat \
  --runThreadN 8 \
  --outSAMtype BAM SortedByCoordinate \
  --outSAMunmapped Within \
  --outFilterMultimapNmax 20 \
  --alignSJoverhangMin 8 \
  --alignSJDBoverhangMin 1 \
  --outFilterMismatchNmax 999 \
  --outFilterMismatchNoverReadLmax 0.04 \
  --alignIntronMin 20 \
  --alignIntronMax 1000000 \
  --alignMatesGapMax 1000000 \
  --quantMode TranscriptomeSAM GeneCounts \
  --twopassMode Basic \
  --outWigType bedGraph \
  --outWigStrand Stranded
```

`--outWigStrand` is the only part of the command the workflow varies: it passes
`Stranded` for `--strandedness reverse|forward` and `Unstranded` for
`--strandedness none`. Everything else is fixed.

The ~32 GB above is the figure for `genomeGenerate`. Alignment is a separate budget: the
`STAR_ALIGN` process is configured for 36 GB (see `nextflow.config`), multiplied by the
attempt number on a retry, up to `--max_memory`.

## Parameters

| Parameter | Value | Notes |
|-----------|-------|-------|
| 2-pass mode | `--twopassMode Basic` | Discovers novel junctions, re-maps in 2nd pass |
| Multi-map max | `--outFilterMultimapNmax 20` | Allow up to 20 alignments (RSEM handles them) |
| SJ overhang | `--alignSJoverhangMin 8` | Minimum overhang for unannotated junctions |
| SJ DB overhang | `--alignSJDBoverhangMin 1` | Minimum overhang for annotated junctions |
| Mismatch rate | `--outFilterMismatchNoverReadLmax 0.04` | Max 4% mismatch rate |
| Intron min | `--alignIntronMin 20` | Minimum intron length |
| Intron max | `--alignIntronMax 1000000` | Maximum intron length (1 Mb) |
| Quant mode | `TranscriptomeSAM GeneCounts` | Output transcriptome BAM + gene counts |
| Signal output | `--outWigType bedGraph --outWigStrand Stranded` | Strand-specific bedGraph (`Unstranded` when `--strandedness none`) |

## Expected Output

All of these are published to `star/`, prefixed with the sample ID:

- `<sample>.Aligned.sortedByCoord.out.bam` (+ `.bai`) -- genome-sorted BAM
- `<sample>.Aligned.toTranscriptome.out.bam` -- transcriptome BAM (input for RSEM)
- `<sample>.ReadsPerGene.out.tab` -- STAR gene counts (column 2: unstranded, 3: sense, 4: antisense)
- `<sample>.SJ.out.tab` -- splice junction table (novel + annotated)
- `<sample>.Signal.UniqueMultiple.str1.out.bg` -- fragments whose read 1 maps to the + strand
- `<sample>.Signal.UniqueMultiple.str2.out.bg` -- the remaining fragments (absent when `--strandedness none`)
- `<sample>.Signal.Unique.str1.out.bg` (and `str2`) -- the same signal counting only
  uniquely mapped reads
- `<sample>.Log.final.out` -- alignment summary statistics

Which transcript strand str1 and str2 correspond to depends on the library; see
`04-signal-tracks.md`. Both the `UniqueMultiple` and the `Unique` bedGraphs are published
to `star/`, but only the `UniqueMultiple` pair is converted to bigWig.

## QC Checkpoints

| Check | Threshold | Action if Failed |
|-------|-----------|------------------|
| Uniquely mapped reads | >=70% | Check genome build, contamination, rRNA |
| Multi-mapped reads | <10% | High multi-map suggests repetitive contamination |
| Unmapped: too short | <10% | High rate suggests aggressive trimming or poor quality |
| Unmapped: too many mismatches | <5% | Check for contamination (wrong organism) |
| Chimeric reads | 0 unless chimeric detection is enabled | The workflow passes no `--chim*` options |
| Splice junctions (novel) | Thousands expected | Very few suggests annotation mismatch |

## Notes
- STAR 2-pass mode is slower than 1-pass but discovers 10-20% more novel splice junctions.
- The transcriptome BAM (`Aligned.toTranscriptome.out.bam`) must NOT be coordinate-sorted;
  RSEM requires it in its native transcript-coordinate order.
- `--outFilterMultimapNmax 20` is essential for RSEM. Do not reduce this before RSEM.
- Chimeric/fusion detection is not part of this workflow. To do it, run STAR by hand with
  `--chimSegmentMin 12 --chimJunctionOverhangMin 8 --chimOutType Junctions WithinBAM SoftClip`.

# Stage 3: Gene and Transcript Quantification

## Tools
- **RSEM 1.3.3**: Gene/transcript quantification (Li & Dewey 2011, ~6,000 citations)
- **Kallisto 0.50.1**: Fast pseudoalignment quantification (Bray et al. 2016, ~4,000 citations)

## RSEM Quantification (Primary)

RSEM uses an expectation-maximization (EM) algorithm to probabilistically assign
multi-mapped reads to genes and transcripts, providing accurate quantification even
for overlapping gene families and repetitive elements.

### RSEM Index Preparation (one-time prep, outside the workflow)

```bash
# Prepare RSEM reference (run once)
rsem-prepare-reference --gtf gencode.v38.primary_assembly.annotation.gtf \
  --star GRCh38.primary_assembly.genome.fa rsem_index/GRCh38
```

The last argument is a **prefix**, not a directory, and it is what you pass as
`--rsem_index` (here `rsem_index/GRCh38`). The workflow stages every file whose name
starts with that prefix. The annotation is fixed at this point; there is no `--gtf`
parameter downstream.

### RSEM Quantification

```bash
rsem-calculate-expression \
  --paired-end \
  --bam \
  --no-bam-output \
  --estimate-rspd \
  --strandedness reverse \
  --num-threads 8 \
  Aligned.toTranscriptome.out.bam \
  rsem_index/GRCh38 \
  sample_name
```

### RSEM Output Files

| File | Contents | Key Columns |
|------|----------|-------------|
| `rsem/<sample>.genes.results` | Gene-level quantification | gene_id, transcript_id(s), length, effective_length, expected_count, TPM, FPKM |
| `rsem/<sample>.isoforms.results` | Transcript-level quantification | transcript_id, gene_id, length, effective_length, expected_count, TPM, FPKM, IsoPct |
| `rsem/<sample>.stat/` | Model and read statistics | Directory; parsed by MultiQC |

`--no-bam-output` is passed, so RSEM writes no BAM of its own.

### RSEM Strandedness Flags

The workflow passes its `--strandedness` value straight through to RSEM, so these are the
same three values:

| Library Type | Pipeline and RSEM value | Description |
|-------------|-----------|-------------|
| dUTP / rf-stranded | `--strandedness reverse` | ENCODE standard, the pipeline default |
| fr-stranded | `--strandedness forward` | Directional ligation |
| Unstranded | `--strandedness none` | SMARTer, SMART-Seq2, older protocols |

## Kallisto Quantification (Optional Fast Alternative)

Kallisto uses pseudoalignment (k-mer matching without full alignment) for ultra-fast
transcript quantification. It runs 10-100x faster than STAR+RSEM but does not produce
BAM files or support fusion detection.

In this workflow Kallisto runs on the trimmed FASTQ files, in parallel with (not instead
of) STAR and RSEM. Skip it with `--skip_kallisto`, which also stops `--kallisto_index`
from being read.

### Kallisto Index (one-time prep, outside the workflow)

```bash
# Build Kallisto index from transcriptome FASTA (run once)
kallisto index -i kallisto_index.idx gencode.v38.transcripts.fa
```

Pass the resulting file as `--kallisto_index`. Build it with kallisto 0.50.1, the version
in the image: 0.50.1 writes index version 13, and an index built with kallisto 0.48 or
earlier is rejected at load time. Rebuild rather than reuse an older `.idx`.

### Kallisto Quantification

```bash
kallisto quant \
  -i kallisto_index.idx \
  -o sample/ \
  --rf-stranded \
  -t 8 \
  R1.fq.gz R2.fq.gz
```

The strand flag follows `--strandedness`: `--rf-stranded` for `reverse`, `--fr-stranded`
for `forward`, and no flag for `none`. With `--single_end` the workflow substitutes the
fixed `--single -l 200 -s 20`; if your fragment length distribution differs, rerun
kallisto by hand with the right values.

### Kallisto Output

| File | Contents | Published |
|------|----------|-----------|
| `kallisto/<sample>/abundance.tsv` | transcript_id, length, effective_length, est_counts, tpm | Yes |
| `kallisto/<sample>/run_info.json` | Run metadata and statistics | Yes |
| `kallisto/<sample>/abundance.h5` | Binary HDF5 format (for sleuth) | Only when the kallisto build has HDF5 support; the process declares it as an optional output |

## TPM vs FPKM vs Raw Counts

| Metric | Definition | Cross-Sample Comparable | Use Case |
|--------|-----------|------------------------|----------|
| **Raw counts** | Number of reads/fragments mapped to gene | No | Input for DESeq2/edgeR differential expression |
| **TPM** (Transcripts Per Million) | Counts normalized by gene length then library size | Yes | Cross-sample expression comparison |
| **FPKM** (Fragments Per Kilobase per Million) | Counts normalized by library size then gene length | No | Legacy; avoid for cross-sample comparison |

### When to Use Each

- **Differential expression**: Use raw `expected_count` from RSEM with DESeq2 or edgeR.
  These tools apply their own normalization (median-of-ratios or TMM).
- **Cross-sample comparison**: Use TPM. It sums to 1M per sample, enabling direct comparison.
- **Single-gene reporting**: TPM is appropriate for reporting expression of individual genes.
- **Avoid FPKM**: FPKM does not sum to a constant across samples, making it unreliable
  for cross-sample comparison. TPM is strictly preferred.

## QC Checkpoints

None of these are computed by the workflow; each is a check to run on the published
quantifications (see `05-qc-metrics.md` for the detected-gene command).

| Check | Threshold | Action if Failed |
|-------|-----------|------------------|
| Detected genes (TPM>1) | >12,000 (human) | Check sequencing depth, RNA quality |
| RSEM mapping rate (from `<sample>.stat/`) | >70% of transcriptome BAM reads | Check strandedness setting |
| TPM correlation between replicates | r > 0.95 (Pearson) | Check batch effects, sample swap |
| Gene count distribution | Log-normal shape expected | Skewed distribution suggests degradation |

# Stage 4: Signal Track Generation

## Tools
- **STAR**: Generates strand-specific bedGraph output during alignment
- **bedGraphToBigWig**: Converts bedGraph to bigWig (UCSC Kent tools)

## Overview

RNA-seq signal tracks display read coverage across the genome for visualization in
genome browsers (UCSC, IGV, WashU Epigenome Browser). For stranded RNA-seq, separate
plus-strand and minus-strand tracks are generated, enabling gene-level visualization
of sense and antisense transcription.

## STAR bedGraph Output

STAR generates bedGraph files directly during alignment. The workflow passes
`--outWigType bedGraph` always, and `--outWigStrand` follows `--strandedness`:

| `--strandedness` | STAR receives | bedGraph files written |
|------------------|---------------|------------------------|
| `reverse`, `forward` | `--outWigStrand Stranded` | `Signal.{Unique,UniqueMultiple}.str1.out.bg` and `...str2.out.bg` |
| `none` | `--outWigStrand Unstranded` | `Signal.{Unique,UniqueMultiple}.str1.out.bg` only |

## What str1 and str2 Mean

STAR assigns a fragment to str1 when **read 1 maps to the + strand** of the genome
(`signalFromBAM.cpp`: `iStrand = ((flag&0x10)>0) == ((flag&0x80)==0)`). str1 is therefore
not "the plus strand" in the transcript sense — which transcript strand it represents
depends on the library:

| `--strandedness` | Library | str1 holds | str2 holds |
|------------------|---------|-----------|-----------|
| `reverse` (dUTP, ENCODE standard) | read 1 is antisense to the transcript | minus-strand transcription | plus-strand transcription |
| `forward` (directional ligation) | read 1 is sense to the transcript | plus-strand transcription | minus-strand transcription |
| `none` | unstranded | all fragments, one track | — |

The workflow applies exactly this mapping, matching ENCODE's `STAR_RSEM.sh`
(`str[1]="-"; str[2]="+"` for the dUTP default), and names the outputs accordingly:

| `--strandedness` | Published tracks |
|------------------|------------------|
| `reverse` | str1 -> `signal/<sample>_minus.bw`, str2 -> `signal/<sample>_plus.bw` |
| `forward` | str1 -> `signal/<sample>_plus.bw`, str2 -> `signal/<sample>_minus.bw` |
| `none` | str1 -> `signal/<sample>_unstranded.bw` |

## Unique vs UniqueMultiple

`Signal.UniqueMultiple.*` includes multi-mapped reads (ENCODE standard) and is what the
workflow converts to bigWig. `Signal.Unique.*` counts only uniquely mapped reads and is a
more conservative estimate; it is published as
`star/<sample>.Signal.Unique.str*.out.bg` but is not converted, so run the commands below
on it yourself if you want that track as a bigWig.

## bedGraph to bigWig Conversion

This is what the workflow runs for a `reverse`-stranded sample:

```bash
# Sort bedGraph (required by bedGraphToBigWig)
sort -k1,1 -k2,2n sample.Signal.UniqueMultiple.str1.out.bg > minus_sorted.bg
sort -k1,1 -k2,2n sample.Signal.UniqueMultiple.str2.out.bg > plus_sorted.bg

# Convert to bigWig
bedGraphToBigWig minus_sorted.bg chrom.sizes sample_minus.bw
bedGraphToBigWig plus_sorted.bg chrom.sizes sample_plus.bw
```

Swap the two strand labels for a `forward`-stranded library.

## Chromosome Sizes File

The workflow uses `--chrom_sizes`, which defaults to `chrNameLength.txt` inside the STAR
index directory. To build one by hand:

```bash
# From the aligned BAM header
samtools view -H sample.Aligned.sortedByCoord.out.bam | \
  grep '@SQ' | awk '{print $2"\t"$3}' | \
  sed 's/SN://;s/LN://' > chrom.sizes

# Or fetch from UCSC
wget https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.chrom.sizes
```

## Signal Normalization

STAR bedGraph output is raw read counts per position, and the workflow publishes the
bigWigs unnormalized. For cross-sample comparison, normalize them yourself afterwards:

| Method | Description | When to Use |
|--------|-------------|-------------|
| Raw | Unnormalized read counts (what this workflow writes) | Single-sample visualization |
| RPM | Reads per million mapped | Cross-sample comparison (simple) |
| RPKM | RPM per kilobase | Length-normalized (rarely needed for signal) |

ENCODE distributes both raw and RPM-normalized bigWig files. For publication, RPM is
standard for cross-sample comparison in browser screenshots.

## Manual Alternative: bamCoverage (deepTools, not run by this workflow)

deepTools is in neither the container image nor `rnaseq-env.yml`, so this needs a
separate install. It gives more control over normalization and strand selection:

```bash
bamCoverage -b sample.Aligned.sortedByCoord.out.bam \
  -o sample.bw \
  --normalizeUsing RPKM \
  --binSize 10 \
  --filterRNAstrand forward \
  --numberOfProcessors 8
```

## Expected Output
- `signal/<sample>_plus.bw` -- plus-strand transcription (stranded runs)
- `signal/<sample>_minus.bw` -- minus-strand transcription (stranded runs)
- `signal/<sample>_unstranded.bw` -- single track when `--strandedness none`

## Notes
- bedGraphToBigWig requires sorted input and a chromosome sizes file.
- For stranded data, always keep separate plus/minus tracks. Combining them
  loses strand information and confounds sense/antisense transcription.
- A swapped plus/minus pair is easy to miss: check a gene you know the orientation of
  in a browser before publishing tracks.
- bigWig files are typically 50-200 MB each, much smaller than BAM files.

# Stage 5: QC Metrics

## Tools
- **RSeQC 5.0.3**: RNA-seq quality control suite (Wang et al. 2012, ~3,500 citations)
- **MultiQC 1.21**: Aggregated QC report generation

The workflow runs four RSeQC modules — `infer_experiment.py`, `read_distribution.py`,
`geneBody_coverage.py` and, for paired-end runs, `inner_distance.py` — all against the
BED12 file given by `--rseqc_bed`, and publishes their output to `qc/rseqc/`. R is not
installed in the image, so the modules that end by calling Rscript write their `.txt` and
`.r` files but no rendered plot.

## RSeQC Modules

### infer_experiment.py (Strandedness Check)

```bash
infer_experiment.py -r hg38_RefSeq.bed -i sample.Aligned.sortedByCoord.out.bam \
  > sample.infer_experiment.txt
```

Determines library strandedness by sampling read orientation relative to annotated
transcripts. This is a post-hoc check: it does not feed back into quantification. If it
disagrees with the `--strandedness` value the run used, rerun the pipeline with the
correct value.

| Output Pattern | Interpretation | `--strandedness` to use |
|----------------|---------------|-------------|
| "1++,1--,2+-,2-+" > 90% | Forward stranded | `forward` |
| "1+-,1-+,2++,2--" > 90% | Reverse stranded (dUTP) | `reverse` |
| ~50/50 split | Unstranded | `none` |

**ENCODE standard**: Expect >90% reverse-stranded reads for dUTP libraries.

### read_distribution.py (Mapping Distribution)

```bash
read_distribution.py -r hg38_RefSeq.bed -i sample.Aligned.sortedByCoord.out.bam \
  > sample.read_distribution.txt
```

Reports fraction of reads mapping to CDS exons, 5' UTR, 3' UTR, introns, and
intergenic regions.

| Region | Expected (mRNA-seq) | Concern Threshold |
|--------|--------------------|--------------------|
| CDS exons | 40-60% | <30% suggests degradation or DNA contamination |
| 5' UTR | 5-10% | <2% suggests 5' degradation |
| 3' UTR | 15-25% | >40% suggests 3' bias (degraded RNA) |
| Introns | 10-25% | >40% suggests DNA contamination or pre-mRNA |
| Intergenic | <5% | >10% suggests DNA contamination |

### geneBody_coverage.py (Gene Body Coverage)

```bash
geneBody_coverage.py -r hg38_RefSeq.bed \
  -i sample.Aligned.sortedByCoord.out.bam -o sample.geneBody_coverage
```

The workflow passes the same `--rseqc_bed` file used by the other modules. A
housekeeping-gene BED (for example `hg38_HouseKeeping.bed`) runs faster and is the usual
choice when calling this module by hand; both are valid inputs.

Reports normalized coverage across gene bodies (5' to 3') in
`<sample>.geneBody_coverage.geneBodyCoverage.txt`. Uniform coverage indicates intact RNA;
strong 3' bias indicates degradation.

| Pattern | Interpretation |
|---------|---------------|
| Uniform (5'/3' ratio 0.7-1.3) | Good RNA quality |
| 3' bias (5'/3' ratio <0.5) | RNA degradation |
| 5' bias (5'/3' ratio >2.0) | Possible oligo-dT priming bias |

### inner_distance.py (Insert Size Distribution)

```bash
inner_distance.py -r hg38_RefSeq.bed \
  -i sample.Aligned.sortedByCoord.out.bam -o sample.inner_distance
```

Paired-end only; skipped when the run uses `--single_end`. Reports the inner distance
between mates across a set of files matching `<sample>.inner_distance.*`. Negative values
indicate overlapping reads (common for short inserts). The peak should match the expected
library insert size (typically 150-300 bp).

### STAR Log Metrics

The `star/<sample>.Log.final.out` from STAR provides critical metrics:

| Metric | Threshold | Notes |
|--------|-----------|-------|
| Uniquely mapped reads % | >=70% | Primary quality indicator |
| Multi-mapped reads % | <10% | High suggests repetitive contamination |
| Unmapped: too short % | <10% | High suggests over-trimming |
| % of reads mapped to multiple loci | <10% | Expected for gene families |
| % of chimeric reads | <1% | Zero unless STAR is run with chimeric detection |
| Number of splices: Total | Millions expected | Low count suggests annotation mismatch |

## Manual Checks (not run by this workflow)

The metrics below are not produced by the pipeline. Run them yourself against the
published BAM and RSEM output when you need them, and label them as separate steps when
reporting.

### rRNA Rate Assessment

```bash
# Count mapped reads overlapping rRNA loci. -L takes a BED file of regions; without it
# samtools would read the file name as a region string and fail.
samtools view -c -F 4 -L rRNA_intervals.bed star/sample.Aligned.sortedByCoord.out.bam
```

Divide by the total mapped count from `Log.final.out` to get the rate.

| rRNA Rate | Interpretation |
|-----------|---------------|
| <5% | Excellent rRNA depletion |
| 5-10% | Acceptable |
| 10-30% | Suboptimal; reduced effective depth |
| >30% | Failed rRNA depletion; consider re-prep |

### Saturation Analysis

```bash
# RSeQC RPKM saturation (installed in the image, but never invoked by the workflow;
# it writes an .r script that needs R to render)
RPKM_saturation.py -r hg38_RefSeq.bed \
  -i star/sample.Aligned.sortedByCoord.out.bam -o sample_saturation
```

Subsamples reads at increasing fractions (5%, 10%, ..., 100%) and measures gene
detection. A plateau indicates sufficient sequencing depth. If the curve is still
rising at 100%, more sequencing is recommended.

### Detected Genes

```bash
# Genes with TPM > 1; the TPM column of RSEM genes.results is column 6
awk 'NR > 1 && $6 > 1 {n++} END {print "Genes with TPM>1:", n}' rsem/sample.genes.results
```

Expect >12,000 for a human sample at ENCODE depth.

### Library Duplication Rate

Picard `MarkDuplicates` and `CollectRnaSeqMetrics` are not in the container image. The
FastQC report does carry a sequence-level duplication estimate, which is not the same
quantity; run Picard separately if you need the alignment-based rate.

## MultiQC Aggregation

```bash
multiqc . -o . -f
```

The workflow feeds MultiQC exactly these inputs: FastQC reports for the raw reads, the
Trim Galore trimming reports, the FastQC reports for the trimmed reads, the STAR
`Log.final.out`, the RSEM `.stat/` directories, and the RSeQC `infer_experiment` and
`read_distribution` outputs. `geneBody_coverage`, `inner_distance` and the Kallisto
output are published but are not part of the report — read those files directly. The
report is written to `qc/multiqc/multiqc_report.html` with the parsed values in
`qc/multiqc/multiqc_data/`.

## Expected Output
- `qc/rseqc/<sample>.infer_experiment.txt` -- strandedness inference
- `qc/rseqc/<sample>.read_distribution.txt` -- mapping distribution by genomic feature
- `qc/rseqc/<sample>.geneBody_coverage.geneBodyCoverage.txt` and `.r` -- gene body coverage
  (RSeQC also tries to draw `.curves.pdf`, which needs R and is therefore absent)
- `qc/rseqc/<sample>.inner_distance.*` -- insert size, paired-end runs only (RSeQC also
  tries to draw a PDF, which needs R)
- `qc/multiqc/multiqc_report.html` -- aggregated QC report

# RNA-seq Pipeline — Literature References

**Last updated:** 2026-03-07
**Purpose:** Reference catalog for the pipeline-rnaseq skill — papers defining ENCODE RNA-seq processing standards, splice-aware alignment, transcript quantification, and quality assessment tools.

---

## ENCODE Pipeline Standards

---

### ENCODE Project Consortium 2020 — Expanded encyclopaedias of DNA elements

- **Citation:** ENCODE Project Consortium et al. Expanded encyclopaedias of DNA elements in the human and mouse genomes. Nature, 583(7818):699-710, 2020.
- **DOI:** [10.1038/s41586-020-2493-4](https://doi.org/10.1038/s41586-020-2493-4)
- **PMID:** 32728249 | **PMC:** PMC7410828
- **Citations:** ~2,500
- **Key findings:** ENCODE Phase 3 paper establishing RNA-seq as a core assay for the encyclopedia. RNA-seq quantification provides the expression data used to classify candidate cis-Regulatory Elements (cCREs) and annotate gene activity across cell types. Established uniform processing standards: stranded library protocol (dUTP), STAR 2-pass alignment, RSEM quantification, minimum 30M uniquely mapped reads per replicate, and GENCODE annotation as the reference gene model. Assembly standards: GRCh38/mm10.

---

### Hitz et al. 2023 — The ENCODE Uniform Analysis Pipelines

- **Citation:** Hitz BC, Lee JW, Jolanki O, et al. The ENCODE Uniform Analysis Pipelines. Research Square (preprint), 2023.
- **DOI:** [10.21203/rs.3.rs-311127/v1](https://doi.org/10.21203/rs.3.rs-311127/v1)
- **Citations:** ~84
- **Key findings:** Documents the official ENCODE RNA-seq pipeline specification: STAR 2-pass mode for splice-aware alignment, RSEM for gene/transcript quantification, strand-specific signal track generation, and comprehensive QC with RSeQC. Pipeline infrastructure uses Docker containers and WDL for reproducible execution on local machines, HPC clusters, or cloud environments.

---

### Conesa et al. 2016 — A survey of best practices for RNA-seq data analysis

- **Citation:** Conesa A, Madrigal P, Tarazona S, Gomez-Cabrero D, Cervera A, McPherson A, Szczesniak MW, Gaffney DJ, Elo LL, Zhang X, Mortazavi A. A survey of best practices for RNA-seq data analysis. Genome Biology, 17:13, 2016.
- **DOI:** [10.1186/s13059-016-0881-8](https://doi.org/10.1186/s13059-016-0881-8)
- **PMID:** 26813401 | **PMC:** PMC4728800
- **Citations:** ~4,500
- **Key findings:** Comprehensive best-practices review covering experimental design, quality control, alignment, quantification, normalization, differential expression, alternative splicing, functional analysis, and visualization. Recommended minimum sequencing depth of 10-30M reads for differential expression and >60M for transcript discovery. Established that RSEM with STAR alignment provides the most accurate quantification for multi-mapped reads through expectation-maximization. Defined key QC metrics: mapping rate >80%, rRNA contamination <10%, replicate Pearson correlation >= 0.9, and 3'/5' coverage bias assessment. This paper provides the scientific basis for many parameter choices in the ENCODE RNA-seq pipeline.

---

## Core Pipeline Tools

---

### Dobin et al. 2013 — STAR: ultrafast universal RNA-seq aligner

- **Citation:** Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics, 29(1):15-21, 2013.
- **DOI:** [10.1093/bioinformatics/bts635](https://doi.org/10.1093/bioinformatics/bts635)
- **PMID:** 23104886 | **PMC:** PMC3530905
- **Citations:** ~12,700
- **Key findings:** Introduced the STAR (Spliced Transcripts Alignment to a Reference) aligner, developed specifically for the ENCODE Transcriptome project to align >80 billion reads. STAR uses sequential maximum mappable seed search in uncompressed suffix arrays followed by seed clustering and stitching. Outperforms other aligners by >50x in mapping speed (550M 2x76bp paired-end reads per hour on 12 cores) while improving alignment sensitivity and precision. Supports de novo detection of canonical and non-canonical splice junctions, chimeric transcripts, and full-length RNA sequences. The ENCODE pipeline uses STAR in 2-pass mode: first pass discovers novel splice junctions, second pass uses the full junction database for improved alignment. Experimentally validated 1,960 novel intergenic splice junctions with 80-90% success rate.

---

### Li & Dewey 2011 — RSEM: accurate transcript quantification

- **Citation:** Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics, 12:323, 2011.
- **DOI:** [10.1186/1471-2105-12-323](https://doi.org/10.1186/1471-2105-12-323)
- **PMID:** 21816040 | **PMC:** PMC3163565
- **Citations:** ~18,600
- **Key findings:** Introduced RSEM (RNA-Seq by Expectation-Maximization), which uses a generative model and EM algorithm to handle multi-mapped reads for accurate gene and isoform quantification. RSEM outputs expected counts, TPM (Transcripts Per Million), and FPKM (Fragments Per Kilobase of transcript per Million mapped reads) along with 95% credibility intervals. Unlike simple counting methods (e.g., featureCounts), RSEM probabilistically assigns ambiguously mapped reads to their most likely transcripts of origin based on the overall expression landscape. Benchmarks showed RSEM slightly outperforms other quantification pipelines (Teng et al. 2016). The ENCODE pipeline uses RSEM with STAR's transcriptome BAM output for gene and transcript quantification.

---

### Bray et al. 2016 — Kallisto: near-optimal RNA-seq quantification

- **Citation:** Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nature Biotechnology, 34(5):525-527, 2016.
- **DOI:** [10.1038/nbt.3519](https://doi.org/10.1038/nbt.3519)
- **PMID:** 27043002
- **Citations:** ~4,000
- **Key findings:** Introduced Kallisto, which uses pseudoalignment to a de Bruijn graph built from a transcriptome reference to quantify transcript abundances without full read alignment. Kallisto processes 30M reads in minutes (vs hours for alignment-based methods) while maintaining accuracy comparable to or better than alignment-based quantification. The ENCODE pipeline includes Kallisto as an optional fast pseudoalignment step that operates directly on FASTQ files without requiring genome alignment, providing an independent quantification for cross-validation with RSEM results.

---

## Quality Control

---

### Wang et al. 2012 — RSeQC: quality control of RNA-seq experiments

- **Citation:** Wang L, Wang S, Li W. RSeQC: quality control of RNA-seq experiments. Bioinformatics, 28(16):2184-2185, 2012.
- **DOI:** [10.1093/bioinformatics/bts356](https://doi.org/10.1093/bioinformatics/bts356)
- **PMID:** 22743226
- **Citations:** ~2,300
- **Key findings:** Comprehensive RNA-seq QC package evaluating sequence quality, GC bias, PCR bias, nucleotide composition bias, sequencing depth saturation, strand specificity, coverage uniformity, and read distribution over genome structure (exonic, intronic, intergenic). Key modules used in this pipeline: infer_experiment.py (strand specificity detection), read_distribution.py (exonic/intronic/intergenic classification), geneBody_coverage.py (3'/5' bias assessment), junction_saturation.py (splice junction discovery saturation), and inner_distance.py (fragment size estimation). RSeQC is the primary RNA-seq-specific QC tool in the ENCODE pipeline.

---

### Ewels et al. 2016 — MultiQC

- **Citation:** Ewels P, Magnusson M, Lundin S, Käller M. MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics, 32(19):3047-3048, 2016.
- **DOI:** [10.1093/bioinformatics/btw354](https://doi.org/10.1093/bioinformatics/btw354)
- **PMID:** 27312411 | **Citations:** ~6,800
- **RNA-seq role:** Aggregates QC metrics from FastQC, STAR, RSEM, RSeQC, and Picard into a unified HTML report for batch-level assessment of RNA-seq experiments.

---

## Gene Annotation

---

### Frankish et al. 2021 — GENCODE 2021

- **Citation:** Frankish A, Diekhans M, Jungreis I, Lagarde J, Loveland JE, Mudge JM, Sisu C, Wright JC, Armstrong J, Barnes I, Berry A, Bignell A, Boix C, Carbonell Sala S, Cunningham F, Di Domenico T, Donaldson S, Fiddes IT, Garcia Giron C, Gonzalez JM, Grego T, Hardy M, Hourlier T, Howe KL, Hunt T, Izuogu OG, Johnson R, Martin FJ, Martinez L, Mohanan S, Muir P, Navarro FCP, Parker A, Pei B, Pozo F, Riera FC, Ruffier M, Schmitt BM, Stapleton E, Suner MM, Sycheva I, Uszczynska-Ratajczak B, Wolf MY, Xu J, Yang YT, Yates A, Zerbino D, Zhang Y, Choudhary JS, Gerstein M, Guigo R, Hubbard TJP, Kellis M, Paten B, Tress ML, Flicek P. GENCODE 2021. Nucleic Acids Research, 49(D1):D916-D923, 2021.
- **DOI:** [10.1093/nar/gkaa1087](https://doi.org/10.1093/nar/gkaa1087)
- **PMID:** 33270111 | **PMC:** PMC7778937
- **Citations:** ~2,500
- **Key findings:** GENCODE provides the comprehensive gene annotation used as the reference for RNA-seq quantification in ENCODE pipelines. Version 41 (GRCh38.p13) annotates 62,764 genes (20,090 protein-coding) and 248,396 transcripts. GENCODE annotations are distinguished from RefSeq by their inclusion of more alternative transcripts, non-coding RNA species, and pseudogenes. The ENCODE pipeline requires GENCODE GTF files for both STAR genome index generation and RSEM reference preparation, ensuring consistent annotation across all quantification steps.

---

## Supplementary Tools (Non-RNA-seq-Specific)

See pipeline-chipseq/references/literature.md for detailed descriptions of shared tools.

---

### Li et al. 2009 — SAMtools

- **DOI:** [10.1093/bioinformatics/btp352](https://doi.org/10.1093/bioinformatics/btp352) | **PMID:** 19505943 | **Citations:** ~53,700
- **RNA-seq role:** BAM sorting, indexing, and alignment statistics. Used for computing mapping rate and rRNA fraction.

---

### Broad Institute — Picard MarkDuplicates

- **URL:** [https://broadinstitute.github.io/picard/](https://broadinstitute.github.io/picard/)
- **RNA-seq role:** PCR duplicate marking. Note: duplicate removal is less critical for RNA-seq than for ChIP-seq because PCR duplicates in expression data do not create false peaks. ENCODE retains duplicates for quantification but marks them for QC statistics.

---

### Andrews 2010 — FastQC

- **URL:** [https://www.bioinformatics.babraham.ac.uk/projects/fastqc/](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)
- **RNA-seq role:** Raw read quality assessment. RNA-seq libraries may show characteristic GC bias from highly expressed transcripts and sequence duplication from abundant mRNAs — these are expected and not concerning.

---

### 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
- **RNA-seq role:** Adapter trimming with Phred quality cutoff of 20. Less critical for RNA-seq than for ATAC-seq since RNA-seq fragments are typically longer than read length, but still recommended for best alignment quality.

---

### Teng et al. 2016 — A benchmark for RNA-seq quantification pipelines

- **Citation:** Teng M, Love MI, Davis CA, Djebali S, Dobin A, Graveley BR, Li S, Mason CE, Olson S, Pervouchine D, Sloan CA, Wei X, Zhan L, Irizarry RA. A benchmark for RNA-seq quantification pipelines. Genome Biology, 17:74, 2016.
- **DOI:** [10.1186/s13059-016-0940-1](https://doi.org/10.1186/s13059-016-0940-1)
- **PMID:** 27107712 | **PMC:** PMC4842274
- **Citations:** ~162
- **Key findings:** Benchmark of seven RNA-seq quantification pipelines using two independent datasets. Found that performance was generally poor across methods, with RSEM slightly outperforming the rest. Established metrics for evaluating quantification accuracy including specificity (proportion of truly non-expressed genes called as non-expressed) and sensitivity (correlation between estimated and true expression). Supports the ENCODE choice of RSEM as the primary quantification tool.