ENCODE RNA-seq Pipeline

When to Use

• User wants to run an RNA-seq processing pipeline from FASTQ to gene quantification
• User asks about "RNA-seq pipeline", "STAR alignment", "RSEM", "gene expression quantification", or "Kallisto"
• User needs to process bulk RNA-seq data with ENCODE-standard 2-pass STAR alignment
• Example queries: "process my RNA-seq FASTQs", "quantify gene expression from RNA-seq", "run STAR and RSEM on my data"

Execute the ENCODE RNA-seq processing pipeline from raw FASTQ files through splice-aware alignment, gene/transcript quantification, and strand-specific signal track generation. This skill provides a complete Nextflow DSL2 implementation following ENCODE uniform analysis standards.

Overview

RNA-seq measures transcriptome-wide gene expression by sequencing cDNA derived from cellular RNA. The ENCODE pipeline processes RNA-seq data through quality control, splice-aware alignment with STAR (2-pass mode), gene and transcript quantification with RSEM, optional fast pseudoalignment with Kallisto, and generation of strand-specific signal tracks as bigWig files.

Key design decisions: STAR 2-pass mode for maximum splice junction sensitivity, RSEM for accurate gene/transcript/isoform quantification including multi-mapped reads, stranded library protocol (dUTP/rf-stranded) as the ENCODE standard, and paired-end sequencing with a minimum of 30 million uniquely mapped reads per replicate.

Key Literature

Pipeline Stages

FASTQ
  ├─> FastQC (raw reads)
  └─> Trim Galore (+ FastQC on the trimmed reads)
        ├─> Kallisto (optional) ──────────> kallisto/<sample>/abundance.tsv
        └─> STAR (2-pass)
              ├─> transcriptome BAM ─> RSEM ─> <sample>.genes.results / .isoforms.results
              ├─> bedGraph str1/str2 ─> bedGraphToBigWig ─> signal/<sample>_{plus,minus}.bw
              └─> genome BAM ─> RSeQC (infer_experiment, read_distribution,
                                       geneBody_coverage, inner_distance)

MultiQC <── FastQC (raw + trimmed), trimming reports, STAR Log.final.out,
            RSEM .stat/, RSeQC infer_experiment + read_distribution
  └─> qc/multiqc/multiqc_report.html

The Kallisto abundances, the gene body coverage and the inner distance files are published but are not part of the MultiQC report; read those files directly.

Stage Summary

Input Requirements

Required Files

• RNA-seq FASTQ: paired-end reads matched by the --reads glob (ENCODE standard;

single-end with --single_end)

• STAR genome index directory (--star_index)
• RSEM reference prefix (--rsem_index) produced by rsem-prepare-reference
• BED12 gene model for RSeQC (--rseqc_bed)
• Kallisto index file (--kallisto_index) built with kallisto 0.50.1, unless

--skip_kallisto is set. kallisto 0.50.1 writes index version 13 and rejects an index built with 0.48 or earlier

There is no sample sheet: samples are the pairs that --reads matches, and the sample ID is the shared prefix of each pair. The gene annotation is not a workflow parameter and there is no --gtf; the GTF is consumed when the STAR and RSEM references are built (references/02 and references/03), so the annotation is fixed by the index you pass in.

Library Strandedness

--strandedness takes one value for the whole run — reverse (default), forward or none — and is validated before the first task. It drives three things at once: the RSEM --strandedness flag, the kallisto strand flag, and which STAR bedGraph becomes which signal track.

There is no per-sample strandedness and the workflow does not detect it. RSeQC infer_experiment.py runs as a post-hoc check and writes qc/rseqc/<sample>.infer_experiment.txt. If the library type is unknown, run a first pass, read that file (references/05 explains the output), and rerun with the correct --strandedness — the RSEM counts, the kallisto abundances and the signal tracks all depend on it, so a wrong value has to be corrected by rerunning, not by post-processing.

QC Thresholds

Not computed by this workflow: rRNA rate, library duplication rate (beyond the sequence-level estimate inside the FastQC report), detected-gene counts, and saturation curves. references/05-qc-metrics.md gives the commands to run those by hand on the published BAM and RSEM output.

Read Depth Guidelines

Execution

The versions the workflow runs are the ones in scripts/Dockerfile: STAR 2.7.11b, RSEM 1.3.3, kallisto 0.50.1, samtools 1.19, RSeQC 5.0.3, Trim Galore 0.6.10, cutadapt 4.6, MultiQC 1.21 and FastQC 0.12.1. The conda environment in bioinformatics-installer (environments/rnaseq-env.yml) is a separate manual route pinned to the same versions of STAR, RSEM, kallisto, samtools, RSeQC, Trim Galore, FastQC and MultiQC; it leaves cutadapt to the Trim Galore package, adds salmon and subread, and does not carry bedGraphToBigWig, which the signal-track step needs.

Every index flag is shown in the examples below because the defaults are bare names resolved in the launch directory (GRCh38_star_index, GRCh38_rsem_index/GRCh38, gencode.v38.kallisto.idx, hg38_RefSeq.bed for --genome GRCh38). The run stops before the first task if any of them is missing. --kallisto_index is the one exception: it is not read when --skip_kallisto is set.

Quick Start (local, Docker)

nextflow run scripts/main.nf -profile local \
    --reads 'fastq/*_R{1,2}.fq.gz' \
    --genome GRCh38 \
    --star_index /ref/GRCh38_star_index \
    --rsem_index /ref/GRCh38_rsem_index/GRCh38 \
    --kallisto_index /ref/gencode.v38.kallisto.idx \
    --rseqc_bed /ref/hg38_RefSeq.bed \
    --strandedness reverse \
    --outdir results/

SLURM HPC

nextflow run scripts/main.nf -profile slurm \
    --container /path/to/pipeline-rnaseq.sif \
    --slurm_queue normal \
    --reads 'fastq/*_R{1,2}.fq.gz' \
    --genome GRCh38 \
    --star_index /ref/GRCh38_star_index \
    --rsem_index /ref/GRCh38_rsem_index/GRCh38 \
    --kallisto_index /ref/gencode.v38.kallisto.idx \
    --rseqc_bed /ref/hg38_RefSeq.bed \
    --outdir results/

Cloud

# Google Cloud Batch
nextflow run scripts/main.nf -profile gcp \
    --container us-docker.pkg.dev/<project>/<repo>/pipeline-rnaseq:1.0.0 \
    --gcp_project <project> \
    --gcp_workdir gs://<bucket>/work \
    --reads 'gs://<bucket>/fastq/*_R{1,2}.fq.gz' \
    --genome GRCh38 \
    --star_index gs://<bucket>/ref/GRCh38_star_index \
    --rsem_index gs://<bucket>/ref/GRCh38_rsem_index/GRCh38 \
    --kallisto_index gs://<bucket>/ref/gencode.v38.kallisto.idx \
    --rseqc_bed gs://<bucket>/ref/hg38_RefSeq.bed \
    --outdir gs://<bucket>/results

# AWS Batch
nextflow run scripts/main.nf -profile aws \
    --container <account>.dkr.ecr.<region>.amazonaws.com/pipeline-rnaseq:1.0.0 \
    --aws_queue <job-queue> \
    --aws_workdir s3://<bucket>/work \
    --reads 's3://<bucket>/fastq/*_R{1,2}.fq.gz' \
    --genome GRCh38 \
    --star_index s3://<bucket>/ref/GRCh38_star_index \
    --rsem_index s3://<bucket>/ref/GRCh38_rsem_index/GRCh38 \
    --kallisto_index s3://<bucket>/ref/gencode.v38.kallisto.idx \
    --rseqc_bed s3://<bucket>/ref/hg38_RefSeq.bed \
    --outdir s3://<bucket>/results

--outdir only sets where results are published; Google Batch and AWS Batch stage every task through the work directory, and the workflow stops with an error if it or the project/queue is missing.

Pipeline Parameters

Infrastructure parameters (nextflow.config)

Cloud Cost Estimates

Memory note: nextflow.config asks for 36 GB for STAR_ALIGN (multiplied by the attempt number on a retry, capped by --max_memory). The local executor refuses the task on a machine with less, so a 32 GB host is not enough. If 36 GB is out of reach, rebuild the STAR index with a larger --genomeSAsparseD (which shrinks the loaded index at some cost in mapping speed) or move to a bigger machine. --limitGenomeGenerateRAM is a genomeGenerate option and is not a parameter of this workflow.

Output Directory Structure

results/
  fastqc/                                        # FastQC HTML/zip for raw and trimmed reads
  trimmed/
    <sample>_R1_val_1.fq.gz, <sample>_R2_val_2.fq.gz
    <sample>_R1.fq.gz_trimming_report.txt        # one per input file
  star/
    <sample>.Aligned.sortedByCoord.out.bam
    <sample>.Aligned.sortedByCoord.out.bam.bai
    <sample>.Aligned.toTranscriptome.out.bam     # RSEM input
    <sample>.Log.final.out
    <sample>.SJ.out.tab
    <sample>.ReadsPerGene.out.tab                # STAR gene counts
    <sample>.Signal.UniqueMultiple.str1.out.bg   # plus str2 unless --strandedness none
    <sample>.Signal.Unique.str1.out.bg           # same signal, unique mappers only
  rsem/
    <sample>.genes.results                       # gene_id, TPM, FPKM, expected_count
    <sample>.isoforms.results                    # transcript_id, TPM, FPKM, IsoPct
    <sample>.stat/                               # RSEM model statistics, read by MultiQC
  kallisto/                                      # absent with --skip_kallisto
    <sample>/abundance.tsv
    <sample>/abundance.h5                        # only from a kallisto build with HDF5 support
    <sample>/run_info.json
  signal/
    <sample>_plus.bw, <sample>_minus.bw          # or <sample>_unstranded.bw
  qc/
    rseqc/
      <sample>.infer_experiment.txt
      <sample>.read_distribution.txt
      <sample>.geneBody_coverage.*
      <sample>.inner_distance.*                  # paired-end only
    multiqc/
      multiqc_report.html
      multiqc_data/
  pipeline_info/
    timeline.html, report.html, trace.txt        # Nextflow execution reports

Common Pitfalls

1. Insufficient Memory for STAR

STAR loads the whole genome index into memory. STAR_ALIGN is configured for 36 GB, so a 32 GB machine will not schedule the task at all under -profile local. On shared HPC systems, check the per-job memory limit before submitting.

2. Wrong Strandedness Setting

Using incorrect strandedness results in near-zero gene counts. If you see uniformly low counts, read qc/rseqc/<sample>.infer_experiment.txt and rerun with the matching --strandedness. ENCODE dUTP libraries are reverse stranded.

3. Using FPKM for Cross-Sample Comparison

FPKM values are not comparable across samples because they depend on total library composition. Use TPM (comparable across samples) or raw counts with DESeq2/edgeR normalization for differential expression.

4. Ignoring Multi-Mapped Reads

RSEM uses an expectation-maximization algorithm to probabilistically assign multi-mapped reads. This is critical for gene families and repetitive elements. Do not pre-filter multi-mappers before RSEM quantification.

5. Assuming rRNA Contamination Was Checked

High rRNA contamination (>10%) indicates failed rRNA depletion and reduces effective sequencing depth. This workflow does not measure it — run the manual check in references/05-qc-metrics.md against star/<sample>.Aligned.sortedByCoord.out.bam before trusting the quantifications.

6. Not Using 2-Pass Mode for Novel Junctions

STAR 1-pass mode only uses annotated splice junctions. 2-pass mode first discovers novel junctions then re-maps, critical for non-model organisms or samples with extensive alternative splicing. This workflow always runs --twopassMode Basic.

Pipeline Scripts

ENCODE Data Integration

After running on your own data, compare with ENCODE reference:

# Find matching ENCODE RNA-seq experiments
encode_search_experiments(
    assay_title="total RNA-seq",
    organ="pancreas",
    biosample_type="tissue"
)

# Download ENCODE gene quantifications for comparison
encode_batch_download(
    download_dir="/data/encode_reference/",
    output_type="gene quantifications",
    assay_title="total RNA-seq",
    organ="pancreas",
    assembly="GRCh38"
)

# Download ENCODE signal tracks for browser visualization
encode_search_files(
    file_format="bigWig",
    assay_title="total RNA-seq",
    organ="pancreas",
    output_type="signal of unique reads"
)

Pitfalls & Edge Cases

• Strandedness must match library prep: --strandedness feeds RSEM, kallisto and the

signal tracks. Using the wrong value can halve gene counts, assign reads to antisense genes, and swap the plus/minus bigWigs.

• rRNA contamination: rRNA >10% wastes sequencing depth. Ribosomal depletion libraries

should have <5%, poly-A selection libraries <1%. The workflow does not measure it; see references/05-qc-metrics.md for a manual count, or run Picard CollectRnaSeqMetrics outside the container.

• STAR 2-pass mode is required: The first pass discovers novel splice junctions; the

second pass uses them. Single-pass STAR misses tissue-specific or rare splicing events, reducing sensitivity for differential exon usage.

• Gene-level vs transcript-level quantification: RSEM provides transcript-level estimates but gene-level aggregation is more robust for differential expression. Transcript-level analysis requires many more replicates (≥6).
• TPM normalization is not for cross-sample comparison: TPM normalizes within a sample but is NOT appropriate for comparing expression across conditions. Use DESeq2 size factors or TMM normalization for differential expression.
• Batch effects in multi-lab data: RNA-seq is highly sensitive to library prep method, sequencer, and lab. Always check for batch effects with PCA before combining datasets from different sources.

Walkthrough: Processing ENCODE RNA-seq from FASTQ to Gene Quantification

Goal: Process raw RNA-seq FASTQ files through the ENCODE pipeline to generate gene expression quantifications (TPM/FPKM) and signal tracks. Context: The ENCODE RNA-seq pipeline uses STAR 2-pass alignment and RSEM quantification, producing both gene-level and transcript-level expression estimates.

Step 1: Find RNA-seq experiment

encode_get_experiment(accession="ENCSR000CPR")

Expected output:

{
  "accession": "ENCSR000CPR",
  "assay_title": "total RNA-seq",
  "biosample_summary": "K562",
  "assembly": ["GRCh38"],
  "bio_replicate_count": 2,
  "status": "released"
}

Step 2: List FASTQ files

encode_list_files(experiment_accession="ENCSR000CPR", file_format="fastq")

Expected output (a JSON array of file records; fields abridged):

[
  {"accession": "ENCFF200RN1", "file_format": "fastq", "output_type": "reads", "file_size_human": "3.2 GB", "biological_replicates": [1], "status": "released"},
  {"accession": "ENCFF201RN2", "file_format": "fastq", "output_type": "reads", "file_size_human": "3.3 GB", "biological_replicates": [1], "status": "released"}
]

Step 3: Download and name the FASTQs so a read-pair glob can find them

encode_download_files(file_accessions=["ENCFF200RN1", "ENCFF201RN2"], download_dir="/data/rnaseq/fastq")

ENCODE names every FASTQ after its accession (ENCFF200RN1.fastq.gz), with no _R1/_R2 in the name, so the two files of a pair share no prefix and the --reads glob cannot pair them. Link them into the shape the glob expects. Which mate an accession is comes from the ENCODE file record on encodeproject.org, which carries paired_end (1 or 2) and paired_with; the MCP file tools do not return those two fields:

cd /data/rnaseq/fastq
ln -s ENCFF200RN1.fastq.gz k562_rep1_R1.fq.gz
ln -s ENCFF201RN2.fastq.gz k562_rep1_R2.fq.gz

Step 4: Run the RNA-seq pipeline

nextflow run scripts/main.nf -profile local \
    --reads '/data/rnaseq/fastq/k562_*_R{1,2}.fq.gz' \
    --genome GRCh38 \
    --star_index /ref/GRCh38_star_index \
    --rsem_index /ref/GRCh38_rsem_index/GRCh38 \
    --kallisto_index /ref/gencode.v38.kallisto.idx \
    --rseqc_bed /ref/hg38_RefSeq.bed \
    --strandedness reverse \
    --outdir results/

Key pipeline steps:

1. FastQC on the raw reads, quality and adapter trimming with Trim Galore (which also runs FastQC on the trimmed reads)
2. STAR 2-pass alignment (splice-aware), writing the genome BAM, the transcriptome BAM, gene counts and bedGraphs
3. RSEM gene and isoform quantification (TPM, FPKM, expected counts)
4. Kallisto transcript quantification (optional, skipped with --skip_kallisto)
5. Signal track generation (bedGraph to bigWig)
6. RSeQC (infer_experiment, read_distribution, geneBody_coverage, inner_distance) and MultiQC

Step 5: Validate output quality

Step 6: Use expression data with ENCODE epigenomic data

Compare gene expression with enhancer marks:

encode_search_experiments(assay_title="Histone ChIP-seq", biosample_term_name="K562", target="H3K27ac", organism="Homo sapiens")

Interpretation: Genes with high TPM AND nearby H3K27ac peaks have validated enhancer-gene connections. Low expression despite nearby enhancer marks suggests poised or tissue-specific regulation.

Integration with downstream skills

• Gene quantifications feed into -> peak-annotation for expression-validated peak targets
• Expression data connects to -> gtex-expression for tissue comparison
• Processed data feeds into -> compare-biosamples for differential expression analysis
• Pipeline provenance logged by -> data-provenance

Code Examples

1. Find RNA-seq experiments for a tissue

encode_search_experiments(assay_title="total RNA-seq", organ="liver", organism="Homo sapiens")

Expected output:

{
  "results": [
    {"accession": "ENCSR300RNA", "assay_title": "total RNA-seq", "biosample_summary": "liver tissue male adult (54 years)", "status": "released"}
  ],
  "total": 35,
  "limit": 25,
  "offset": 0,
  "has_more": true,
  "next_offset": 25
}

2. Check for existing gene quantifications

encode_list_files(experiment_accession="ENCSR300RNA", file_format="tsv", output_type="gene quantifications", assembly="GRCh38")

Expected output (a JSON array of file records; fields abridged):

[
  {"accession": "ENCFF400GEQ", "file_format": "tsv", "output_type": "gene quantifications", "assembly": "GRCh38", "file_size_human": "5.2 MB", "status": "released"}
]

3. Download expression data

encode_download_files(file_accessions=["ENCFF400GEQ"], download_dir="/data/rnaseq/quantification")

Expected output (fields abridged):

{
  "downloaded": [
    {"accession": "ENCFF400GEQ", "file_path": "/data/rnaseq/quantification/ENCFF400GEQ.tsv", "file_size_human": "5.2 MB", "success": true, "md5_verified": true}
  ],
  "errors": [],
  "summary": {"total_requested": 1, "successful": 1, "failed": 0, "total_size_human": "5.2 MB"}
}

Integration

Related Skills

• pipeline-guide (parent): General pipeline selection and resource assessment
• quality-assessment: Deep-dive QC analysis beyond basic metrics
• integrative-analysis: Combine RNA-seq with ChIP-seq/ATAC-seq for regulatory inference
• compare-biosamples: Compare expression profiles across cell types
• single-cell-encode: For scRNA-seq data processing (different pipeline)
• pipeline-chipseq: Sibling pipeline for ChIP-seq data
• pipeline-atacseq: Sibling pipeline for ATAC-seq data
• publication-trust: Verify literature claims backing analytical decisions

Presenting Results

When reporting RNA-seq pipeline results:

• Mapping rate: Report the STAR uniquely mapped rate (>70% expected), multi-mapped rate (<10%), and unmapped rate from star/<sample>.Log.final.out
• Quantification paths: Provide paths to rsem/<sample>.genes.results (TPM, FPKM, expected_count), rsem/<sample>.isoforms.results, and kallisto/<sample>/abundance.tsv when Kallisto ran
• Strandedness: Confirm the orientation reported in qc/rseqc/<sample>.infer_experiment.txt matches the --strandedness value the run used; if it does not, the run has to be repeated
• Key QC metrics: Present the exonic rate from read_distribution.txt and the gene body coverage uniformity from geneBody_coverage.geneBodyCoverage.txt in a summary table, alongside the MultiQC report at qc/multiqc/multiqc_report.html
• Derived metrics: Detected-gene counts (TPM>1), rRNA rate, library duplication rate and saturation are not produced by this workflow. Compute them separately if they are needed, and say so when reporting
• Signal tracks: Provide paths to signal/<sample>_plus.bw and signal/<sample>_minus.bw (or signal/<sample>_unstranded.bw for an unstranded run)
• Next steps: Suggest integrative-analysis to combine RNA-seq with ChIP-seq/ATAC-seq for regulatory inference, or compare-biosamples for cross-tissue expression comparison

For the request: "$ARGUMENTS"