Quantifying alternative splicing from RNA-seq data

Quantifying alternative splicing from RNA-seq data#

Quantifies alternative splicing from aligned RNA-seq with LeafCutter and psichomics, producing the intron-usage and PSI tables that become splicing phenotypes.

Overview#

This module turns aligned RNA-seq into splicing phenotypes. It expects BAM files already mapped with STAR using the WASP option, following the GTEx/TOPMed RNA-seq pipeline; the choice of modules is supported by internal, unpublished benchmarks from the GTEx group.

Two tools are provided because they measure different things, in the sense described by Wang et al. (2008) and Park et al. (2018). They do not depend on each other and can run in parallel.

leafcutter quantifies the usage of alternatively excised introns. Junctions are extracted from each BAM with regtools, then clustered across samples: introns sharing a splice site are grouped, and each intron is reported as a fraction of the reads in its cluster. That single measure collectively captures skipped exons, alternative 5-prime and 3-prime splice site usage, and more complex events, without needing to name the event type Li et al. (2018). The clustering defaults, --min-clu-ratio 0.001 --max-intron-len 500000 --min-clu-reads 30, are those of the GTEx sQTL discovery pipeline, section 3.4.3. This is the approach previously applied to ROSMAP data for Brain xQTL version 2.0.
psichomics quantifies named event types instead, reporting percent-spliced-in (PSI) for skipped exons, mutually exclusive exons, alternative 5-prime and 3-prime splice sites, and alternative first and last exons Agostinho et al. (2019). It reads the SJ.out.tab junction counts STAR emits rather than the BAMs, and needs a splicing annotation to say which junctions constitute which event.
leafcutter_preprocessing is not a third method but the first half of leafcutter. The two share one section, declared as [leafcutter_1, leafcutter_preprocessing_1], and this workflow stops once the list of per-sample .junc files is written, which is exactly the file leafcutter would go on to cluster. Run it only when the clustering is done elsewhere.

Two psichomics annotations are distributed with the example data. hg38_suppa.rds was built with SUPPA from the GTF used throughout this protocol; modified_psichomics_hg38_splicing_annotation.rds is the package default Human hg38 (2018-04-30) annotation, described in the psichomics methods, with gene names replaced by Ensembl IDs so that identifiers are consistent across the protocol.

When to run it. After STAR alignment, and before splicing_normalization, which turns the intron-usage and PSI tables produced here into the splicing phenotype matrices used for sQTL association.

Input#

--samples – the sample manifest, white-space delimited with a header and four columns: sample ID, RNA strandness, the BAM path used by leafcutter, and the SJ.out.tab path used by psichomics. Strandness is rf, fr, or strand_missing:

sample_id strand bam_list SJ_list
sample_1 rf sample_1.Aligned.sortedByCoord.out.bam sample_1.SJ.out.tab
sample_2 fr sample_2.Aligned.sortedByCoord.out.bam sample_2.SJ.out.tab
sample_3 strand_missing sample_3.Aligned.sortedByCoord.out.bam sample_3.SJ.out.tab

--data-dir – directory holding the files named in the manifest; defaults to the directory containing --samples. Every file listed must be found there. If .bam.bai indexes already sit beside the BAMs, LeafCutter reuses them rather than re-indexing.
--splicing-annotation – required by psichomics: the annotation R object naming which junction combinations form which event, e.g. Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.SUPPA_annotation.rds.
--cwd – output directory, default output.
--container – the tool images, oras://ghcr.io/statfungen/leafcutter_apptainer:latest and oras://ghcr.io/statfungen/psichomics_apptainer:latest. The two workflows need different containers, so they are run as separate commands.

Junction extraction (leafcutter, leafcutter_preprocessing):

--anchor-len – minimum anchor length on each side of a junction, default 8.
--min-intron-len – shortest intron considered, default 50.
--max-intron-len – longest intron considered, default 500000.

Intron clustering (leafcutter):

--min-clu-reads – minimum reads in a cluster, default 30.
--min-clu-ratio – minimum fraction of cluster reads supporting a junction, default 0.001.

Reference data and example input for this module are not in the repository. The reference data is built by reference_data_preparation, and the example inputs, both psichomics annotations, a manifest and a LeafCutter blacklist-chromosome file, are on Google Drive.

Runtime:

--numThreads, --job-size, --walltime, --mem – threads (default 8) and cluster resources.
--modular-script-dir – location of the shell and R drivers, default code/script.

Output#

<sample>.junc – per-sample junction files extracted from the BAMs by regtools, written by both LeafCutter workflows.
<samples>_intron_usage_perind.counts.gz – the LeafCutter splicing phenotype: one row per intron, keyed chrom:start:end:clu_N, giving the fraction of that cluster reads supporting the intron in each sample.
<samples>_intron_usage_perind_numers.counts.gz – the same rows as raw supporting-read counts rather than ratios.
<samples>_pooled and <samples>_refined – intermediate cluster definitions written during clustering.
<samples>_intron_usage_perind.junc – the list of per-sample .junc files. It is the final output of leafcutter_preprocessing, and leafcutter writes the same file on its way to clustering, although it declares only the counts file as its output, so SoS does not track this one.
psichomics_junctions.txt – the junction quantification table. When the manifest spans several input directories, one psichomics_junctions_<directory>.txt is written per directory and then merged into this file.
psi_raw_data.tsv – percent-spliced-in values per event and sample, the psichomics splicing phenotype. Mitochondrial, unplaced and random-contig junctions are dropped before quantification.

Each step also writes .stdout and .stderr beside its output. There is no example output for this module in the repository, because its inputs are aligned BAM files that are not distributed here.

Both the intron-usage and the PSI table feed splicing_normalization, which builds the final splicing phenotype matrices.

Minimal Working Example#

LeafCutter intron usage#

Extract junctions from every BAM and cluster them into intron-excision ratios. Around 30 minutes on a full sample set.

Timing: TBD (on toy dataset)

sos run pipeline/splicing_calling.ipynb leafcutter \
    --cwd output/leafcutter \
    --samples path/to/sample_list.txt \
    --data-dir path/to/star_output_wasp \
    --container oras://ghcr.io/statfungen/leafcutter_apptainer:latest

LeafCutter junction extraction only#

The first half of the leafcutter command above, stopping once the list of .junc files is written. leafcutter writes that same list itself, so this workflow is worth running on its own only when the clustering happens elsewhere.

Timing: TBD (on toy dataset)

sos run pipeline/splicing_calling.ipynb leafcutter_preprocessing \
    --cwd output/leafcutter \
    --samples path/to/sample_list.txt \
    --data-dir path/to/star_output_wasp \
    --container oras://ghcr.io/statfungen/leafcutter_apptainer:latest

psichomics PSI quantification#

Convert the STAR junction counts into a psichomics junction table and quantify percent-spliced-in against the splicing annotation. Around 30 minutes on a full sample set.

Timing: TBD (on toy dataset)

sos run pipeline/splicing_calling.ipynb psichomics \
    --cwd output/psichomics \
    --samples path/to/sample_list.txt \
    --data-dir path/to/star_output_wasp \
    --splicing-annotation path/to/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.SUPPA_annotation.rds \
    --container oras://ghcr.io/statfungen/psichomics_apptainer:latest

Command Interface#

sos run pipeline/splicing_calling.ipynb -h

usage: sos run pipeline/splicing_calling.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:
  leafcutter
  leafcutter_preprocessing
  psichomics

Global Workflow Options:
  --modular-script-dir code/script (as path)
  --cwd output (as path)
                        The output directory for generated files.
  --samples VAL (as path, required)
                        Sample meta data list
  --data-dir  path(f"{samples:d}")

                        Raw data directory, default to the same directory as
                        sample list
  --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
  leafcutter_1, leafcutter_preprocessing_1:
    Workflow Options:
      --anchor-len 8 (as int)
                        anchor length (default 8)
      --min-intron-len 50 (as int)
                        minimum intron length to be analyzed (default 50)
      --max-intron-len 500000 (as int)
                        maximum intron length to be analyzed (default 500000)
  leafcutter_2:
    Workflow Options:
      --min-clu-reads 30 (as int)
                        minimum reads in a cluster (default 50 reads)
      --max-intron-len 500000 (as int)
                        maximum intron length to be analyzed (default 500000)
      --min-clu-ratio 0.001 (as float)
                        minimum fraction of reads in a cluster that support a
                        junction (default 0.001)
  leafcutter_preprocessing_2:
  psichomics_1:
  psichomics_2:
    Workflow Options:
      --splicing-annotation VAL (as path, required)
                        splicing annotation for psichomics (raw psichomics
                        annotation, e.g. loadAnnotation output)

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")
# Sample meta data list
parameter: samples = path
# Raw data directory, default to the same directory as sample list
parameter: data_dir = path(f"{samples:d}")

# 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}')

# Build the input file lists from the sample manifest (basic SoS input plumbing; no pandas so
# it runs in the base env, no def/handles left in globals so the lists pickle to task workers).
# regtools strandness codes: rf->RF, fr->FR, strand_missing->XS. Blank/"NA" drops that sample.
_lines = open(samples).read().splitlines()
_hdr = _lines[0].split('\t')
_rows = [dict(zip(_hdr, _ln.split('\t'))) for _ln in _lines[1:] if _ln.strip()]
_sm = {'rf': 'RF', 'fr': 'FR', 'strand_missing': 'XS'}
sample_id  = [(_r.get('sample_id') or 'NA') for _r in _rows]
strandness = [_sm.get((_r.get('strand') or 'NA'), (_r.get('strand') or 'NA')) for _r in _rows]
bam_data   = [f'{data_dir}/{_r.get("coord_bam_list")}' for _r in _rows if (_r.get('coord_bam_list') or 'NA') != 'NA']
SJtab_data = [f'{data_dir}/{_r.get("SJ_list")}' for _r in _rows if (_r.get('SJ_list') or 'NA') != 'NA']
del _lines, _hdr, _rows, _sm

LeafCutter: junction extraction and clustering#

Documentation: LeafCutter; the regtools parameter choices are discussed here. leafcutter_1 extracts junctions from each BAM, leafcutter_2 clusters them into intron-usage ratios, and leafcutter_preprocessing_2 writes the junction list instead of clustering. The underlying clustering script has further options this module does not expose, among them reusing an existing cluster file, skipping the chromosome-name check, and including constitutive introns.

[leafcutter_1, leafcutter_preprocessing_1]
# anchor length (default 8)
parameter: anchor_len = 8
# minimum intron length to be analyzed (default 50)
parameter: min_intron_len = 50
# maximum intron length to be analyzed (default 500000)
parameter: max_intron_len = 500000
input: bam_data, group_by = 1, group_with = "strandness"
output: f'{cwd}/{_input:bn}.junc' 
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
    bash ${modular_script_dir}/molecular_phenotypes/calling/regtools_junctions.sh \
        --bam ${_input} \
        --output ${_output} \
        --min-anchor ${anchor_len} \
        --min-intron ${min_intron_len} \
        --max-intron ${max_intron_len} \
        --strandness ${_strandness}

[leafcutter_2]
# minimum reads in a cluster (default 50 reads)
parameter: min_clu_reads = 30 
# maximum intron length to be analyzed (default 500000)
parameter: max_intron_len = 500000 
# minimum fraction of reads in a cluster that support a junction (default 0.001)
parameter: min_clu_ratio = 0.001
input: group_by = 'all'
output: f'{cwd}/{samples:bn}_intron_usage_perind.counts.gz'
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
    rm -f ${_output:nn}.junc
    for i in ${_input:r}; do
    echo $i >> ${_output:nn}.junc ; done
    Rscript ${modular_script_dir}/molecular_phenotypes/calling/leafcutter_cluster_regtools.R \
        --juncfiles ${_output:nn}.junc \
        --outprefix ${f'{_output:bnn}'.replace("_perind","")} \
        --rundir ${cwd} \
        --minclureads ${min_clu_reads} \
        --maxintronlen ${max_intron_len} \
        --mincluratio ${min_clu_ratio}

[leafcutter_preprocessing_2]
input: group_by = 'all'
output: f'{cwd}/{samples:bn}_intron_usage_perind.junc'
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
    rm -f ${_output:r}
    for i in ${_input:r}; do
    echo $i >> ${_output:r} ; done

psichomics: junction quantification and PSI#

Documentation: psichomics. psichomics_1 turns the STAR junction counts into a psichomics junction table and psichomics_2 quantifies PSI against the annotation. The underlying quantifySplicing call also accepts a restriction to particular event types, a minReads floor below which a value is returned as NA (default 10), and a gene subset; this module quantifies every event type across all genes.

[psichomics_1]
input: SJtab_data
output: f'{cwd}/psichomics_junctions.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/calling/splicing_calling.R \
        --step junction_quant \
        --inputs ${_input} \
        --output ${_output}

[psichomics_2]
# splicing annotation for psichomics (raw psichomics annotation, e.g. loadAnnotation output)
parameter: splicing_annotation = path
output: f'{cwd}/psi_raw_data.tsv'
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/calling/splicing_calling.R \
        --step quantify_psi \
        --junctions ${cwd}/psichomics_junctions.txt \
        --splicing-annotation ${splicing_annotation:a} \
        --output ${_output}

Quantifying alternative splicing from RNA-seq data

Contents

Quantifying alternative splicing from RNA-seq data#

Overview#

Input#

Output#

Minimal Working Example#

LeafCutter intron usage#

LeafCutter junction extraction only#

psichomics PSI quantification#

Command Interface#

Workflow implementation#

LeafCutter: junction extraction and clustering#

psichomics: junction quantification and PSI#