Quantifies alternative splicing from aligned RNA-seq with LeafCutter, producing the intron-usage 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.
Splicing is quantified as alternatively excised intron usage, in the sense described by Wang et al. (2008) and Park et al. (2018).
leafcutterquantifies 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.leafcutter_preprocessingis not a separate method but the first half ofleafcutter. The two share one section, declared as[leafcutter_1, leafcutter_preprocessing_1], and this workflow stops once the list of per-sample.juncfiles is written, which is exactly the fileleafcutterwould go on to cluster. Run it only when the clustering is done elsewhere.
When to run it. After STAR alignment, and before splicing_normalization, which turns the intron-usage 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 byleafcutter, and theSJ.out.tabpath. Strandness isrf,fr, orstrand_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.baiindexes already sit beside the BAMs, LeafCutter reuses them rather than re-indexing.--cwd-- output directory, defaultoutput.--container-- the tool image,oras://ghcr.io/statfungen/leafcutter_apptainer:latest.
Junction extraction (leafcutter, leafcutter_preprocessing):
--anchor-len-- minimum anchor length on each side of a junction, default8.--min-intron-len-- shortest intron considered, default50.--max-intron-len-- longest intron considered, default500000.
Intron clustering (leafcutter):
--min-clu-reads-- minimum reads in a cluster, default30.--min-clu-ratio-- minimum fraction of cluster reads supporting a junction, default0.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, a manifest and a LeafCutter blacklist-chromosome file, are on Google Drive.
Runtime:
--numThreads,--job-size,--walltime,--mem-- threads (default8) and cluster resources.--modular-script-dir-- location of the shell and R drivers, defaultcode/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, keyedchrom: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>_pooledand<samples>_refined-- intermediate cluster definitions written during clustering.<samples>_intron_usage_perind.junc-- the list of per-sample.juncfiles. It is the final output ofleafcutter_preprocessing, andleafcutterwrites 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.
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.
The intron-usage table feeds 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:latestCommand Interface¶
sos run pipeline/splicing_calling.ipynb -husage: 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
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: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']
del _lines, _hdr, _rows, _smLeafCutter: 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- Wang, E. T., Sandberg, R., Luo, S., Khrebtukova, I., Zhang, L., Mayr, C., Kingsmore, S. F., Schroth, G. P., & Burge, C. B. (2008). Alternative isoform regulation in human tissue transcriptomes. Nature, 456(7221), 470–476. 10.1038/nature07509
- Park, E., Pan, Z., Zhang, Z., Lin, L., & Xing, Y. (2018). The Expanding Landscape of Alternative Splicing Variation in Human Populations. The American Journal of Human Genetics, 102(1), 11–26. 10.1016/j.ajhg.2017.11.002
- Li, Y. I., Knowles, D. A., Humphrey, J., Barbeira, A. N., Dickinson, S. P., Im, H. K., & Pritchard, J. K. (2017). Annotation-free quantification of RNA splicing using LeafCutter. Nature Genetics, 50(1), 151–158. 10.1038/s41588-017-0004-9