Bulk RNA-seq counts normalization

Bulk RNA-seq counts normalization#

Filters out lowly expressed genes and normalizes the RNA-seq count matrix into the molecular phenotype BED file used for QTL analysis.

Overview#

The single normalize workflow follows the GTEx pipeline. Genes are kept only when they clear both expression thresholds in a minimum fraction of samples: TPM at or above --tpm-threshold (0.1 by default) in at least --sample-frac-threshold of samples (0.2), and raw read count at or above --count-threshold (6) in at least that same fraction. Both conditions must hold, so a gene that is highly expressed in a handful of samples and absent elsewhere is dropped.

What survives is then rescaled. --normalization-method selects between TMM followed by CPM through voom (tmm_cpm_voom, the default), TMM followed by CPM through edgeR (tmm_cpm_edger), and quantile normalization across samples (qn). Separately, --quantile-normalize is on by default and applies a per-gene inverse-normal transform after rescaling, which is what makes the distribution of each gene comparable across genes rather than across samples. Turning it off adds no_qnorm to the output file name, so the two variants never collide.

Sample identifiers are translated on the way out: --sample-participant-lookup maps expression sample names to the genotype sample names, so the resulting BED can be matched against genotypes directly. Gene IDs are stripped of their Ensembl version suffix.

When to run it. After bulk_expression_QC, which produces the filtered and outlier-removed TPM and count matrices this step consumes, and before covariate preprocessing and the association scan.

Input#

  • --tpm-gct – the TPM matrix in RNA-SeQC GCT format, from bulk_expression_QC, e.g. output/rnaseq/protocol_example.low_expression_filtered.outlier_removed.tpm.gct.gz. The first two lines are the GCT version and matrix dimensions, the third the sample header:

    #1.2
# 300 54
gene_ID	SAMPLE_001	SAMPLE_002	SAMPLE_003	SAMPLE_004	SAMPLE_005	SAMPLE_006	SAMPLE_007	SAM

  • --counts-gct – the raw read-count matrix in the same GCT format and from the same step, e.g. output/rnaseq/protocol_example.low_expression_filtered.outlier_removed.geneCount.gct.gz. It must cover the same genes and samples as the TPM matrix.

  • --annotation-gtf – the collapsed gene model, e.g. reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.collapse_only.gene.ERCC.gtf, from reference_data_preparation. It supplies the chromosome and position written into the output BED, so its gene names must match the matrices and use the chr prefix.

  • --sample-participant-lookup – tab-delimited table with a header and two columns, sample_id and participant_id, mapping expression sample names to genotype sample names; the two may be identical. Every sample in the matrices has to appear here:

    sample_id	participant_id
SAMPLE_001	SAMPLE_001
SAMPLE_002	SAMPLE_002

  • --cwd – output directory, default output.

Filtering and normalization:

  • --tpm-threshold – minimum TPM, default 0.1.

  • --count-threshold – minimum read count, default 6. The example below lowers it to 1 because the toy matrix is small.

  • --sample-frac-threshold – fraction of samples that must clear both thresholds, default 0.2.

  • --normalization-method – tmm_cpm_voom (default), tmm_cpm_edger, or qn.

  • --quantile-normalize / --no-quantile-normalize – per-gene inverse-normal transform after rescaling, on by default.

Runtime:

  • --numThreads, --job-size, --walltime, --mem – threads (default 20) and cluster resources.

  • --container, --entrypoint – software environment.

  • --modular-script-dir – location of the R driver, default code/script.

Output#

  • <tpm_gct>.<normalization_method>.expression.bed.gz – the normalized expression matrix in BED layout: #chr, start, end and gene_id, then one column per sample, named with the genotype IDs from the look-up table. no_qnorm is inserted before expression when --no-quantile-normalize is used. This file is the molecular phenotype consumed by the QTL and TWAS workflows:

    #chr	start	end	gene_id	SAMPLE_001	SAMPLE_002	SAMPLE_003	SAMPLE_004	SAMPLE_005	SAMPLE_006
chr22	10939387	10961338	ENSG00000283047	2.092837798505774	-1.4557760251170175	0.16019524
chr22	15528191	15529139	ENSG00000130538	-1.6022926552158765	-0.2065470244188349	0.498428

The step also writes .stdout and .stderr beside its output. Example results for the toy data are under output/rnaseq/.

Minimal Working Example#

Filter and normalize the QC-ed matrices#

The TPM and count matrices are the output of bulk_expression_QC. --count-threshold 1 relaxes the default of 6, which few genes would clear in a toy matrix this small.

Timing: ~1-2 min (on toy dataset)

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

Command Interface#

sos run pipeline/bulk_expression_normalization.ipynb -h

/home/caox2/.pixi/envs/python/lib/python3.12/site-packages/sos/targets.py:22: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
  import pkg_resources
usage: sos run bulk_expression_normalization.ipynb
               [workflow_name | -t targets] [options] [workflow_options]
  workflow_name:        Single or combined workflows defined in this script
  targets:              One or more targets to generate
  options:              Single-hyphen sos parameters (see "sos run -h" for details)
  workflow_options:     Double-hyphen workflow-specific parameters

Workflows:
  normalize

Global Workflow Options:
  --cwd output (as path)
                        Work directory & output directory
  --counts-gct VAL (as path, required)
                        gene count table
  --tpm-gct VAL (as path, required)
                        gene TPM table
  --annotation-gtf VAL (as path, required)
                        gene gtf annotation table
  --sample-participant-lookup VAL (as path, required)
                        A file to map sample ID from expression to genotype,must
                        contain two columns, sample_id and participant_id,
                        mapping IDs in the expression files to IDs in the
                        genotype (these can be the same).
  --tpm-threshold 0.1 (as float)
  --count-threshold 6 (as int)
  --sample-frac-threshold 0.2 (as float)
  --normalization-method 'tmm_cpm_voom'
                        Normalization method: TMM + CPM(voom) (tmm_cpm_voom),
                        TMM + CPM(edgeR) (tmm_cpm_edger), or quantile
                        normalization (qn)
  --[no-]quantile-normalize (default to True)
  --job-size 1 (as int)
                        For cluster jobs, number commands to run per job
  --walltime 5h
                        Wall clock time expected
  --mem 16G
                        Memory expected
  --numThreads 20 (as int)
                        Number of threads
  --container ''


Sections
  normalize:

usage: sos run pipeline/bulk_expression_normalization.ipynb
               [workflow_name | -t targets] [options] [workflow_options]
  workflow_name:        Single or combined workflows defined in this script
  targets:              One or more targets to generate
  options:              Single-hyphen sos parameters (see "sos run -h" for details)
  workflow_options:     Double-hyphen workflow-specific parameters

Workflows:
  normalize

Global Workflow Options:
  --modular-script-dir code/script (as path)
  --cwd output (as path)
                        Work directory & output directory
  --counts-gct VAL (as path, required)
                        gene count table
  --tpm-gct VAL (as path, required)
                        gene TPM table
  --annotation-gtf VAL (as path, required)
                        gene gtf annotation table
  --sample-participant-lookup VAL (as path, required)
                        A file to map sample ID from expression to genotype,must
                        contain two columns, sample_id and participant_id,
                        mapping IDs in the expression files to IDs in the
                        genotype (these can be the same).
  --tpm-threshold 0.1 (as float)
  --count-threshold 6 (as int)
  --sample-frac-threshold 0.2 (as float)
  --normalization-method 'tmm_cpm_voom'
                        Normalization method: TMM + CPM(voom) (tmm_cpm_voom),
                        TMM + CPM(edgeR) (tmm_cpm_edger), or quantile
                        normalization (qn)
  --[no-]quantile-normalize (default to True)
                        Quantile Normalization after rescale
  --job-size 1 (as int)
                        For cluster jobs, number commands to run per job
  --walltime 5h
                        Wall clock time expected
  --mem 16G
                        Memory expected
  --numThreads 20 (as int)
                        Number of threads
  --container ''
  --entrypoint  ('micromamba run -a "" -n' + ' ' + re.sub(r'(_apptainer:latest|_docker:latest|\.sif)$', '', container.split('/')[-1])) if container else ""


Sections
  normalize:

Workflow implementation#

The SoS workflow definitions below are unchanged from the original protocol.

[global]
parameter: modular_script_dir = path('code/script')  # override with --modular-script-dir
# Work directory & output directory
parameter: cwd = path("output")
#  gene count table
parameter: counts_gct = path
#  gene TPM table
parameter: tpm_gct = path
#  gene gtf annotation table
parameter: annotation_gtf = path
# A file to map sample ID from expression to genotype,must contain two columns, sample_id and participant_id, mapping IDs in the expression files to IDs in the genotype (these can be the same).
parameter: sample_participant_lookup = path
parameter: tpm_threshold = 0.1
parameter: count_threshold = 6
parameter: sample_frac_threshold = 0.2
# Normalization method: TMM + CPM(voom) (tmm_cpm_voom), TMM + CPM(edgeR) (tmm_cpm_edger), or quantile normalization (qn)
parameter: normalization_method = 'tmm_cpm_voom'
# Quantile Normalization after rescale
parameter: quantile_normalize = True
# For cluster jobs, number commands to run per job
parameter: job_size = 1
# Wall clock time expected
parameter: walltime = "5h"
# Memory expected
parameter: mem = "16G"
# Number of threads
parameter: numThreads = 20
parameter: container = ""
import re
parameter: entrypoint= ('micromamba run -a "" -n' + ' ' + re.sub(r'(_apptainer:latest|_docker:latest|\.sif)$', '', container.split('/')[-1])) if container else ""

[normalize]
# Path to the input molecular phenotype data, should be a processd and indexed bed.gz file, with tabix index.
input: tpm_gct, counts_gct, annotation_gtf, sample_participant_lookup
output: f'{cwd:a}/{_input[0]:bnnn}.{normalization_method}.{"no_qnorm." if not quantile_normalize else ""}expression.bed.gz'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime,  mem = mem, tags = f'{step_name}_{_output[0]:bn}'  
bash: expand= "${ }", stderr = f'{_output[0]:nn}.stderr', stdout = f'{_output[0]:nn}.stdout', container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/molecular_phenotypes/QC/bulk_expression_normalization.R \
        --step normalize \
        --cwd "${cwd}" \
        --tpm-gct "${_input[0]}" \
        --counts-gct "${_input[1]}" \
        --annotation-gtf "${_input[2]}" \
        --sample-participant-lookup "${_input[3]}" \
        --tpm-threshold ${tpm_threshold} \
        --count-threshold ${count_threshold} \
        --sample-frac-threshold ${sample_frac_threshold} \
        --normalization-method "${normalization_method}" \
        ${"--quantile-normalize" if quantile_normalize else ""} \
        --numThreads ${numThreads}