Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

RNA-seq expression

Process bulk RNA-seq FASTQ files into aligned reads and expression estimates, then perform cohort-level expression quality control and normalization.

Miniprotocol Timing

Timing: TBD

Overview

This mini-protocol walks through bulk RNA-seq expression processing. RNA_calling.ipynb performs read-level quality assessment, optional adapter trimming, STAR alignment, and gene- or transcript-level quantification. bulk_expression_QC.ipynb filters low-expression features and detects sample outliers, and bulk_expression_normalization.ipynb produces normalized expression matrices for downstream xQTL analysis.

The commands form selectable routes rather than one mandatory chain. Gene-level RNA-SeQC and transcript-level RSEM quantification are alternatives after alignment. The bundled matrix-QC example starts from existing count and TPM matrices, so steps 6–7 can be run independently of the toy FASTQ example in steps 1–5.

Steps

Choose a route before running commands.

Analysis goalCommands to run, in orderInputs
Inspect raw-read quality only1tests/fixtures/rna_calling/protocol_example.rnaseq.fastq.list.txt
input/rnaseq/protocol_example.rnaseq.fastq/
Align untrimmed reads with STAR1 → 3FASTQ inputs above
input/reference_data/STAR_Index/
input/reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.gtf
input/reference_data/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta
input/reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.ref.flat
Trim adapters before STAR alignment1 → 2 → 3Same inputs as the STAR route
Estimate gene-level expression1 → 3 → 4STAR-route inputs
input/rnaseq/protocol_example.rnaseq.bam/sample_bam_list.txt
input/reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.collapse_only.gene.gtf
Estimate transcript-level expression1 → 3 → 5STAR-route inputs
input/rnaseq/protocol_example.rnaseq.bam/sample_bam_list.txt
input/reference_data/RSEM_Index/
QC and normalize existing expression matrices6 → 7tests/fixtures/bulk_expression_normalization/protocol_example.rnaseq.tpm.gct.gz
tests/fixtures/bulk_expression_normalization/protocol_example.rnaseq.geneCount.gct.gz
tests/fixtures/bulk_expression_normalization/protocol_example.rnaseq.sample_participant_lookup.txt
input/reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.collapse_only.gene.ERCC.gtf

Run only the row matching the intended analysis goal. Adapter trimming is optional; use it when read-level QC or library preparation indicates adapter contamination.

1. Assess FASTQ read quality

What it does: Generate per-sample FastQC reports before alignment.

sos run pipeline/RNA_calling.ipynb fastqc \
    --cwd output/rnaseq/fastqc \
    --sample-list tests/fixtures/rna_calling/protocol_example.rnaseq.fastq.list.txt \
    --data-dir tests/fixtures/rna_calling/fastq

2. Trim sequencing adapters

What it does: Remove configured adapter sequences and write trimmed FASTQ files for alignment.

sos run pipeline/RNA_calling.ipynb fastp_trim_adaptor \
    --cwd output/rnaseq --sample-list tests/fixtures/rna_calling/protocol_example.rnaseq.fastq.list.txt \
    --data-dir tests/fixtures/rna_calling/fastq --STAR-index <path/to/STAR_Index> \
    --gtf <path/to/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.gtf> \
    --reference-fasta <path/to/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta> \
    --ref-flat <path/to/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.ref.flat>

3. Align reads with STAR and run Picard QC

What it does: Map reads to the reference genome and calculate alignment-level RNA-seq quality metrics.

sos run pipeline/RNA_calling.ipynb STAR_align \
    --cwd output/rnaseq/bam --sample-list tests/fixtures/rna_calling/protocol_example.rnaseq.fastq.list.txt \
    --data-dir tests/fixtures/rna_calling/fastq --STAR-index <path/to/STAR_Index> \
    --gtf <path/to/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.gtf> \
    --reference-fasta <path/to/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta> \
    --ref-flat <path/to/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.ref.flat> \
    --chimSegmentMin 0 \
    -J 50 --mem 200G --numThreads 8

4. Quantify gene-level expression with RNA-SeQC

What it does: Summarize aligned reads into gene-level expression measurements.

sos run pipeline/RNA_calling.ipynb rnaseqc_call \
    --cwd output/rnaseq/bam \
    --sample-list tests/fixtures/rna_calling/protocol_example.rnaseq.fastq.list.txt \
    --data-dir tests/fixtures/rna_calling/fastq \
    --gtf <path/to/Homo_sapiens.GRCh38.103.chr.reformatted.collapse_only.gene.gtf> \
    --reference-fasta <path/to/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta> \
    --bam_list output/rnaseq/fastqc/sample_bam_list.txt

5. Quantify transcript-level expression with RSEM

What it does: Estimate transcript- and gene-level abundance using the RSEM reference index.

sos run pipeline/RNA_calling.ipynb rsem_call \
    --cwd output/rnaseq/bam \
    --sample-list tests/fixtures/rna_calling/protocol_example.rnaseq.fastq.list.txt \
    --data-dir tests/fixtures/rna_calling/fastq \
    --STAR-index <path/to/STAR_Index> \
    --gtf <path/to/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.gtf> \
    --reference-fasta <path/to/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta> \
    --ref-flat <path/to/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.ref.flat> \
    --bam_list output/rnaseq/fastqc/sample_bam_list.txt \
    --RSEM-index <path/to/RSEM_Index>

6. Perform cohort-level expression QC

What it does: Filter low-expression features and identify expression outlier samples across the cohort.

sos run pipeline/bulk_expression_QC.ipynb qc \
    --cwd output/rnaseq \
    --tpm-gct tests/fixtures/bulk_expression_normalization/protocol_example.rnaseq.tpm.gct.gz \
    --counts-gct tests/fixtures/bulk_expression_normalization/protocol_example.rnaseq.geneCount.gct.gz

7. Normalize the QC-passed expression matrices

What it does: Normalize the retained count and TPM matrices and write association-ready expression phenotypes.

sos run pipeline/bulk_expression_normalization.ipynb normalize \
    --cwd output/rnaseq \
    --tpm-gct output/rnaseq/protocol_example.low_expression_filtered.outlier_removed.tpm.gct.gz \
    --counts-gct output/rnaseq/protocol_example.low_expression_filtered.outlier_removed.geneCount.gct.gz \
    --annotation-gtf tests/fixtures/gene_annotation/Homo_sapiens.GRCh38.103.collapse_only.gene.chr22.gtf.gz  \
    --count-threshold 1 --sample_participant_lookup tests/fixtures/bulk_expression_normalization/protocol_example.rnaseq.sample_participant_lookup.txt

Output

StepOutput filename and relative pathDescription
1output/rnaseq/fastqc/*_fastqc.html
output/rnaseq/fastqc/*_fastqc.zip		Per-sample FastQC report and archive
2output/rnaseq/*.trimmed.fq.gzAdapter-trimmed FASTQ files
3output/rnaseq/bam/*.Aligned.sortedByCoord.out.bam
output/rnaseq/bam/*.alignment_summary_metrics
output/rnaseq/bam/*.bam_file_list		STAR alignment, Picard metrics and BAM manifest
4output/rnaseq/bam/*.rnaseqc.gene_tpm.gct.gzGene-level RNA-SeQC TPM matrix
5output/rnaseq/bam/*.rsem.isoforms.results
output/rnaseq/bam/*.rsem.genes.results
output/rnaseq/bam/*.rsem_transcripts_expected_count.txt.gz		RSEM transcript- and gene-level estimates
6output/rnaseq/protocol_example.low_expression_filtered.tpm.gct.gz
output/rnaseq/protocol_example.low_expression_filtered.outlier_removed.tpm.gct.gz
output/rnaseq/protocol_example.low_expression_filtered.outlier_removed.geneCount.gct.gz		Filtered TPM and count matrices after sample QC
7output/rnaseq/protocol_example.low_expression_filtered.outlier_removed.tpm.gct.<normalization_method>.expression.bed.gzFinal normalized BED phenotype; the method name is part of the filename
Example QC figureswebsite/nature_protocol/PCC_sample_list_subset.rnaseqc.low_expression_filtered.outlier_removed.tpm.gct.D_stat_hist.png
website/nature_protocol/PCC_sample_list_subset.rnaseqc.low_expression_filtered.outlier_removed.tpm.gct.RLEplot.png
website/nature_protocol/PCC_sample_list_subset.rnaseqc.low_expression_filtered.outlier_removed.tpm.gct.preQC_cluster.png		Tracked cohort-level examples shown below; the toy-output directory does not currently contain these plots

Anticipated Results

The calling routes produce aligned BAM files and either gene-level or transcript-level expression estimates. The cohort-level route produces QC-passed and normalized expression matrices suitable for phenotype annotation and xQTL association testing.

The following QC figures summarize distributional and sample-level checks applied before normalization.

Figure 1A. Bulk RNA-seq quality-control D-statistic distribution.

Figure 1A. Bulk RNA-seq quality-control D-statistic distribution.

Figure 1B. Bulk RNA-seq quality-control relative log-expression residuals.

Figure 1B. Bulk RNA-seq quality-control relative log-expression residuals.

Figure 1C. Bulk RNA-seq quality-control Mahalanobis-distance p-value clustering.

Figure 1C. Bulk RNA-seq quality-control Mahalanobis-distance p-value clustering.

Command Interface

List the workflows and parameters available in each module used by this mini-protocol.

sos run pipeline/RNA_calling.ipynb -h
sos run pipeline/bulk_expression_QC.ipynb -h
sos run pipeline/bulk_expression_normalization.ipynb -h
Reference data
RSS LD Sketch Pipeline
Molecular Phenotypes
Quantifying expression from RNA-seq data