. Author manuscript; available in PMC: 2022 May 10.

Published in final edited form as: Nat Methods. 2022 Mar 3;19(3):307–310. doi: 10.1038/s41592-022-01400-x

The SpliZ generalizes ‘percent spliced in’ to reveal regulated splicing at single-cell resolution

Julia Eve Olivieri ^1,^2,⁴, Roozbeh Dehghannasiri ^1,^3,⁴, Julia Salzman ^1,³

PMCID: PMC9089759 NIHMSID: NIHMS1792899 PMID: 35241832

Abstract

Detecting single-cell-regulated splicing from droplet-based technologies is challenging. Here, we introduce the splicing Z score (SpliZ), an annotation-free statistical method to detect regulated splicing in single-cell RNA sequencing. We applied the SpliZ to human lung cells, discovering hundreds of genes with cell-type-specific splicing patterns including ones with potential implications for basic and translational biology.

Splicing is a core function of eukaryotic genomes that generates proteins with diverse and even opposite functions from a single gene¹, changes translation efficiency², controls localization³, and generates noncoding RNAs⁴. So far, very few genes’ splicing programs have been characterized at single-cell resolution, and the function of splicing remains a critical open problem in biology⁵. Constant advances in single-cell RNA-sequencing (scRNA-seq) technology now provide an unprecedented opportunity to understand how splicing is regulated at the single-cell level. However, the enormous complexity of splicing in eukaryotic genomes and the low sequencing depth per cell in scRNA-seq experiments, especially in droplet-based data, makes it extremely challenging to precisely quantify RNA isoform expression and its differential regulation in single cells.

With few exceptions⁶, the field has typically attempted to either estimate isoform expression using model-based approaches^7,8 and then perform differential splicing analysis, or directly quantify exon inclusion^9,10 using percent spliced in (PSI). Annotation-based isoform quantification methods give notoriously unstable point estimates in the presence of nonuniform read sampling⁷, incomplete annotations or low-depth and 3′-biased sequencing. The second approach, which is to use PSI, quantifies the fraction of transcripts skipping the exon¹¹. However, tests based on PSI must proceed exon by exon, requiring thousands of tests, and cannot detect splicing events beyond simple exon skipping. These issues have led to the view that robust differential splicing analysis in droplet-based scRNA-seq is out of reach^5,11,12.

Here, we introduce the splicing Z score (SpliZ), a scalar value assigned to each gene–cell pair that quantifies how deviant a cell’s splicing is compared to a population average. The SpliZ can be applied to compiled junctional read counts from any plate- or droplet-based single-cell data, aligned by STAR¹³ or SICILIAN¹⁴. The SpliZ integrates all nonconstitutive spliced reads on a per-gene basis to detect deviant splicing patterns in single cells under low and biased sampling, while remaining uncorrelated with gene expression (Fig. 1a and Extended Data Fig. 1). It also requires one test per gene, greatly reducing the number of tests required by PSI and has the power to detect isoform expression changes beyond exon skipping.

Fig. 1 ∣ — a, The SpliZ takes spliced reads and metadata and returns genes described as differentially spliced by cell type. b, Illustration of the SpliZ calculation as described in the Methods. c, SpliZ scores can be aggregated for each cell type, identifying cell types with statistically deviant splicing. The cell types with significant enrichment for short and long introns have large negative and positive median SpliZ values, respectively. d, Simulations of the SpliZ and PSI: reads of two cell types (20 cells each), each cell type with a unique fractional isoform abundance shown (1,000 trials). At each read depth n, Poisson(n) reads are sampled for each cell in proportion to isoform fractions. Cell populations with identical isoform expressions are used to estimate calls under the null hypothesis ( Methods). SpliZsites identified with the SVD are marked with asterisks and coincide with simulated alternative splice sites. Error bars are 95% binomial confidence intervals. The SpliZ increases power over PSI to detect cell types expressing different proportions of *CD47* isoforms. e, A representation of a splicing matrix, along with the first eigenvector from its SVD decomposition, with the SpliZsites marked with an asterisk. f, Simulation regime identical to d except for a different ground truth *TPM2*-like exon structure and isoform proportion. Error bars are 95% binomial confidence intervals. The SpliZVD calls differential alternative splicing between the two cell types with different proportions of *TPM2*-like isoforms with power greater than or equal to PSI, while the SpliZ does not.

The SpliZ quantifies splicing for each cell–gene pair by (1) assigning a rank to each read aligning to a splice junction based on the relative size of the intron compared to the set of observed introns for that splice site; (2) converting the rank to a residual measuring its statistical deviation compared to ‘the typical’ intron length rank and (3) statistically grounded scaling and summing of these values to a single SpliZ value per gene and cell (Fig. 1b). A cell–gene pair has a large negative (respectively positive) SpliZ, if, on average, introns have significantly smaller (respectively longer) intron length compared to the population of cells profiled in the experiment (Extended Data Fig. 2). Note that this definition of the SpliZ allows its values to change not only based on exon skipping, but also other events such as intronic alternative polyadenylation that could shift read coverage within the gene. The SpliZ is a mathematical generalization of PSI while increasing power in cases of more complicated exon skipping (Methods). Also, the theoretical properties of the SpliZ allow it to be efficiently integrated into significance testing: under the null hypothesis, the SpliZ has median 0 for every cell type, while under the alternative hypothesis the median SpliZ value per cell type is not necessarily zero (Fig. 1c and Supplementary Information).

As predicted by theory, the SpliZ has higher power than PSI in simulation when there are multiple isoforms for a gene, each of which includes different exons (Extended Data Fig. 3 and Fig. 1d). The SpliZ builds strength across these isoforms to identify real differences between cell types at lower read depths than PSI, while maintaining the same type II error rate (Fig. 1d). The Singular Value Decomposition (SVD) of the residual matrix is also used for biological interpretability. Splice sites corresponding to up to the three largest magnitude components of the first eigenvector are nominated as the statistically most variable splice sites or SpliZsites (Fig. 1e and Methods).

To increase the power of the SpliZ to detect alternative isoform expression in cases where the SpliZ has low power, we also developed the SpliZVD. The SpliZVD modifies the SpliZ by computing by projecting splicing residuals onto eigenvectors of its SVD (we consider only the first projection in this work) (Fig. 1f and Extended Data Fig. 4 and Methods). The SpliZVD maintains the power of PSI in simulation with two cassette exons (Fig. 1f), a situation where the SpliZ has no power to detect differences because both isoforms contain one ‘long’ and one ‘short’ intron. In these simulations, the splice sites involved in differential alternative splicing were correctly identified as SpliZsites (Extended Data Fig. 5).

We applied the SpliZ(VD) to 60,550 carefully annotated human lung cells across two individuals from the Human Lung Cell Atlas (HLCA)¹⁵. The HLCA, sequenced by the 10X (53,469 cells) and Smart-seq2 (7,081 cells) platforms, includes 57 annotated cell types with highly variable depth of sampling per cell type (Extended Data Fig. 6). Differential splicing analysis was performed on the 1,754 genes with a computable SpliZ in at least ten cells in one of the individuals (11,640 for Smart-seq2 data) (Methods and Supplementary Table 1).

The SpliZ and SpliZVD identified hundreds of differential alternative splicing events between 10X cell types in each individual. A total of 207 in individual 1 (respectively 219 in individual 2) genes were described as having significant cell-type-specific splicing by the SpliZ, and 135 (142) genes were called by the SpliZVD, with 88 (91) genes called by both (P value < 0.05, Supplementary Table 2). Also, 175 genes were called by either the SpliZ or SpliZVD in both individuals’ 10X data (P = 1.11 × 10⁻¹⁶, Extended Data Fig. 7). Restricting analysis to the genes that are significant in both 10X and Smart-seq2 data in the same individual, there is a positive correlation between median SpliZ scores for the same gene and cell type between 10X and Smart-seq2 data (Pearson correlation of 0.315 and 0.650, Extended Data Fig. 8), suggesting that SpliZ results are reproducible between technologies.

Genes with the most deviant SpliZ values include ATP5F1C, a gene encoding a core component of the mitochondrial ATP-synthase machinery, and MYL6, a gene encoding an essential component of the actin cytoskeleton with partially characterized splicing (Extended Data Fig. 9a,b and Supplementary Table 3). Other examples of highly cell-type-specifically spliced genes include LMO7, a gene encoding an emerin-binding protein with alternative splicing affecting a protein domain of unknown function (Fig. 2a), and CD47 in which we identified differential alternative splicing involving three isoforms affecting the cytoplasmic tail¹⁶ (Fig. 2b and Supplementary Table 3). CD47 encodes an immune-regulatory membrane protein that has recently been identified as a therapeutic target in a set of myeloid malignancies¹⁷ but has undescribed cell-type-specific splicing programs. One of the genes with the highest magnitude median SpliZ is PPP1R12A, a gene that encodes a protein phosphatase regulatory subunit, where the SpliZsite is not annotated as being alternatively spliced (Extended Data Fig. 9c). PPP1R12A is also called by the SpliZVD in Smart-seq2. This supports the use of the SpliZ to not only re-identify known alternative splicing events but also discover cell-type-specific splicing events through a purely statistical and annotation-free approach.

Fig. 2 ∣ — a, *LMO7* has significantly differential exon skipping affecting a domain of unknown function within epithelial cell types from splice site 75817161. Read fractions from two 10X and Smart-seq2 datasets are shown by dots and squares, respectively. Dot size represents the fraction of reads from each 5′ splice site to the 3′ site at 75817161. For each dot, the outer ring (in white) shows the upper confidence interval and the inner ring (color-coded) shows the lower confidence interval. Box plots show the distribution of the average alternative splice site per cell across all technologies and individuals. Each box shows 25–75% quantiles of average splice site per cell. The mid line of each box is the median and whiskers extend to 1.5 times the interquartile range. All points outside of this range are plotted individually. The cell types are grouped in two sets according to the sign of the median SpliZ for *LMO7*. The sashimi arcs for each cell-type group show the average splice site use across all cell types within that group and also across technologies and individuals. Thicker arcs correspond to higher fractional use measured in both 10X and Smart-seq2. b, Differential splicing of *CD47* is tissue-compartment-specific and affects the cytoplasmic tail protein domain. The alternative splicing involves one 3′ splice site (at 108047292) and four 5′ splice sites, including one unannotated splice junction (specified by vertical dashed lines). Cell types have highly tissue-compartment-specific inclusion of 5′ splice sites.

We compared the findings from the SpliZ analysis to Leafcutter¹⁸, a method designed to detect differential splicing in bulk RNA-seq. The intersection between genes called differentially spliced by the SpliZ and Leafcutter was greater than would be expected by chance (P < 10 × 10⁻⁶, hypergeometric test), providing orthogonal support for the validity of SpliZ calls (Methods and Supplementary Table 5). However, in both HLCA individuals, Leafcutter called fewer genes (92 compared to 207 for individual 1, 58 compared to 219 for individual 2, Methods). Drawbacks of using Leafcutter for scRNA-seq analysis are discussed in the Supplementary Information. This supports the idea that Leafcutter and similar bulk methods are generally not as appropriate for scRNA-seq analysis: in addition to their requirement of laborious preprocessing steps for pseudobulking, they require multiple pseudobulked samples from each cell type for reliable inference, which makes them statistically unfit for a wide range of scRNA-seq studies with only one or two 10X samples (even when each containing tens of thousands of cells).

Because the SpliZ has a tractable statistical distribution and is single-cell-resolved, it enables integrative biological analysis, including clustering approaches based on the SpliZ alone or in combination with gene expression. In addition, the SpliZ can be correlated with any other numerical phenotype, such as pseudotime, cell cycle trajectory, gene expression, and spatial transcriptomics¹⁹. Limitations of using the SpliZ on droplet-based data include restricted discovery due to 3′ bias (Supplementary Information). Creating scRNA-seq splicing methods that can discover other splicing changes, such as intron retention, is the subject of future work. In summary, by deconvolving technical and biological noise in splicing, the SpliZ provides a framework to identify and prioritize splicing events that are regulated and functional.

Methods

Preparing splice junction input files.

The scRNA-seq datasets were mapped to the reference human genome (hg38) using STAR v.2.7.5a with default parameters¹³.

We used SICILIAN¹⁴ with default settings for calling splice junctions from the STAR BAM files. SICILIAN is a statistical wrapper that can be applied to the alignment output file from a spliced aligner and can distinguish false positive junction calls from true positives via assigning a statistical score to each splice junction reported by the aligner.

SpliZ calculation.

The SpliZ score computation for a gene consists of two parts: one relative to the 3′ splice sites in the gene (the 3′ splice site SpliZ) and one relative to the 5′ splice sites (the 5′ splice site SpliZ). We first explain how to calculate the 3′ splice site SpliZ for one gene (suppressing the notation specifying the gene for simplicity) and the 5′ splice site SpliZ can be computed similarly. Let i specify the 3′ splice site, j specify the 5′ splice site, k specify the cell and ℓ specify the read. Therefore, each junctional read for the gene in the dataset is specified by a unique combination of ijkℓ. Note that we only consider junctions for which the 3′ splice site has multiple 5′ splice sites in the dataset.

The 3′ splice site SpliZ score calculation proceeds by treating each 3′ splice site separately. We will consider a plus strand gene and assume that 3′ splice site i has multiple 5′ splice sites across the whole dataset (otherwise we filter it out). We rank these 5′ splice sites in order from closest to the farthest from the 3′ splice site in question i (that is, from lowest to the highest genomic coordinate). For example, if there were four 5′ splice sites partnered with i across the whole dataset, we would rank them 1, 2, 3 and 4 in the order of their genomic coordinates. For genes that are on the minus strand, the 5′ splice sites are ranked differently as described in the Supplementary Information. Let r_ijkℓ denote the 5′ splice site rank for the read specified by ijkℓ. If r_ijkℓ = 1, it means the 5′ splice site has the smallest genomic coordinate among all 5′ splice sites in the dataset.

Now, let N_i be the number of junctional reads observed across all cells for 3′ splice site i. We can compute $\bar{r_{i}}$ , the average rank of the 5′ splice sites for 3′ splice site i, as

\bar{r_{i}} = \frac{\sum_{j, k, ℓ} r_{i j k ℓ}}{N_{i}} .

We can also find the variance of the ranks as:

{\hat{σ}}_{{\bar{r}}_{i}}^{2} = \frac{\sum_{j, k, ℓ} (r_{i j k ℓ} - \bar{r_{i}})^{2}}{N_{i}} .

Now, we renormalize the rank r_ijkℓ using the sample variance and mean as:

S_{i j k ℓ} = \frac{r_{i j k ℓ} - \bar{r_{i}}}{{\hat{σ}}_{\bar{r_{i}}}} .

We can see that E(S_ijkℓ) = 0 and var(S_ijkℓ) = 1 as we subtract the mean and divide by the standard deviation. The closer the 5′ splice site of a read is to the 3′ splice site, the smaller its corresponding S_ijkℓ value would be. In practice, we truncate the S_ijkℓ values at the tenth and 90th quantiles across all genes to avoid effects from extreme outliers.

We now aggregate these zero-mean and unit-variance variables for all 3′ splice sites in the gene to compute the 3′ splice site $SpliZ (z_{k}^{d})$ in the kth cell as:

z_{k}^{d} = \frac{\sum_{i, j, ℓ} S_{i j k ℓ}}{M_{k}},

where M_k is the number of junctional reads mapping to the gene in the kth cell. It is straightforward to see that E(z_k) = 0. Note that under the alternative hypothesis, for a given cell type E(S_ijkℓ) = μ ≠ 0. For a cell of this cell type,

E [z_{k}^{d}] = E [\frac{\sum_{i, j, ℓ} S_{i j k ℓ}}{M_{k}}] = \frac{\sum_{i, j, ℓ} E [S_{i j k ℓ}]}{M_{k}} = \frac{M_{k} μ}{M_{k}} = μ,

meaning that the SpliZ is not correlated with read depth M_k.

Knowing that the variance of the sum of independent random variables is the sum of their variances:

var (z_{k}^{d}) \approx \frac{\sum_{i j ℓ} var (S_{i j k ℓ})}{M_{k}^{2}} = \frac{1}{M_{k}} .

The approximation is due to the fact that the S_ijkℓs are not necessarily independent but we expect them to be close enough to independent.

Similarly, we can compute the 5′ splice site SpliZ $z_{k}^{a}$ . We average the two scores $z_{k}^{a}$ and $z_{k}^{d}$ to compute the SpliZ z_k for the gene in cell k:

z_{k} = (z_{k}^{d} - z_{k}^{a}) ∕ \sqrt{2} .

These values are subtracted to correct for signs, such that short introns correspond to small values and long introns correspond to high values for both. Division by $\sqrt{2}$ ensures that the variance will be comparable between averaged SpliZ values and SpliZ values for which only one of $z_{k}^{d}$ and $z_{k}^{a}$ is calculable (in which case that is the SpliZ value for the cell and there is no averaging).

For computing SpliZ, we only consider cells with at least six junctional reads mapping to the gene, and the junctional reads for which the 3′ splice site has more than one 5′ splice site observed in the dataset and vice versa.

Equivalence of the SpliZ and PSI.

We will prove the equivalence of the SpliZ and PSI under a specific scenario: when only one of the junctional reads of the exon inclusion event is measured; this could correspond to different ending exons of the transcript.

Let n₁ be the number of reads for a given cell mapping to the junction between the 3′ splice site and the first exon, and let n₂ be the number of reads mapping between the 3′ splice site and the second exon for that cell. Then the value of PSI for the cell is $ψ = \frac{n_{1}}{n_{1} + n_{2}}$ .

To calculate the SpliZ value z, we need to calculate the rank mean and standard deviation across the population of cells. Assume in the overall population the fraction of junctions including the exon is f. Then the mean rank $\bar{r} = \frac{1 \times (f) N + 2 \times (1 - f) N}{N} = 2 - f$ , where N is the total number of reads mapping to these two junctions across all cells. The variance σ² is given by

σ^{2} = \frac{(f) N (1 - (2 - f))^{2} + (1 - f) N (2 - (2 - f))^{2}}{N} = f (1 - f) .

Then the SpliZ value is given by

z = \frac{n_{1} (\frac{1 - (2 - f)}{\sqrt{f (1 - f)}}) + n_{2} (\frac{2 - (2 - f)}{\sqrt{f (1 - f)}})}{n_{1} + n_{2}} = \frac{1}{\sqrt{f (1 - f)}} (f - \frac{n_{1}}{n_{1} + n_{2}}) .

Therefore, z can be written in terms of ψ and f:

z = \frac{1}{\sqrt{f (1 - f)}} (f - ψ),

implying that z and ψ are equivalent in this case. For example, if f = 0.5 then z = 1 − 2ψ.

SpliZVD calculation.

Let M be an n × p matrix, where n is the number of cells and p is the number of splice sites for the given gene. Matrix entries are defined by:

M_{k i} = \frac{\sum_{ℓ \in L_{i k}} \tilde{S} (ℓ)}{∣ L_{i k} ∣} .

Here, L_ik is the set of reads using splice site i in cell k. $\tilde{S} (ℓ)$ is the normalized residual of the rank of read ℓ, defined as follows. For a given gene, let G be the set of spliced reads that map to the gene across the dataset. Let S(ℓ) be the residual of read l. Then let

μ = \frac{\sum_{ℓ \in G} S (ℓ)}{∣ G ∣},

and

σ^{2} = \frac{\sum_{ℓ \in G} (S (ℓ) - μ)^{2}}{∣ G ∣} .

then $\tilde{S} (ℓ) = \frac{S (ℓ) - μ}{σ}$ , so $E [\tilde{S} (ℓ)] = 0$ and $var \tilde{S} (ℓ) = 1$ . Let α^(j) be the jth eigenvector of M. Then the jth SpliZVD score for cell k is given by <m_k,α^(j)> (note that for this paper we only consider the SpliZVD score based on the first eigenvector, <m_k,^α⁽¹⁾>, and refer to this single score as the SpliZVD though more components can be considered). This score has mean 0 and variance $\frac{1}{∣ L_{i k} ∣}$ (j suppressed for simplicity in these calculations):

E [< m_{k}, α >] = E [\sum_{i = 1}^{p} α_{i} (\frac{\sum_{ℓ \in L_{i k}} \tilde{S} (ℓ)}{∣ L_{i k} ∣})] = \sum_{i = 1}^{p} α_{i} \frac{\sum_{ℓ \in L_{i k}} E [\tilde{S} (ℓ)]}{∣ L_{i k} ∣} = 0,

and

var (< m_{k}, α >) = var (\sum_{i = 1}^{p} (α_{i} \frac{\sum_{ℓ \in L_{i k}} \tilde{S} (ℓ)}{∣ L_{i k} ∣})) = \sum_{i = 1}^{p} α_{i}^{2} \frac{\sum_{ℓ \in L_{i k}} var (\tilde{S} (ℓ))}{∣ L_{i k} ∣^{2}} = \frac{1}{∣ L_{i k} ∣} \sum_{i = 1}^{p} α_{i}^{2} = \frac{1}{∣ L_{i k} ∣} .

Note, this is assuming the residuals are uncorrelated. Under the alternative hypothesis for a cell type, as with the SpliZ, for an individual cell type if $E [\tilde{S} (ℓ)] = μ \neq 0$ :

E [< m_{k}, α >] = \sum_{i = 1}^{p} α_{i} \frac{\sum_{ℓ \in L_{i k}} E \tilde{S} (ℓ)}{∣ L_{i k} ∣} = \sum_{i = 1}^{p} α_{i} \frac{μ ∣ L_{i k} ∣}{∣ L_{i k} ∣} = \sum_{i = 1}^{p} α_{i} μ .

Therefore, the SpliZVD value is not dependent on read depth.

Calling differentially alternatively spliced genes by cell type.

We perform the following test independently on each gene for each individual, technology and score (SpliZ and SpliZVD). The procedure is the same for both SpliZ and SpliZVD, but we will only discuss the SpliZ here for simplicity. For a given gene, we only consider cell types with more than ten cells with SpliZ values for that gene.

Significance of alternative splicing between cell types in a gene is determined by calculating P values using a two-step procedure based on the work in ref. ²⁰. Consider the following equation from ref. ²⁰:

T_{n, 1} = \sum_{i = 1}^{k} \frac{n_{i}}{{\hat{σ}}_{n, i}^{2}} {[{\hat{θ}}_{n, i} - \frac{\sum_{i = 1}^{k} n_{i} {\hat{θ}}_{n, i} ∕ {\hat{σ}}_{n, i}^{2}}{\sum_{i = 1}^{k} n_{i} ∕ {\hat{σ}}_{n, i}^{2}}]}^{2} .

We can compute T_n,1 based on the SpliZ values for all cells with splicing expression for the given gene as follows: the k samples represent the k cell types with splicing expression for the given gene. Then ${\hat{θ}}_{n, i} = median (X_{i, 1}, \dots, (X_{i, n_{i}})$ is our test statistic, where X_i,1, …, X_{i,n_i} are the SpliZ values for the n_i cells from cell type i with splicing expression for the gene. $σ_{n, i}^{2}$ is the sample standard deviation of the SpliZ values for cell type i; in practice, we inflate the variance to $σ_{n, i}^{2} + 0.1$ for robustness. Our null hypothesis is that all cell types have the same median SpliZ for this gene.

Performing all permutations of cell types to cells and recalculating T_n,1 yields the permutation distribution (a subset of the permutations yields an approximation of the permutation distribution). Theorem 3.1 in ref. ²⁰ states that under some assumptions, this permutation distribution converges to the χ² distribution with k − 1 degrees of freedom. Also, if the sample distributions do not have different medians, the probability that the permutation test rejects the null hypothesis tends to nominal level α.

It is a lot quicker in practice compare to the χ² distribution than compute permutations; therefore, we first calculate the P value based on comparing to the χ² distribution (P_χ²). Then if Pχ² < 0.05, we compute the permutation P value P_perm by permuting the assignments of cell types to cells for the gene and calculating $T_{n, 1}^{(j)}$ based on this permuted data where j is the current permutation. This results in a permutation null distribution $T_{n, 1}^{(1)}, \dots, T_{n, 1}^{(J)}$ where J is the number of permutations performed. Then

cdfperm = \frac{\sum_{j = 1}^{J} I {T_{n, 1}^{(j)} < T_{n, 1}}}{J},

and P_perm = 2 min(cdf_perm, 1 − cdf_perm) (quantifies whether the real value is extreme in either direction). We then adjust the P values for multiple hypothesis testing using the Benjamini–Hochberg procedure. Genes are called significant if their adjusted P values are less than 0.05.

Genes with the ‘most deviant SpliZ values’ are defined as follows: find the largest magnitude median SpliZ value Z_i across all cell types for each gene i. Then genes with higher Z_i values have more ‘deviant’ SpliZ values.

Finding the SpliZsites driving the alternative splicing of a gene.

To provide further biological interpretation for the SpliZ and SpliZVD scores and automatically identify the SpliZsites driving the alternative splicing in a significantly regulated gene, we take the SVD of the residual matrix and use its eigenvalues and eigenvectors to select the SpliZsites. To do so, we select up to three eigenvectors depending on the values of their corresponding eigenvalues and then in each chosen eigenvector we define up to three splice sites that correspond to up to three elements of the eigenvector with the highest absolute value as the SpliZsite. If the first eigenvalue is greater than 0.7, we select only the first eigenvector and otherwise we select the second and third eigenvectors if their corresponding eigenvalues are at least 0.2. When an eigenvector is selected, we choose the three splice sites corresponding to the top entries that have at least 10% of the eigenvector loadings based on the L² norm of the eigenvector.

For each driving splice site, we additionally report its annotation status, including whether it is annotated as an alternatively spliced exon. To determine whether an exon is known to be involved in alternative splicing, we extracted the splice sites from the hg38 annotation GTF file and then considered those 3′ splice sites (respectively 5′ splice sites) that are observed to be spliced to more than one distinct 5′ splice site (respectively 3′ splice site) across the extracted junctions as known alternatively spliced sites.

Simulation methods.

For each simulation, two cell types are simulated with 20 cells in each. They have predefined proportions of each isoform as described in Fig. 1d,f. At each mean read depth n, 100 trials are performed, each of which proceeds as follows. First, a read depth for each cell k is sampled independently from Poisson(n). Then for each cell of cell type c, the distribution of the k reads among splice junctions is drawn from a multinomial distribution, for which each junction from isoform i has the probability $f_{i}^{(c)} ∕ m_{i}$ , where $f_{i}^{(c)}$ is the underlying fraction of isoform i in cell type c and m_i is the number of exons in isoform i.

Next, the SpliZ and SpliZVD are calculated as described above for this gene and this population of 40 cells. PSI is calculated on a cell-by-cell basis for each cassette exon by dividing the number of reads including that exon over all of the reads that either skip or include the exon in the given cell (note that some reads may not either skip or include the exon; these are disregarded in the PSI calculation).

Then P values are calculated separately for the SpliZ, SpliZVD and each PSI value as described above, except χ² filtering is not used and variances are not inflated. Benjamini–Hochberg multiple hypothesis testing correction is not used (there is only one score calculated for each of the SpliZ and SpliZVD; if the multiple PSI values were corrected, it would only cause them to be less significant). The gene is described as having differential alternative splicing between cell types if the calculated P <0.05. SpliZsites are calculated for each by simulating the two cell types at a read depth of 20 and performing the procedure described in the methods.

For each simulation, a ‘null’ simulation is also included, which follows the same setup except both cell types have the same distribution among the isoforms as the first cell type in the original simulation. This allows us to estimate the type I error rate.

Smart-seq2 data processing.

To make the Smart-seq2 data most suitable for validation, we subsetted the dataset to only include junctions that were ever found in 10X data and then ran it through the pipeline (Extended Data Fig. 10). To calculate the correlation between the median SpliZ scores we subsetted each 10X and Smart-seq2 dataset separately for each individual, only including junctions shared by both datasets for that individual, and only cell types shared by both technologies for that individual. Note that the proportions of cell types were not matched.

Concordance between datasets.

We test whether the number of genes significant in both 10X datasets is more than expected under the null hypothesis, which is that genes found significant in both datasets are due to chance. Let s be the number genes that are called significant in both individual 1 and 2’s 10X data. The probability that a given gene is significant in both datasets is:

P (gene significant in both) = P (gene significant in individual 1) \times P (gene significant in individual 2 ∣ gene significant in individual 1)

Under the null hypothesis, assume that a gene being significant in one individual is independent from it being significant in the other individual. Therefore, under the null hypothesis

P (gene significant in both) = P (gene significant in individual 1) \times P (gene significant in individual 2)

We can estimate the quantities on the right-hand side of the question for each individual:

P_{i} = (no genes called in individual i) ∕ (total num genes with computable SpliZ in individual i)

Therefore, under the null hypothesis, the probability that ≥x genes are significant in both individuals out of n genes is given by

1 - binom_cdf (x, n, P_{1} P_{2}) .

Correlation between datasets is calculated by considering only gene/cell-type pairs with at least ten cells in both datasets, restricting to genes called significant in both datasets and calculating the Pearson correlation.

Protein domain annotation.

Protein domains for analyzed genes were determined by finding whichever protein domains overlapped with the chromosome and coordinates of the splicing in the Pfam file (location downloaded described in the Data Availability section). The CD47 cytoplasmic domain was added based on^16,21.

Extended Data

Extended Data Fig. 3 ∣ — In both simulations, two cell types with 20 cells respectively are simulated, each type having a different proportion of isoforms (1000 trials for a, 100 trials for b). At each read depth, Poisson(n) reads are sampled in proportion to isoform abundances. Null values are calculated from cell populations with identical isoform expression. SpliZsites based on the SVD described in the methods are starred by asterisks, and coincide with all splice sites with differential partner exon usage between the two cell types. a. The SpliZ, SpliZVD, and PSI all have the same power in the case of exon skipping. b. The simulation from Fig. 1d was modified, changing the fractions of isoform abundance as described in the figure such that when exon 2 is present, exon 3 is included 99% of the time. While this splice site was identified as a SpliZsite in Fig. 1d, here it is not identified (red X), as would be expected for a splice junction with 99% constitutive splicing. This simulation setting shows that the SpliZ can correctly identify only alternative exons as SpliZsites.

Extended Data Fig. 4 ∣ — The SpliZVD is the projection of the matrix of splicing residuals onto its first eigenvector. The residual matrix’s SVD is used to identify the most variably alternatively spliced sites. The top left shows a gene structure of *TPM2*, with reads aligning to different junctions above. Each read is assigned a SpliZ residual and SpliZVD residual, the latter of which is based on the first eigenvector (shown below the gene annotations). Ovals representing different cells are colored based on the sign of their SpliZ or SpliZVD values, showing that in this case the SpliZVD is able to distinguish differences in splicing where the SpliZ cannot.

Extended Data Fig. 5 ∣ — SVD decompositions of the splicing residual matrix based on the simulations in Fig. 1d,f, and Extended Data Fig. 3 at average read depth 20. Rows correspond to splice sites and columns correspond to eigenvectors. For each splice site, the value from the SVD that causes the splice site to be picked as SpliZsite is boxed in green. a,b,c. As expected, all non-constitutive splice sites are selected as SpliZsites. d. As expected, the three most variable splice sites are chosen as SpliZsites, while the 99% constitutive splice site is not chosen as SpliZsite.

Extended Data Fig. 6 ∣ — a. Box plot of the number of reads per cell for each cell type in both individuals’ 10x data. Numbers to the right of the plot indicate the number of cells plotted for each cell type over the given experiment. The middle line of each box is the median, and each box extends from the first to third quartile. Whiskers extend to 1.5 times the interquartile range. All points outside of this range are plotted individually. b. Bar plots of the number of cells per cell type for each individuals’ 10x data.

Extended Data Fig. 7 ∣ — a. Intersection of genes called by the SpliZ between 10X and subsetted Smart-seq2 data. b. Intersection of genes called by the SpliZVD between 10X and subsetted Smart-seq2 data. c. Intersection between genes called by the SpliZ and SpliZVD in 10x data.

Extended Data Fig. 8 ∣ — Correlation between median SpliZ values per matched gene and cell type a. Between individuals for 10X data. b. Between technologies for each HLCA individual (datasets subsetted to shared junctions and cell types per individual before running).

Extended Data Fig. 9 ∣ — The compartment-specific regulated alternative splicing of (a) *ATP5FC1*, (b) *MYL6*, and, (c) *PPP1R12A*. For *ATP5F1C*, endothelial and stromal cells have a higher fraction of junctional reads for the exon exclusion event compared to other compartments. For *MYL6*, immune cells have lower fraction of junctional reads for the exon inclusion event compared to other compartments. For *PPP1R12A*, immune cells have a higher fraction of junctional reads for the exon inclusion event. The plots include the cell types in which the splice site has at least 20 junctional reads in at least 10 cells in at least one of the 4 datasets (two individuals and two technologies: 10X and SS2) were chosen. Dots represent the fraction of junctional reads for each alternative splice site in the celltype. For each dot, the outer ring (in white) shows the upper CI and the inner ring (color-coded) shows the lower CI. Box plots show the average alternative splice site for each cell across all technologies and individuals. Each box shows 25-75% quantiles of average splice site per cell. Arcs represent the average fractions at the compartment level. There is no dot when the alternative splice site has 0 junctional reads.

Extended Data Fig. 10 ∣ — a. Upset plot showing the number of cell types sequenced in each HLCA individual and technology. b. Bar plot showing the number of cells sequenced for each HLCA individual and technology. c. Upset plot showing the number of alternative junctions (junctions for which at least one splice site has at least one other partner in the dataset) for each HLCA individual and technology. d. The SpliZ is calculated independently for Smart-seq2 data restricted to junctions detected by 10X to measure technology-dependence of results.

Supplementary Material

Supplementary Information

NIHMS1792899-supplement-Supplementary_Information.pdf^{(545.5KB, pdf)}

Supplementary Tables 1--3

NIHMS1792899-supplement-Supplementary_Tables_1--3.xlsx^{(2.2MB, xlsx)}

Acknowledgements

We thank P. Wang and S. Quake for insightful and instrumental comments during the development of the method, E. Meyer and R. Bierman for comments on the manuscript, J. Klein for creating parts of Fig. 1, and K. Travaglini and M. Krasnow for providing advanced access to the HLCA data before its publication. J.O. is supported by the National Science Foundation Graduate Research Fellowship under grant no. DGE-1656518, a Stanford Graduate Fellowship and a Lieberman Fellowship. R.D. is supported by the Cancer Systems Biology Scholars Program grant no. R25 CA180993 and Clinical Data Science Fellowship grant no. T15 LM7033-36. J.S. is supported by the National Institute of General Medical Sciences grant nos. R01 GM116847 and R35 GM139517 and the National Science Federation Faculty Early Career Development Program Award no. MCB1552196.

Footnotes

Competing interests

The authors declare no competing interests.

Extended data are available for this paper at https://doi.org/10.1038/s41592-022-01400-x.

Supplementary information The online version contains supplementary material available at https://doi.org/10.1038/s41592-022-01400-x.

Reporting Summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Code availability

The SpliZ code along with the code used for data analysis and to create the figures are available through a GitHub repository https://github.com/juliaolivieri/SpliZ_pipeline/. This repository is archived with Zenodo under the following https://doi.org/10.5281/zenodo.5781783 (ref. ²³). The pipeline was written in Python (v.3.6.7), and installed package versions are the following (also available in an environment. yml file on github): matplotlib (v.2.2.3); numpy (v.1.18.4); pandas (v.1.0.4); pyarrow (v.0.15.1); scipy (v.1.4.1); snakemake-minimal (v.5.4.5); statsmodels (v.0.11.1) and tqdm (v.4.46.0). We used Leafcutter (https://github.com/davidaknowles/leafcutter) and regtools (https://github.com/griffithlab/regtools).

Online content

Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41592-022-01400-x.

Data availability

HLCA data was downloaded from the European Genome-Phenome Archive at accession number EGAS00001004344 (ref. ²²). We refer to patient 2 in HLCA as individual 1 and patient 3 as individual 2 in our analysis. Our cell-type definition is based on concatenating the ‘compartment’ and ‘free annotation’ columns from the HLCA metadata and only considering lung cells (not blood). SpliZ scores and Leafcutter results, as well as the original data needed to reproduce these results, are available at FigShare: https://doi.org/10.6084/m9.figshare.14378819.v1. Human RefSeq hg38 annotation file was downloaded from ftp://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/annotation/GRCh38_latest/refseq_identifiers/GRCh38_latest_genomic.gff.gz. The UCSC Pfam database for the hg38 genome assembly was downloaded from http://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/ucscGenePfam.txt.gz.

References

1.Shao Y et al. Alternative splicing-derived intersectin1-l and intersectin1-s exert opposite function in glioma progression. Cell Death Dis. 10, 431 (2019). [DOI] [PMC free article] [PubMed] [Google Scholar]
2.Nakka K, Kovac R, Wong MM-K & Dilworth FJ Intron retained, transcript detained: intron retention as a hallmark of the quiescent satellite cell state. Dev. Cell 53, 623–625 (2020). [DOI] [PubMed] [Google Scholar]
3.Oleynikov Y & Singer RH RNA localization: different zipcodes, same postman? Trends Cell Biol. 8, 381–383 (1998). [DOI] [PubMed] [Google Scholar]
4.Yang Y & Carstens RP Alternative splicing regulates distinct subcellular localization of epithelial splicing regulatory protein 1 (esrp1) isoforms. Sci. Rep 7, 3848 (2017). [DOI] [PMC free article] [PubMed] [Google Scholar]
5.Arzalluz-Luque Á & Conesa A Single-cell RNAseq for the study of isoforms—how is that possible? Genome Biol. 19, 110 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
6.Vaquero-Garcia J et al. A new view of transcriptome complexity and regulation through the lens of local splicing variations. eLife 5, e11752 (2016). [DOI] [PMC free article] [PubMed] [Google Scholar]
7.Salzman J, Jiang H & Wong WH Statistical modeling of RNA-seq data. Statistical Sci. 10.1214/10-STS343 (2011). [DOI] [PMC free article] [PubMed] [Google Scholar]
8.Li JJ, Jiang C-R, Brown JB, Huang H & Bickel PJ Sparse linear modeling of next-generation mRNA sequencing (RNA-seq) data for isoform discovery and abundance estimation. Proc. Natl Acad. Sci. USA 108, 19867–19872 (2011). [DOI] [PMC free article] [PubMed] [Google Scholar]
9.Trincado JL et al. Suppa2: fast, accurate, and uncertainty-aware differential splicing analysis across multiple conditions. Genome Biol. 19, 40 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
10.Shen S et al. rmats: robust and flexible detection of differential alternative splicing from replicate RNA-seq data. Proc. Natl Acad. Sci. USA 111, E5593–E5601 (2014). [DOI] [PMC free article] [PubMed] [Google Scholar]
11.Najar CFBA, Yosef N & Lareau LF Coverage-dependent bias creates the appearance of binary splicing in single cells. eLife 9, e54603 (2020). [DOI] [PMC free article] [PubMed] [Google Scholar]
12.Westoby J, Artemov P, Hemberg M & Ferguson-Smith A Obstacles to detecting isoforms using full-length scRNA-seq data. Genome Biol. 21, 74 (2020). [DOI] [PMC free article] [PubMed] [Google Scholar]
13.Dobin A et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013). [DOI] [PMC free article] [PubMed] [Google Scholar]
14.Dehghannasiri R, Olivieri JE, Damljanovic A & Salzman J Specific splice junction detection in single cells with SICILIAN. Genome Biol. 22, 219 (2021). [DOI] [PMC free article] [PubMed] [Google Scholar]
15.Travaglini KJ et al. A molecular cell atlas of the human lung from single-cell RNA sequencing. Nature 587, 619–625 (2020). [DOI] [PMC free article] [PubMed] [Google Scholar]
16.Hayat SMG et al. Cd47: role in the immune system and application to cancer therapy. Cell. Oncol 43, 19–30 (2020). [DOI] [PMC free article] [PubMed] [Google Scholar]
17.Chao MP et al. Therapeutic targeting of the macrophage immune checkpoint CD47 in myeloid malignancies. Front. Oncol 9, 1380 (2020). [DOI] [PMC free article] [PubMed] [Google Scholar]
18.Li YI et al. Annotation-free quantification of RNA splicing using leafcutter. Nat. Genet 50, 151–158 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
19.Olivieri JE et al. RNA splicing programs define tissue compartments and cell types at single cell resolution. eLife 10, e70692 (2021). [DOI] [PMC free article] [PubMed] [Google Scholar]
20.Chung E & Romano JP Exact and asymptotically robust permutation tests. Ann. Stat 41, 484–507 (2013). [Google Scholar]
21.Li YI et al. Annotation-free quantification of RNA splicing using leafcutter. Nat. Genet 50, 151–158 (2018). [DOI] [PMC free article] [PubMed] [Google Scholar]
22.Travaglini KJ et al. A molecular cell atlas of the human lung from single-cell RNA sequencing. Nature 587, 619–625 (2020). [DOI] [PMC free article] [PubMed] [Google Scholar]
23.Olivieri JE juliaolivieri/SpliZ_pipeline: v1.0. Zenodo; 10.5281/zenodo.5781783 (2021). [DOI] [Google Scholar]

Associated Data

This section collects any data citations, data availability statements, or supplementary materials included in this article.

Supplementary Materials

Supplementary Information

NIHMS1792899-supplement-Supplementary_Information.pdf^{(545.5KB, pdf)}

Supplementary Tables 1--3

NIHMS1792899-supplement-Supplementary_Tables_1--3.xlsx^{(2.2MB, xlsx)}

PERMALINK

The SpliZ generalizes ‘percent spliced in’ to reveal regulated splicing at single-cell resolution

Julia Eve Olivieri

Roozbeh Dehghannasiri

Julia Salzman

Abstract

Fig. 1 ∣. The SpliZ outperforms PSI in simulation.

Fig. 2 ∣. The SpliZ detects cell-type-specific splicing in HLCA dataset.

Methods

Preparing splice junction input files.

SpliZ calculation.

Equivalence of the SpliZ and PSI.

SpliZVD calculation.

Calling differentially alternatively spliced genes by cell type.

Finding the SpliZsites driving the alternative splicing of a gene.

Simulation methods.

Smart-seq2 data processing.

Concordance between datasets.

Protein domain annotation.

Extended Data

Extended Data Fig. 1 ∣. The SpliZ is not correlated with gene expression.

Extended Data Fig. 2 ∣. SpliZ toy example.

Extended Data Fig. 3 ∣. Simulations of the SpliZ, SpliZVD, and PSI.

Extended Data Fig. 4 ∣. SpliZVD calculation example.

Extended Data Fig. 5 ∣. Splicing residual matrices.

Extended Data Fig. 6 ∣. Cell type summary in HLCA dataset.

Extended Data Fig. 7 ∣. Called gene intersections.

Extended Data Fig. 8 ∣. SpliZ correlations.

Extended Data Fig. 9 ∣. Compartment-specific splicing.

Extended Data Fig. 10 ∣. 10X vs Smart-seq2 comparison.

Supplementary Material

Acknowledgements

Footnotes

Data availability

References

Associated Data

Supplementary Materials

Data Availability Statement

ACTIONS

PERMALINK

RESOURCES

Similar articles

Cited by other articles

Links to NCBI Databases