RNA-seq expression

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 goal

Commands to run, in order

Inputs

Inspect raw-read quality only

1

input/rnaseq/protocol_example.rnaseq.fastq.list.txt
input/rnaseq/protocol_example.rnaseq.fastq/

Align untrimmed reads with STAR

1 → 3

FASTQ 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 alignment

1 → 2 → 3

Same inputs as the STAR route

Estimate gene-level expression

1 → 3 → 4

STAR-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 expression

1 → 3 → 5

STAR-route inputs
input/rnaseq/protocol_example.rnaseq.bam/sample_bam_list.txt
input/reference_data/RSEM_Index/

QC and normalize existing expression matrices

6 → 7

input/rnaseq/protocol_example.rnaseq.tpm_matrix.bed
input/rnaseq/protocol_example.rnaseq.count_matrix.bed
input/rnaseq/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 input/rnaseq/protocol_example.rnaseq.fastq.list.txt \
    --data-dir input/rnaseq/protocol_example.rnaseq.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 input/rnaseq/protocol_example.rnaseq.fastq.list.txt \
    --data-dir input/rnaseq/protocol_example.rnaseq.fastq --STAR-index input/reference_data/STAR_Index/ \
    --gtf input/reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.gtf \
    --reference-fasta input/reference_data/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta \
    --ref-flat input/reference_data/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 input/rnaseq/protocol_example.rnaseq.fastq.list.txt \
    --data-dir input/rnaseq/protocol_example.rnaseq.fastq --STAR-index input/reference_data/STAR_Index/ \
    --gtf input/reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.gtf \
    --reference-fasta input/reference_data/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta \
    --ref-flat input/reference_data/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 input/rnaseq/protocol_example.rnaseq.fastq.list.txt \
    --data-dir input/rnaseq/protocol_example.rnaseq.fastq \
    --gtf input/reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.collapse_only.gene.gtf \
    --reference-fasta input/reference_data/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta \
    --bam_list input/rnaseq/protocol_example.rnaseq.bam/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 input/rnaseq/protocol_example.rnaseq.fastq.list.txt \
    --data-dir input/rnaseq/protocol_example.rnaseq.fastq \
    --STAR-index input/reference_data/STAR_Index/ \
    --gtf input/reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.gtf \
    --reference-fasta input/reference_data/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta \
    --ref-flat input/reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.ref.flat \
    --bam_list input/rnaseq/protocol_example.rnaseq.bam/sample_bam_list.txt \
    --RSEM-index input/reference_data/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 input/rnaseq/protocol_example.rnaseq.tpm_matrix.bed \
    --counts-gct input/rnaseq/protocol_example.rnaseq.count_matrix.bed

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 input/reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.collapse_only.gene.ERCC.gtf  \
    --count-threshold 1 --sample_participant_lookup input/rnaseq/protocol_example.rnaseq.sample_participant_lookup.txt

Output Files#

Step

Output filename and relative path

Description

1

output/rnaseq/fastqc/*_fastqc.html
output/rnaseq/fastqc/*_fastqc.zip

Per-sample FastQC report and archive

2

output/rnaseq/*.trimmed.fq.gz

Adapter-trimmed FASTQ files

3

output/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

4

output/rnaseq/bam/*.rnaseqc.gene_tpm.gct.gz

Gene-level RNA-SeQC TPM matrix

5

output/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

6

output/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

7

output/rnaseq/protocol_example.low_expression_filtered.outlier_removed.tpm.gct.<normalization_method>.expression.bed.gz

Final normalized BED phenotype; the method name is part of the filename

Example QC figures

website/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