Splicing QC and Normalization

Splicing QC and Normalization#

Quality control and normalisation of alternative-splicing quantifications from leafcutter and psichomics.

Overview#

Quality control and normalization of alternative-splicing quantifications produced by the leafcutter and psichomics tools. Raw outputs are converted to BED format, features with high per-sample missingness (default 40%) are dropped, low-variance introns/events (default min variance 0.001) are removed, values are quantile-normalized, and remaining missing values are mean-imputed. The normalized tables are ready to use as molecular phenotypes for tensorQTL.

Intron excision ratios are not comparable across samples as they come out of the caller. They are compositional - each intron is a fraction of its cluster - and sparse, because an intron unused in a sample has no reads rather than a zero. Left alone, both properties show up in association testing as signal that tracks sequencing depth and sample handling rather than genotype. The steps here filter clusters, impute the missing entries and quantile normalise, which is what makes the phenotype usable downstream.

Two quantification tools are supported and they need different handling, hence the separate workflows: leafcutter_norm and leafcutter_qqnorm for LeafCutter output, psichomics_norm for PSI values from psichomics.

When to run it. After splicing calling, before phenotype preprocessing and association testing.

Input#

  • --ratios: the splicing ratio matrix to normalise, and the main input to all three workflows. Which form it takes depends on the workflow: leafcutter_norm expects the leafcutter per-cluster intron usage counts, leafcutter_qqnorm an already-QCed numeric intron-ratio table, and psichomics_norm a psichomics PSI matrix keyed by splicing-event ID. Example input/rnaseq/protocol_example.psichomics.psi_raw_data.tsv:

    ID                                                    SAMPLE_001  SAMPLE_002  SAMPLE_003
SE_22_+_42438332_42438832_42439332_42439832_ENSG...   NA          0.3414      0.1858
SE_22_-_42440332_42440832_42441332_42441832_ENSG...   0.7942      0.1614      0.2509

    The leafcutter equivalents are protocol_example.leafcutter.intron_usage_perind.counts.gz (counts, count/total per sample) and ..._raw_data.txt (the QC-ready numeric table).

  • --qced-data: where step 3 looks for its input. A directory, the default ., is scanned for files ending _raw_data.txt. A single file is used directly, and that case also forces --mean_impute on to fill any remaining NAs.

  • --no_norm: stop leafcutter_norm before step 3, leaving the QCed table un-normalised. This is how the three-step route gets a QC-only result to hand to the imputer.

  • --mean_impute: off by default. Fills missing ratios by the mean; step 3 turns it on by itself when --qced-data names a file.

  • --cwd: the directory outputs are written to.

  • --container: the Apptainer or Docker image to run in, for example oras://ghcr.io/cumc/leafcutter_apptainer:latest.

  • --numThreads, --mem and --walltime: job resources.

  • --junc-path: the directory Junc_list globs for per-sample *.junc files. Required; it has no default.

  • --file-suffix: the extension Junc_list looks for, junc by default.

  • --dataset: names the junction list Junc_list writes, ROSMAP_DLPFC by default.

  • --sample-table: the sample lookup table Jointcall_samples reads, previewed above.

  • --junc-list: the junction list Jointcall_samples reads, normally the file Junc_list wrote.

  • --leafcutter-version: the subdirectory both steps write into, leafcutter2 by default.

Route B step 2 calls phenotype_imputation.ipynb, not this notebook, so the flags in that command (--phenoFile for the QCed table from step 1, and --prior and --varType for the EBMF model) are that module’s options and are documented there. They are not parameters of this notebook.

Output#

  • {cwd}/{leafcutter_version}/{dataset}_intron_usage_perind.junc - the per-cluster junction file assembled by Junc_list.

  • {cwd}/{leafcutter_version}/{sample_table}.rnaseq - the joint-called sample table from Jointcall_samples.

  • {ratios}_phenotype_file_list.txt - the phenotype file list written by the first leafcutter_norm step. Example input/rnaseq/protocol_example.leafcutter.intron_usage_perind.counts.gz_phenotype_file_list.stdout sits beside it.

  • {input}_raw_data.txt - the QC’d numeric intron-ratio table from the second leafcutter_norm step, which is also what leafcutter_qqnorm takes as its starting point.

  • {input}.qqnorm.txt - the quantile-normalised phenotype table, from either leafcutter_norm or the standalone leafcutter_qqnorm. Example input/rnaseq/protocol_example.leafcutter.intron_usage_perind.counts.gz_raw_data.qqnorm.txt, with a bgzipped and tabix-indexed .qqnorm.bed.gz beside it. Example input/rnaseq/protocol_example.leafcutter.intron_usage_perind.counts.gz_raw_data.qqnorm.txt, 24 columns and 112 rows, first 5 columns:

    #Chr  start  end  ID  SAMPLE_001
chr22  16671178  16671518  chr22:16671178:1667151  -0.25101591868327155
chr22  16671588  16673600  chr22:16671588:1667360  -1.261282157758418

  • {cwd}/psichomics_raw_data_bedded.txt and {input}.qqnorm.txt - the psichomics PSI matrix parsed into BED format and then quantile-normalised.

    Example output/splicing/psichomics_raw_data_bedded.txt, 24 columns and 31 rows:

    #Chr  start  end  ID
chr22  42438332  42439832  SE_22_+_42438332_424
chr22  42440332  42441832  SE_22_-_42440332_424

    and output/splicing/psichomics_raw_data_bedded.qqnorm.txt, 24 columns and 31 rows:

    #Chr  start  end  ID
chr22  42438332  42439832  SE_22_+_42438332_424
chr22  42440332  42441832  SE_22_-_42440332_424

Example output/splicing/leafcutter2/protocol_example_intron_usage_perind.junc from Junc_list, 2 lines, one path per sample:

output/splicing/junc/SAMPLE_001.junc
output/splicing/junc/SAMPLE_002.junc

Example output/splicing/leafcutter2/protocol_example.rnaseq.sample_participant_lookup.txt.rnaseq from Jointcall_samples, 3 lines:

sample_id	participant_id
SAMPLE_001	SAMPLE_001
SAMPLE_002	SAMPLE_002

Each workflow ends in a quantile-normalised phenotype table in BED-like form, #Chr start end ID plus one column per sample, ready to use directly as a molecular phenotype input to TensorQTL.

Minimal Working Example#

Build the junction list#

RNA-seq sample IDs differ from the WGS sample IDs, so the joint-call sample lookup table is first subset to the RNA-seq samples used in the sQTL analysis, matching the sample set used elsewhere in the protocol. The table pairs each sample ID with its participant ID. Example input/rnaseq/protocol_example.rnaseq.sample_participant_lookup.txt, 2 columns and 61 rows:

Junc_list globs --junc-path for per-sample *.junc files, resolves them to absolute paths and drops the all entry, leaving one list file. No directory of .junc files ships here, so --junc-path is a placeholder.

Timing: TBD (on toy dataset)

sos run pipeline/splicing_normalization.ipynb Junc_list \
    --cwd output/splicing \
    --junc-path <dir of per-sample *.junc files> \
    --dataset protocol_example \
    --leafcutter-version leafcutter2

Joint-call the samples#

Jointcall_samples takes the junction list from the previous step together with the sample lookup table and writes the .rnaseq table that the leafcutter routes consume. The two chain by hand: --junc-list is the file Junc_list wrote.

Timing: TBD (on toy dataset)

sos run pipeline/splicing_normalization.ipynb Jointcall_samples \
    --cwd output/splicing \
    --sample-table input/rnaseq/protocol_example.rnaseq.sample_participant_lookup.txt \
    --junc-list output/splicing/leafcutter2/protocol_example_intron_usage_perind.junc \
    --leafcutter-version leafcutter2

Leafcutter#

This workflow runs the full leafcutter chain in one command: it builds the phenotype table, applies an autosome filter and QC, and quantile-normalizes the intron-usage ratios. The --mean_impute option fills missing ratios with the per-intron (row) mean. See the leafcutter documentation; the choice of regtool parameters is discussed here.

Parameters. chr_blacklist is a file of blacklisted chromosomes to exclude from analysis, one per line; if none is provided, no chromosomes are excluded.

Timing: TBD (on toy dataset)

sos run pipeline/splicing_normalization.ipynb leafcutter_norm \
    --cwd output/splicing \
    --ratios input/rnaseq/protocol_example.leafcutter.intron_usage_perind.counts.gz \
    --mean_impute

Psichomics PSI values#

This workflow parses psichomics splicing-event IDs into BED coordinates, applies NA and variance filters, and then quantile-normalizes each event. See the psichomics documentation. The only QC applied to PSI values here is NA removal and a minimal variance filter; the psichomics maintainers suggest some further QC, discussed here. For reference, the default minimal variance used in leafcutter QC is 0.001. Because the alternative-splicing event types in psichomics outputs differ, normalization is performed per event (rows) rather than per sample (columns).

Timing: TBD (on toy dataset)

sos run pipeline/splicing_normalization.ipynb psichomics_norm \
    --cwd output/splicing \
    --ratios input/rnaseq/protocol_example.psichomics.psi_raw_data.tsv

Leafcutter with EBMF imputation#

This is the leafcutter route split into three commands so that EBMF imputation can be used instead of the default mean imputation. Run it instead of, not after, the single leafcutter_norm command above.

The steps are chained by hand rather than by SoS, so the placeholders in the commands below must be filled in with real paths:

  1. QC - leafcutter_norm --no_norm stops after the raw-data table, writing {ratios}_raw_data.txt.

  2. Imputation - phenotype_imputation.ipynb EBMF takes that table through --phenoFile (<output_from_step1>). This step runs a different module.

  3. Normalization - leafcutter_qqnorm takes the imputed matrix through --qced-data (<output_from_step2>). Because that argument names a file rather than a directory, this step also switches on mean imputation for any NAs EBMF left behind.

To perform leafcutter analysis with QC (setting cluster counts of 0 as NA), imputation (flashR), and normalization, use the following commands (RECOMMENDED):

Step 1. QC#

Timing: TBD (on toy dataset)

sos run pipeline/splicing_normalization.ipynb leafcutter_norm \
    --cwd ../../output_test/leafcutter/normalize \
    --ratios ../../output_test/leafcutter/PCC_sample_list_subset_leafcutter_intron_usage_perind.counts.gz \
    --container oras://ghcr.io/cumc/leafcutter_apptainer:latest \
    --no_norm # add no norm to skip last step (qqnorm) in leafcutter_norm

Step 2. Imputation#

Timing: TBD (on toy dataset)

sos run pipeline/phenotype_imputation.ipynb EBMF \
    --phenoFile <output_from_step1> \
    --cwd ../../output_test/leafcutter/imputation \
    --prior ebnm_point_laplace --varType 1 \
    --container oras://ghcr.io/cumc/factor_analysis_apptainer:latest \
    --mem 40G --numThreads 20 --walltime 100h

Step 3. Normalization#

Timing: TBD (on toy dataset)

sos run pipeline/splicing_normalization.ipynb leafcutter_qqnorm \
    --cwd output/splicing \
    --qced-data input/rnaseq/protocol_example.leafcutter.intron_usage_perind.counts.gz_raw_data.txt \
    --container oras://ghcr.io/cumc/leafcutter_apptainer:latest

Or if you are going to use the default mean imputation method, you can run it with one-step command (remove –no_norm in the first step):

Timing: TBD (on toy dataset)

sos run pipeline/splicing_normalization.ipynb leafcutter_norm \
    --cwd ../../output_test/leafcutter/normalize \
    --ratios ../../output_test/leafcutter/PCC_sample_list_subset_leafcutter_intron_usage_perind.counts.gz \
    --container oras://ghcr.io/cumc/leafcutter_apptainer:latest \
    --mean_impute

Command Interface#

sos run pipeline/splicing_normalization.ipynb -h

usage: sos run pipeline/splicing_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:
  Junc_list
  Jointcall_samples
  leafcutter_norm
  leafcutter_qqnorm
  psichomics_norm

Global Workflow Options:
  --modular-script-dir code/script (as path)
  --cwd output (as path)
                        The output directory for generated files.
  --chr-blacklist . (as path)
                        intron usage ratio file wiht samples after QC optional
                        parameter black list if user want to blacklist some
                        chromosomes and not to analyze
  --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
  --container ''
                        Software container option

Sections
  Junc_list:
    Workflow Options:
      --junc-path VAL (as path, required)
      --file-suffix junc
      --leafcutter-version leafcutter2
      --dataset 'ROSMAP_DLPFC'
  Jointcall_samples:
    Workflow Options:
      --sample-table VAL (as path, required)
      --junc-list VAL (as path, required)
      --leafcutter-version leafcutter2
  leafcutter_norm_1:
    Workflow Options:
      --ratios VAL (as path, required)
      --[no-]pseudo-ratio (default to False)
  leafcutter_norm_2:
    Workflow Options:
      --[no-]autosomes (default to True)
  leafcutter_norm_3, leafcutter_qqnorm:
    Workflow Options:
      --[no-]no-norm (default to False)
      --[no-]mean-impute (default to False)
      --na-rate 0.4 (as float)
                        minimal NA rate within sample values for an event to be
                        kept (default 0.4, leafcutter default)
      --min-variance 0.001 (as float)
                        minimal variance across samples for an event to be kept
                        (default 0.001, psichomics suggested)
      --qced-data . (as path)
      --[no-]bgzip (default to False)
  psichomics_norm_1:
    Workflow Options:
      --ratios VAL (as path, required)
  psichomics_norm_2:
    Workflow Options:
      --na-rate 0.4 (as float)
                        minimal NA rate within sample values for an event to be
                        kept (default 0.4, leafcutter default)
      --min-variance 0.001 (as float)
                        minimal variance across samples for an event to be kept
                        (default 0.001, psichomics suggested)

Workflow implementation#

[global]
parameter: modular_script_dir = path('code/script')  # override with --modular-script-dir
# The output directory for generated files. 
parameter: cwd = path("output")
# intron usage ratio file wiht samples after QC
# optional parameter black list if user want to blacklist some chromosomes and not to analyze
parameter: chr_blacklist = path(".")
# 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
# Software container option
parameter: container = ""
from sos.utils import expand_size
cwd = path(f'{cwd:a}')

[Junc_list]
parameter: junc_path = path
parameter: file_suffix = "junc"
parameter: leafcutter_version = "leafcutter2"
parameter: dataset = "ROSMAP_DLPFC"
output: f'{cwd}/{leafcutter_version}/{dataset}_intron_usage_perind.junc'
bash: container=container, expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stderr' 
    cd ${cwd}
    realpath  ${junc_path}/*${file_suffix} > ${_output}.tmp
    grep -v "all" ${_output}.tmp > ${_output}
    rm ${_output}.tmp

[Jointcall_samples]
parameter: sample_table = path
parameter: junc_list = path
parameter: leafcutter_version = "leafcutter2"
input: sample_table, junc_list
output: f'{cwd}/{leafcutter_version}/{sample_table:b}.rnaseq'
bash: container=container, expand = "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stderr'
    Rscript ${modular_script_dir}/molecular_phenotypes/QC/splicing_normalization.R \
        --step jointcall_samples \
        --sample-table ${_input[0]} \
        --junc-list ${_input[1]} \
        --output ${_output}

[leafcutter_norm_1]
parameter: ratios = path
parameter: pseudo_ratio = False 
import os
if os.path.isfile(f'{ratios:dd}/black_list.txt'):
    chr_blacklist = f'{ratios:dd}/black_list.txt'
input: ratios, group_by = 'all'
output: f'{ratios}_phenotype_file_list.txt'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout', container=container
    Rscript ${modular_script_dir}/molecular_phenotypes/QC/splicing_normalization.R \
        --step prepare_phenotype \
        --input ${ratios} \
        --chr-blacklist ${chr_blacklist} \
        ${('--pseudo-ratio' if pseudo_ratio else '')}

[leafcutter_norm_2]
parameter: autosomes = True
# Build the per-chromosome file list for concatenation (no pandas; runs in the base env).
_lines = open(f'{_input[0]}').read().splitlines()
_hdr = _lines[0].split('\t')
_rows = [dict(zip(_hdr, _ln.split('\t'))) for _ln in _lines[1:] if _ln.strip()]
if autosomes:
    print("Analyzing autosomes 1 to 22...")
    _rows = [_r for _r in _rows if _r['#chr'].startswith('chr') and _r['#chr'][3:].isdigit() and 1 <= int(_r['#chr'][3:]) <= 22]
file_inv = [_r['#dir'] for _r in _rows]
del _lines, _hdr, _rows
input: file_inv # This design is necessary to avoid using for_each, as sos can not take chr number as an input.
output: f'{_input[0]:n}_raw_data.txt'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container=container
    head -1 ${_input[0]:r}  > ${_output}
    cat ${_input:r} | grep -v "#Chr" >> ${_output}

[leafcutter_norm_3,leafcutter_qqnorm]
parameter: no_norm = False
stop_if(no_norm)
parameter: mean_impute = False
# minimal NA rate within sample values for an event to be kept (default 0.4, leafcutter default)
parameter: na_rate = 0.4
# minimal variance across samples for an event to be kept (default 0.001, psichomics suggested)
parameter: min_variance = 0.001
parameter: qced_data = path(".")
parameter: bgzip=False
import os
if qced_data.is_dir():
    print('Target Qced file with suffix "_raw_data.txt"')
    input_file = paths([os.path.join(qced_data, x) for x in os.listdir(qced_data) if x.endswith('_raw_data.txt')])
    bgzip=True
else:
    print('Loading QCed data as input file; enabling mean imputation for any remaining NAs.')
    input_file = paths(qced_data)
    mean_impute = True
input: input_file
output: f'{_input:n}.qqnorm.txt'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout', container=container
    Rscript ${modular_script_dir}/molecular_phenotypes/QC/splicing_normalization.R \
        --step qqnorm \
        --input ${_input} \
        --output ${_output} \
        --na-rate ${na_rate} \
        --min-variance ${min_variance} \
        ${('--mean-impute' if mean_impute else '')}
    if ${'true' if bgzip else 'false'}; then
        bgzip -c "${_output}" > "${_output:n}.bed.gz"
        tabix -p bed "${_output:n}.bed.gz"
    fi

[psichomics_norm_1]
parameter: ratios = path
input: ratios, group_by = 'all'
output: f'{cwd}/psichomics_raw_data_bedded.txt' 
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', container=container
    Rscript ${modular_script_dir}/molecular_phenotypes/QC/splicing_normalization.R \
        --step psichomics_bed \
        --input ${_input} \
        --output ${_output}

[psichomics_norm_2]
# minimal NA rate within sample values for an event to be kept (default 0.4, leafcutter default)
parameter: na_rate = 0.4
# minimal variance across samples for an event to be kept (default 0.001, psichomics suggested)
parameter: min_variance = 0.001
output: f'{_input:n}.qqnorm.txt'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout', container=container
    Rscript ${modular_script_dir}/molecular_phenotypes/QC/splicing_normalization.R \
        --step qqnorm \
        --input ${_input} \
        --output ${_output} \
        --na-rate ${na_rate} \
        --min-variance ${min_variance}