Sample-level RNA-seq quality control

Sample-level RNA-seq quality control#

Sample-level quality control for a bulk expression matrix, applied before normalisation.

Overview#

Before normalisation it is worth asking whether every sample in the matrix belongs there. This step applies sample-level quality control to a bulk expression matrix - flagging samples whose overall expression profile is inconsistent with the rest of the cohort, which usually means a library preparation or labelling problem rather than biology.

Dropping such samples here matters more than it looks: they inflate the residual variance in the association scan and can distort the hidden-factor covariates estimated downstream.

Sample-level QC of a bulk RNA-seq expression matrix, following the GTEx V8 approach. It runs in two parts: qc_1 removes low-expression genes (genes expressed below the TPM threshold in too many samples), and qc_2 detects and removes outlier samples using three checks (a D-statistic histogram, a Relative Log Expression (RLE) plot, and hierarchical clustering). A raw count matrix, if provided, is filtered to the same genes and samples as the QC’d TPM matrix.

When to run it. After expression quantification, before normalisation and covariate preprocessing.

Input#

  • --tpm-gct input/rnaseq/protocol_example.rnaseq.tpm_matrix.bed (gene expression in TPM; gene ID in the first column, one column per sample. Required.)

gene_ID	SAMPLE_001	SAMPLE_002	SAMPLE_003	SAMPLE_004	SAMPLE_005	SAMPLE_006	SAMPLE_007	...
ENSG00000284070	4.46908614068004	3.6324384389436	0	0.926449274182614	8.44864560635242	16.3456237133014	4.07740959667351	...
ENSG00000273342	0	19.892726241898	13.6287618206549	28.5428919916226	10.1099209273507	17.0947587639983	0	...

  • --counts-gct input/rnaseq/protocol_example.rnaseq.count_matrix.bed (raw gene count matrix. Optional; if given it is filtered to match the QC’d TPM genes and samples.)

gene_ID	SAMPLE_001	SAMPLE_002	SAMPLE_003	SAMPLE_004	SAMPLE_005	SAMPLE_006	SAMPLE_007	...
ENSG00000284070	197	193	122	180	217	257	195	...
ENSG00000273342	80	274	243	318	226	260	133	...

  • --cwd output/rnaseq (working directory all outputs are written under. Defaults to output.)

  • --container (an empty string runs the workflow with the locally installed tools instead of a container)

The filtering thresholds all have defaults and only need setting to override them:

  • --low-expr-TPM 0.1 (a gene must exceed this TPM to count as expressed)

  • --low-expr-TPM-percent 0.2 (a gene is kept only if it is expressed in at least this fraction of samples)

  • --RLEFilterPercent 0.05 (fraction of samples dropped by the relative log expression filter)

  • --DSFilterPercent 0.05 (fraction of samples dropped by the D-statistic filter. For a small cohort raise it, e.g. 0.1, so a single sample is not over-flagged.)

  • --topk-genes 100 (number of top-expressed genes used for the sample clustering check)

  • --cluster-percent 0.6 (fraction of the top genes a sample must share to join a cluster)

  • --pvalue-cutoff 0.05 (significance threshold for the outlier test)

  • --cluster-level 5 (tree depth at which clusters are cut)

Output#

  • protocol_example.low_expression_filtered.tpm.gct.gz (qc_1) (TPM matrix after low-expression gene filtering. All outputs go to --cwd; the basename comes from the TPM input name truncated at the first dot, so here it is protocol_example.)

  • protocol_example.low_expression_filtered.outlier_removed.tpm.gct.gz (qc_2) (TPM matrix after outlier-sample removal)

  • protocol_example.low_expression_filtered.outlier_removed.geneCount.gct.gz (qc_3) (raw count matrix filtered to the same genes and samples)

Two diagnostic plots are also written but are not declared in any step output: statement, so SoS does not track them:

  • protocol_example.low_expression_filtered.outlier_removed.tpm.gct.RLEplot.pdf (relative log expression plot)

  • protocol_example.low_expression_filtered.outlier_removed.tpm.gct.D_stat_hist.pdf (D-statistic histogram; not present in the example output)

The qc_3 step additionally:

  1. Filter out the geneCount table based on TPM table.

  2. Adds two comment lines above the header of TPM and geneCount table to mimick the original output from RNASeQC.

Minimal Working Example#

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

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

Command Interface#

sos run pipeline/bulk_expression_QC.ipynb -h

usage: sos run pipeline/bulk_expression_QC.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:
  qc

Global Workflow Options:
  --modular-script-dir code/script (as path)
  --tpm-gct VAL (as path, required)
                        Required input is TPM file
  --counts-gct . (as path)
                        Raw counts file is optional and if available, it will be
                        filtered to match with the TPM file sample and genes
  --cwd output (as path)
  --container ''
  --entrypoint  ('micromamba run -a "" -n' + ' ' + re.sub(r'(_apptainer:latest|_docker:latest|\.sif)$', '', container.split('/')[-1])) if container else ""

  --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 8 (as int)
                        Number of threads

Sections
  qc_1:
    Workflow Options:
      --low-expr-TPM 0.1 (as float)
      --low-expr-TPM-percent 0.2 (as float)
  qc_2:
    Workflow Options:
      --RLEFilterPercent 0.05 (as float)
      --DSFilterPercent 0.05 (as float)
      --topk-genes 100 (as int)
      --cluster-percent 0.6 (as float)
      --pvalue-cutoff 0.05 (as float)
      --cluster-level 5 (as int)
  qc_3:

Workflow implementation#

[global]
parameter: modular_script_dir = path('code/script')  # override with --modular-script-dir
# Required input is TPM file
parameter: tpm_gct = path
# Raw counts file is optional and if available, it will be filtered to match with the TPM file sample and genes
parameter: counts_gct = path()
parameter: cwd = path("output")
parameter: container = ""
import re
parameter: entrypoint= ('micromamba run -a "" -n' + ' ' + re.sub(r'(_apptainer:latest|_docker:latest|\.sif)$', '', container.split('/')[-1])) if container else ""
cwd = path(f'{cwd:a}')
# 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 = 8

[qc_1 (basic check and low expression filtering)]
parameter: low_expr_TPM = 0.1
parameter: low_expr_TPM_percent = 0.2
input: tpm_gct
output: f'{cwd}/{_input:bnnn}.low_expression_filtered.tpm.gct.gz'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand= "${ }", stderr = f'{_output:nnn}.stderr', stdout = f'{_output:nnn}.log', container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/molecular_phenotypes/QC/bulk_expression_QC.R \
        --step qc_1 \
        --cwd "${cwd}" \
        --tpm-gct "${_input}" \
        --low-expr-TPM ${low_expr_TPM} \
        --low-expr-TPM-percent ${low_expr_TPM_percent} \
        --numThreads ${numThreads}

Caution. The offset used for the log transformation on TPM is unresolved. GTEx suggests an offset of 1; the recommendation kept here follows these authors, where both 0.0001 and 1 are used at different steps.

[qc_2 (remove outliers)]
parameter: RLEFilterPercent = 0.05
parameter: DSFilterPercent = 0.05
parameter: topk_genes = 100
parameter: cluster_percent = 0.6
parameter: pvalue_cutoff = 0.05
parameter: cluster_level = 5
output: f'{cwd}/{_input:bnnn}.outlier_removed.tpm.gct.gz'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand= "${ }", stderr = f'{_output:nnn}.stderr', stdout = f'{_output:nnn}.log', container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/molecular_phenotypes/QC/bulk_expression_QC.R \
        --step qc_2 \
        --cwd "${cwd}" \
        --tpm-gct "${_input}" \
        --RLEFilterPercent ${RLEFilterPercent} \
        --DSFilterPercent ${DSFilterPercent} \
        --topk-genes ${topk_genes} \
        --cluster-percent ${cluster_percent} \
        --pvalue-cutoff ${pvalue_cutoff} \
        --cluster-level ${cluster_level} \
        --numThreads ${numThreads}

Remove genes and samples from raw counts#

[qc_3 (remove gene and samples from raw counts)]
stop_if(not counts_gct.is_file())
output: f'{cwd}/{_input:bnnn}.geneCount.gct.gz'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand= "${ }", stderr = f'{_output:nn}.stderr', stdout = f'{_output:nn}.stdout', container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/molecular_phenotypes/QC/bulk_expression_QC.R \
        --step qc_3 \
        --cwd "${cwd}" \
        --tpm-gct "${_input}" \
        --counts-gct "${counts_gct}" \
        --numThreads ${numThreads}