junctionCounts

junctionCounts is a tool intended to identify and quantify binary alternative splicing events from RNA-seq data. junctionCounts currently consists of two primary utilities: infer_pairwise_events.py, which identifies binary alternative splicing events (e.g. skipped exon events), and junctionCounts.py, which quantifies those events (i.e. provides a percent-spliced-in, or PSI value) from RNA-seq data (specifically from a BAM file). These utilities can be used independent of one another, and events identified independently by infer_pairwise_events.py and the transcripts with which they are associated (as marked in the .ioe files - a format introduced by the authors of SUPPA) may be useful to a variety of researchers, regardless of their interest in performing short read quantification. While junctionCounts.py could in principle be run without the use of infer_pairwise_events.py, it may be easiest to use both tools together for quantification due to the specific reqirements of the input GTF for the quantification. For predicting the functional outcomes of alternative splicing events, including NMD, NSD and coding-to-noncoding switch events, check out junctionCounts' partner utilities at https://github.com/ajw2329/cds_insertion.

Dependencies

python3
ncls (https://github.com/hunt-genes/ncls)
numpy
pysam splice_lib (github.com/ajw2329/splice_lib/ - install using python setup.py install)
bedtools
R (≥ 3.5.0) dplyr DEXSeq pacman optparse

Note that ncls, numpy, and pysam are all python packages that can be installed using pip install package_name.

Benchmarking

Panels A-K: We evaluated junctionCounts on its performance relative to four established splicing analysis tools: MAJIQ v2.5.6.dev1+g8423f68b (https://doi.org/10.1038/s41467-023-36585-y), rMATS-turbo v4.3.0 (https://doi.org/10.1038/s41596-023-00944-2), splAdder v3.0.4 (https://doi.org/10.1093/bioinformatics/btw076) and Whippet v1.6.2 (https://doi.org/10.1016/j.molcel.2018.08.018). We generated four paired-end simulated datasets from real RNA-seq data using polyester v1.36.0 (https://doi.org/10.1093/bioinformatics/btv272) with the arguments: read.length = 100, fragment.length.min = 100, fragment.length.max = 500, fragment.length.mean = 180, fragment.length.sd = 40 and simulated.sequencing.error = TRUE. Three of the simulated datasets: 25M, 50M and 75M, were used to evaluate performance at different library depths (25, 50 and 75 million reads per library respectively). These datasets were modeled on mouse cerebellum and liver RNA-seq data in triplicate for each cell type from Vaquero-Garcia et al. (2016) (https://doi.org/10.7554/eLife.11752), who replicated the experiments in Zhang et al. (2014) (https://doi.org/10.1073/pnas.1408886111). The fourth simulated dataset was modeled on human RNA-seq data from HeLa cells treated with either spliceostatin A or DMSO (https://doi.org/10.1038/s41467-017-02007-z) at 50 million reads per library with triplicates for each condition. This dataset introduced a larger pool of potential AS events to test compared to the murine simulated datasets, and importantly increased the total number of MR events tested.

Panels L-N: To test junctionCounts' NMD event predictions, we procured RT-PCR data for well-characterized NMD targets and accompanying RNA-seq data from HEK-293 cells upon UPF1 knockdown versus non-targeting siRNA (https://doi.org/10.15252/embj.2021108898). We then ran junctionCounts and its partner utilities with default settings on the RNA-seq data to predict NMD switch events within the NMD targets verified with RT-PCR in the original study. Out of the 13 predicted NMD switch events associated with the NMD targets, 11 had the expected dPSI directionality (Panel L). We then used the same approach on a dataset in which treatment with emetine, a translation elongation inhibitor, was reported to increase the abundance of NMD substrates in HEK-293 T cells (https://doi.org/10.1093/nar/gkw1109). cdsInsertion and findSwitchEvents identified 636 potential NMD switch events with significant changes in splicing (|dPSI| ≥ 0.1 and Q-value ≤ 0.05) upon emetine treatment, out of which 472 (74.2%) exhibited the expected dPSI directionality (Panels M, N).

infer_pairwise_events.py

Description

infer_pairwise_events.py takes as input a transcriptome GTF file and performs a pairwise comparison of all overlapping transcripts to identify the set of (one or more) minimal binary alternative splicing (or other forms of co/post-transcriptional processing) events that differentiate each pair. infer_pairwise_events.py considers an alternative event to be a collection of non-shared (i.e. unique to one transcript in the considered pair) exon boundaries that is bounded either a) by a set of shared outer exon boundaries, or b) by transcript termini. This definition includes all canonical binary alternative splicing events (of which the author is aware such as skipped exons, alternative donors/acceptors, alternative first/last exons) while also including non-canonical alternative splicing events. A schematic illustrating the events in question is below. Note that certain event "types" (specifically CO, CF, and CL) are generic terms that stand for 'complex', 'complex first', and 'complex last' respectively. These terms do not relate to modern concepts of alternative splicing complexity (i.e. non-binary events) but rather are buckets for events that both a) do not fit into the other, more specific categories, and b) do not involve the transcript termini (CO), involve different 5'-termini (CF), or involve different 3'-termini (CL). This also applies to event types MF, ML, MS, MR (multiple first, last, multiple skipped exon, and multiple retained intron respectively) - these events are highly analogous to AF, AL and SE events, but are not restricted to one pair of alternative exons in the case of AF, AL (for example, an MF or ML event could have 3 exons unique to the included form and two unique to the excluded form) or a single alternative exon in the case of SE.

The meaning of each abbreviation (bearing in mind the above discussion) is as follows:
SE - skipped exon
MS - multiple skipped exon
MX - mutually exclusive exons
A5 - alternative 5'-splice site
A3 - alternative 3'-splice site
RI - retained intron
MR - multiple retained intron
AF - alternative first exon
AL - alternative last exon
MF - multiple first exon
ML - multiple last exon
CF - complex first
CL - complex last

Basic usage

python /path/to/infer_pairwise_events.py --transcript_gtf /path/to/transcriptome.gtf --outdir /path/to/output/dir/

Input file description

transcriptome GTF

Output file descriptions

splice_lib_events.gtf

splice_lib_events.bed

splice_lib_events.ioe

splice_lib_events.gff3

junctioncounts_unfriendly_events.gtf

junctioncounts_unfriendly_events.bed

junctioncounts_unfriendly_events.ioe

junctioncounts_unfriendly_events.gff3

Help statement

usage: infer_pairwise_events.py [-h] [--transcript_gtf TRANSCRIPT_GTF]
                                --outdir OUTDIR
                                [--min_exon_length MIN_EXON_LENGTH]
                                [--min_intron_length MIN_INTRON_LENGTH]
                                [--max_exon_length MAX_EXON_LENGTH]
                                [--max_intron_length MAX_INTRON_LENGTH]
                                [--dump_pkl_file]
                                [--bedtools_path BEDTOOLS_PATH]
                                [--suppress_output]
                                [--min_AP_AT_dist MIN_AP_AT_DIST]

optional arguments:
  -h, --help            show this help message and exit
  --transcript_gtf TRANSCRIPT_GTF
                        Full transcript gtf file. Not required, but if not
                        provided a transcript dict must be passed as a
                        parameter to the main function.
  --outdir OUTDIR       Path to output directory
  --min_exon_length MIN_EXON_LENGTH
                        Minimum allowable exon length in input gtf.
                        Transcripts with shorter exons will be filtered.
                        (default: 3)
  --min_intron_length MIN_INTRON_LENGTH
                        Minimum allowable intron length in input gtf.
                        Transcripts with shorter introns will be filtered.
                        (default: 20)
  --max_exon_length MAX_EXON_LENGTH
                        Maximum allowable exon length in input gtf.
                        Transcripts with longer exons will be filtered
                        (default = 35000)
  --max_intron_length MAX_INTRON_LENGTH
                        Maximum allowable intron length in input gtf.
                        Transcripts with longer introns will be fitlered
                        (default = 1000000)
  --dump_pkl_file       If set, program will dump pickle file of event dict.
  --bedtools_path BEDTOOLS_PATH
                        Path to bedtools executable (default = 'bedtools')
  --suppress_output     If set, GTF, GFF3, and IOE files will not be written.
  --min_AP_AT_dist MIN_AP_AT_DIST
                        Specifies minimum distance between alternative
                        polyadenylation sites, TSS in order for AP, AT events
                        to be retained. Default is 100

junctionCounts.py

Description

Basic usage

python /path/to/junctionCounts.py --event_gtf /path/to/splice_lib_events.gtf --bam /path/to/input_file.bam --outdir /path/to/output/dir/ --sample_name sample_name

Output file descriptions

h9esc_cyto_jd004_count_psi_outfile.tsv

Field descriptions:

sample_name : sample name provided by user via --sample_name
event_id : unique name for event
event_type : type of event as discussed above
min_ijc : (integer) minimum read count of all junctions in the included form
min_ejc : (integer) minimum read count of all junctions in the excluded (skipped) form
avg_psi : (float) average of all PSI values computed all pairs of included, excluded junction counts
max_gene_frac : (float) max(min_ijc, min_sjc)/max_gene_jc where max_gene_jc is the maximum junction count for any junction in the event's gene. Only computed if --calc_gene_frac is passed along with an ioe file
all_ijc : (comma-sep list of integers) list of junction counts for all included form junctions
all_ejc : (comma-sep list of integers) list of junction counts for all excluded form junctions
avg_ijc : (float) mean junction count for included form junctions
avg_ejc : (float) mean junction count for excluded form junctions
span_psi : (float) max_psi - min_psi, i.e. the difference between the largest and smallest PSI values computed using any pair of included, excluded form junction counts. Large spans may indicate unreliability of the PSI estimate for a number of reasons such as partial junction overlap with another event combined with an abundance of multiply counted reads.
min_psi : (float) minimum of all PSI values computed all pairs of included, excluded junction counts
ijc_min_psi : (float) included form junction count corresponding to the min_psi calculation
ejc_min_psi : (float) excluded form junction count corresponding to the min_psi calculation
max_psi : (float) maximum of all PSI values computed all pairs of included, excluded junction counts
ijc_max_psi : (float) included form junction count corresponding to the max_psi calculation
ejc_max_psi : (float) excluded form junction count corresponding to the max_psi calculation
mid_psi : (float) midpoint between min_psi and max_psi

Help statement

usage: junctionCounts.py [-h] [--event_gtf EVENT_GTF] --bam BAM [--se]
                         [--forward_read FORWARD_READ] --outdir OUTDIR
                         --sample_name SAMPLE_NAME [--dump_pkl_dict]
                         [--gzipped] [--event_ioe EVENT_IOE]
                         [--calc_gene_frac] [--suppress_output]
                         [--enable_edge_use] [--turn_off_no_ends]
                         [--suppress_eij_use] [--disable_ri_extrapolation]
                         [--n_bootstraps N_BOOTSTRAPS]

optional arguments:
  -h, --help            show this help message and exit
  --event_gtf EVENT_GTF
                        GTF describing pairwise events
  --bam BAM             BAM read file for counting junction reads
  --se                  BAM file is single-ended. Default assumes paired-end
                        reads.
  --forward_read FORWARD_READ
                        Specify read that matches the 'forward' strand.
                        Options are 'R1','R2' or 'unstranded' if the library
                        is unstranded. Unstranded use is not currently
                        recommended. default = 'R2'
  --outdir OUTDIR       Path for output files
  --sample_name SAMPLE_NAME
                        Sample name
  --dump_pkl_dict       Dumps a pkl file of the splice event dict with all
                        quantifications
  --gzipped             Output files will be gzipped if set
  --event_ioe EVENT_IOE
                        Event ioe file - used to recover event-gene
                        association
  --calc_gene_frac      Requires IOE file to be passed to --event_ioe. If set,
                        a maximal event fraction of gene expression will be
                        estimated by max(min_excluded,
                        min_included)/max(gene_junctions)
  --suppress_output     Suppresses output files if set
  --enable_edge_use     Use individual exon edges that are unique to form in
                        quantification
  --turn_off_no_ends    Disable the exclusion of transcript-termini from
                        isoform-specific exon edge quantification (default is
                        to exclude ends)
  --suppress_eij_use    Don't use exon-overlapping exon-intron junctions for
                        quantification except for RI events.
  --disable_ri_extrapolation
                        Disable the use of RI/MR included forms to inform
                        quantification of other events that contain these
                        exons
  --n_bootstraps N_BOOTSTRAPS
                        Number of bootstraps. default = 0

DEXSeq_comparison.R

Description

DEXSeq_comparison.R statistically tests changes in PSI between two conditions with at least 2 replicates per condition.

Basic usage

Rscript /path/to/DEXSeq_comparison.R -e /path/to/infer_pairwise_events -p /path/to/jc_psi_files -s /path/to/sampleTable.csv 

Example sample_table.csv

Help statement

usage: DEXSeq_comparison.py [-h] [-e, --events INFER_PAIRWISE_EVENTS_DIR]
                                [-p, --psi_dir PSI_TSV_DIR]
                                [-s, --sample_table SAMPLE_TABLE]
                                [--min_jc MINIMUM_JC]
                                [--min_psi MINIMUM_PSI]

optional arguments:
  -h, --help            show this help message and exit
  --min_jc MINIMUM_JC   Minimum number of total mean junction read suport across
                        conditions to report an event. (Default: 15)
  --min_psi MINIMUM_PSI Minimum max(mean(PSI)) across conditions to report an
                        event. (Default: 0.1)
  --ri_span RI_SPAN_PSI Maximum max(mean(span_PSI)) across conditions to report
                        an RI/MR event. (Default: 0.03)

Output file descriptions

npc_esc_dpsi.csv

Field descriptions:

gene : gene symbol
event_id : unique name for event
event_type : type of event
chr : chromosome
start : start coordinate of event
end : end coordinate of event
strand : strand
conditionA_mean_ijc : (float) mean read count of all junctions in the included form across replicates in conditionA
conditionA_mean_ejc : (float) mean read count of all junctions in the excluded (skipped) form across replicates in conditionA
conditionA_mean_psi : (float) mean of PSI values computed for all pairs of included, excluded junction counts across replicates in condtionA
conditionB_mean_ijc : (float) mean read count of all junctions in the included form across replicates in conditionB
conditionB_mean_ejc : (float) mean read count of all junctions in the excluded (skipped) form across replicates in conditionB
conditionB_mean_psi : (float) mean of PSI values computed for all pairs of included, excluded junction counts across replicates in condtionB
dpsi : (float) difference of mean(PSI) between conditions
event_qval : (float) Q-value describing the positive false discovery rate of the difference of mean(PSI) between conditions.
sig : 1 = event is significantly more included (dpsi ≥ 0.1, event_qval ≤ 0.05) in condition B relative to condition A.
-1 = event is significantly less included (dpsi ≤ -0.1, event_qval ≤ 0.05) in condition B relative to condition A.
0 = event is not significantly different between conditions (|dpsi| < 0.1, event_qval > 0.05).