Splicing QC and Normalization - FunGen-xQTL Consortium

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

Overview¶

Quality control and normalization of alternative-splicing quantifications produced by the leafcutter tool. Raw outputs are converted to BED format, features with high per-sample missingness (default 40%) are dropped, low-variance introns (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.

leafcutter_norm runs the whole route; leafcutter_qqnorm is the normalisation half on its own, for when QC and imputation happened elsewhere.

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 both workflows. Which form it takes depends on the workflow: leafcutter_norm expects the leafcutter per-cluster intron usage counts, protocol_example.leafcutter.intron_usage_perind.counts.gz (counts, count/total per sample), and leafcutter_qqnorm an already-QCed numeric intron-ratio table, ..._raw_data.txt.
--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
```

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 tests/fixtures/bulk_expression_normalization/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 tests/fixtures/bulk_expression_normalization/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 output/leafcutter/normalize/protocol_example.leafcutter.intron_usage_perind.counts.gz \
    --mean_impute

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:

QC - leafcutter_norm --no_norm stops after the raw-data table, writing {ratios}_raw_data.txt.
Imputation - phenotype_imputation.ipynb EBMF takes that table through --phenoFile (<output_from_step1>). This step runs a different module.
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/leafcutter/normalize \
    --ratios output/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/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 output/leafcutter/normalize/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/leafcutter/normalize \
    --ratios output/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

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)
      --qced-data . (as path)
      --[no-]bgzip (default to False)

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