Abstract
Introduction:
Alternative splicing (AS) of mRNAs is a highly conserved mechanism which can greatly expand the functional diversity of the transcriptome. Aberrant splicing underpins many diseases, and a better understanding of AS can provide insights regarding the molecular mechanisms involved. Importantly, AS can be affected by genetic variants and several studies have indicated large numbers of splicing quantitative trait loci (sQTL). With the advance of single-cell technology, expression QTL studies have been expanded to identify cell type level variants.
Methods:
We collected eight full-length scRNA-seq pancreatic islet datasets. Genotyping for each individual was done by the CTAT pipeline and Streka2. The isoform quantification was done by RSEM. Finally, sQTL was obtained by sQTLseeker2.
Results:
As a result, we identified 228 cell type level sQTLs for alpha and beta cells across 152 genes. In particular, our study highlights four variants affecting CDC42, a gene related to cell morphology, which have not been observed from bulk sQTL analysis.
Discussion:
Our results provide a proof of concept that it is possible to identify cell type level sQTLs, and we envision that better powered studies will allow us to further uncover the genetic regulation of splicing.
1 Introduction
Alternative splicing (AS) is pervasive in the human genome, and it allows for a remarkable expansion of protein diversity (; ). AS is a highly conserved mechanism, and aberrant splicing has been implicated in a wide range of diseases (Lopez-Bigas et al., 2005; ; Wang and Lee, 2018; ). A better understanding of AS will not only allow us to appreciate the role of transcriptome diversity, but it could also inform the development of novel therapies.
One successful strategy for identifying functional mechanisms and potential drug targets has been through naturally occurring genetic variants. Specific variants can be linked to a trait through genome-wide association studies (GWAS) (Manolio, 2010; Sollis et al., 2023). In addition, variants have also been linked to molecular phenotypes and expression quantitative trait loci (eQTL) analysis has been widely used to show a link between SNPs and gene expression (; ). Similarly, splicing QTL (sQTL) analysis has been used to link SNPs to the alteration of splicing events (Monlong et al., 2014; Garrido-Martin et al., 2021; Yamaguchi et al., 2022). Previously, more than 30 studies with 127 datasets were reported in the 2023 eQTL catalogue (Kerimov et al., 2023) which contains information from >7,500 donors and >1.7 million fine mapped associations. There are also several splicing studies, including GTEx (Zhang et al., 2020; Garrido-Martin et al., 2021; ; ; Qi et al., 2022; Yamaguchi et al., 2022) which contains >210,000 sQTLs from >600 donors.
With the advent of single-cell technologies, eQTL studies have been expanded to reveal cell type specific variants (van der Wijst et al., 2020; ; Neavin et al., 2021; Nathan et al., 2022; Yazar et al., 2022). Neavin et al. (2021) studied single-cell eQTLs in fibroblasts and fibroblast derived iPSCs from a total of 110 donors, and they found 71,926 reprogramming related eQTLs. They observed numerous cell type-specific eQTLs, while only 41.1% were observed in bulk-eQTL results from GTEx. Yazar et al. (2022) conducted single-cell eQTL analysis with PBMCs collected from 982 donors, and found loci related to the autoimmune disease at the cellular level. They also observed immune cell-specific eQTLs, of which only 40%–63% were observed in bulk PBMC data. Nathan et al. (2022) studied single-cell eQTLs in memory T cells with 259 individuals, and they also reported cell state specific eQTLs, of which 80%–90% were observed from bulk data.
However, cell-type level sQTLs have not yet received the same amount of attention. The main reason why this has not been attempted previously is because sQTL studies require both a sufficient number of samples as well as good coverage of splice junctions. A sufficient sample size is critical for those QTL analyses since previous analyses using bulk data in the GTEX database (Shan et al., 2019) have shown that the number of significant eQTLs is strongly correlated with the sample size (Supplementary Figure S1A), indicating that the number of donors is a key factor for sQTL analysis. Although there are studies profiling a large number of individuals, they employ single cell technologies that detect only the 3′- or 5′-ends of each transcript, and thus there is not enough information to quantify splice junction usage (; Westoby et al., 2018). Despite this limitation, a recent study reported splicing QTL at the single-cell level using 5′-end 10X scRNA-seq for PBMC data. The authors focused on ancestry-related sQTL,the reads that were biased to the 5′ region, which may lead to poor coverage for detecting other splicing events, something the authors mentioned as a limitation (Tian et al., 2024). Since the platforms that are capable of profiling full-length transcripts are more expensive, there are not many datasets suitable for sQTL analyses available. An additional obstacle for sQTL analyses comes from the fact that donor genotypes are required. Even though the gold standard is to use WGS or a chip to profile each donor, genotyping by using full-length scRNA-seq was recently reported, and benchmarks demonstrated an overall F1-score >0.95 (Liu et al., 2019). Moreover, it has been shown that cell type specific splicing events can be reliably identified via full-length scRNA-seq (Westoby et al., 2018). Combining these approaches, we collected eight pancreatic islet datasets to obtain a sample size sufficient for performing single-cell sQTL analysis. From a previous bulk RNAseq study of sQTLs in pancreatic islets which used 399 human donors, 4,858 sQTLs from 2,058 sGenes were reported (). Based on the association between those sQTLs and known T1D or T2D risk variants, the authors concluded that there is a functional association of T1D or T2D disease due to aberrant splicing events. The authors also observed an association between sGenes and fatty acid biosynthesis () or mTORC1 pathway (; ; ). These implications are important since these pathways had previously been associated with cell types in pancreatic islets. By using single-cell level sQTL analysis, we can further link those sQTLs in a cell type-specific manner.
2 Methods
2.1 Data acquisition
The dataset for islet atlas was obtained from Camunas-Soler et. al. (islet1) (), Segerstolpe et. al. (islet2) (Segerstolpe et al., 2016), Enge et. al. (islet3) (), Xin et. al. (islet4) (Xin et al., 2016), Dorajoo et. al. (islet5) (), Lawlor et. al. (islet6) (Lawlor et al., 2017), Wang et. al. (islet7) (Wang et al., 2016), and Marquez-Curtis et. al. (islet8) (Marquez-Curtis et al., 2022). For islet1, we excluded 3 T1D patients due to the small sample size. We excluded the R269 sample since the raw data was not available. For islet3, we only used adult samples. For islet5, we excluded “subject 3” sample since the raw data was not available. For islet7, we excluded HP15041 due to genotyping problems. We used cell-type information from the original papers. For islet6, we used marker genes threshold for cell typing based on Supplementary Table S3 in the original paper. We only considered major islet cell types: alpha, beta, gamma, delta, acinar, and ductal cells.
2.2 Genotyping
We used the Trinity Cancer Transcriptome Analysis Toolkit (CTAT) which is based on the GATK Best Practices pipeline (Liu et al., 2019; ) for the variant calling of single cells in scRNA-seq data (). We ran the pipeline with either “single-end” or “paired-end” for each corresponding dataset. We set the parameter “boosting_method”, which speeds up during the variant calling, to be “none”. We excluded SNVs from the known RNA-editing sites from RADAR (Ramaswami and Li, 2014) and RediProtal (Picardi et al., 2017) databases, which were provided by the CTAT pipeline. During the process, 13/1,464 islet1 cells, 24/2,038 islet2 cells, 13/1,578 islet3 cells, 13/1,507 islet4 cells, 222/333 islet5 cells, 9/593 islet6 cells, 59/268 islet7 cells, and 1/147 islet8 cells were excluded due to lack of confident SNVs or absence of raw data from a certain sample. After obtaining all the SNVs from single cells, we merged per donor using bcftools (1.13) (). Since somatic mutations could exist, we ran a germline mutation calling pipeline from Strelka2 (2.9.2) (Kim et al., 2018) with merged bam files (using Samtools; v1.11) () for each donor to identify germline mutations. We only keep the “passed” germline mutation from Strelka2 overlapping the merged mutation profiles at the donor level from the CTAT pipeline. We only keep the genotype information if there is a read at that region. We excluded X, Y, and M chromosomes. We excluded variants with minor allele frequency (MAF) < 0.1 (genomicMateSelectR v0.2.0) (Wolfe, 2022) from both healthy individuals and T2D patients. Locational information for SNP was obtained from CTAT result.
2.3 Transcript and gene quantification
We obtained the transcript and gene expression matrix by running RSEM (Li and Dewey, 2011; Westoby et al., 2018) (1.3.0) using rsem-calculate-expression--estimate-rspd--no-bam-output -p 16 with GRCh37 version 19 (Ensembl 74) (transcript file: *.isoforms.results, gene file: *.genes.results). As a result, we could align reads to 38,381 transcripts and 20,392 genes. The mean number of expressed genes per cell was 6,196 and the mean number of expressed transcripts per cell was 6,835.
2.4 Genomic location of transcript
We extracted genomic region of each gene by GenomicFeatures package (v 3.17) with [type = = “gene”, gene_status ! = “PUTATIVE”, no duplicated gene name] (Lawrence et al., 2013).
2.5 scRNA-seq processing
We used log2 (TPM + 1) matrix of transcript from RSEM for running Seurat (v 4.1.0) ()pipeline (FindVariableFeatures, ScaleData, RunPCA, RunUMAP [dims = 1:30]). We corrected batch effect of different dataset by running Harmony (v 0.1.0) (Korsunsky et al., 2019). We performed clustering by FindNeighbor (reduction = “harmony”, dims = 1:30) and FindClusters (resolution = 1.2). We excluded misaligned cell type from each cluster to remove batch outliers. We generated pseudobulk cell type-specific transcript matrix by averaging the transcript expression level per sample. After that, we used inverse log transformation (2^x −1) for sQTLseeker2 input since they suggested using a non-log transformed matrix to better calculate Hellinger distance for transcript expression.
2.6 sQTL analysis
We used the sQTLseekeR2 (v1.1.0) (Garrido-Martin et al., 2021) package for cell type level sQTL. We filtered out the lowly expressed transcript and gene by “min.transcript.exp = 0.1” and “min.gene.exp = 0.1” using prepare.trans.exp function. We used sqtl.seeker function to obtain nominal p-value with “svQTL = TRUE” (to remove a potential false positive sQTL by inhomogeneity in the variance from different genotype groups), “genic.window = 50,000 (bp)”, “min.nb.ind.geno = 5 (at least five samples for a certain genotype)”, “ld.filter = 0.8 (Linkage disequilibrium threshold)”. Age and gender were used for covariables to regress the transcript expression level. For multiple hypotheses correction, we used sqtls function with “FDR = 0.1” and “FDR.svQTL = 0.01” with default parameters (md.min: the minimum MD (Maximum Difference) in relative expression: 0.05).
2.7 Differential expressed gene analysis
For DEG analysis, we normalized the gene expression from the Seurat package (v 4.1.0) () and performed the FindMarkers function with logfc.threshold = 0.25 and min.pct = 0.1. We further used p.adjusted < 0.05.
2.8 eQTL analysis
We used FastQTL (v2.184) (Ongen et al., 2016) to obtain cell type level eQTL with only the SNPs tested for sQTL analysis. The same covariables were used with “--window 1e7” parameter setting. The gene expression matrix was obtained by the same process as the transcript matrix. We only kept eQTLs with q-value < 0.1.
2.9 RBP binding site
We obtained RNA binding protein (RBP) site from POSTAR3 (http://111.198.139.65/index.html) (Zhao et al., 2022) with only the “experimental results”. We obtained RBP overlapping SNP if the SNP is located within the RBP binding region. We only considered an RBP if it is expressed in the cell type (log2 (exp+1) > 0.25 and fraction of cells expressing > 0.1).
2.10 GWAS analysis
To analyze the association between GWAS related to pancreas disease and SNPs from sQTL, we collected “type 2 diabetes mellitus” (EFO ID: MONDO_0005148), “type 1 diabetes mellitus” (EFO ID: MONDO_0005147), and “pancreatitis” (EFO ID: EFO_0000278) from GWAS catalog (https://www.ebi.ac.uk/gwas/) (Sollis et al., 2023). We obtained “glycemic traits” from Ji et al. (). Since GWAS SNPs were in hg38 coordinates, we converted hg38 genome coordinate into hg19 by LiftOver (https://liftover.broadinstitute.org/) (Hinrichs et al., 2006). Unknown alleles were subsequently removed. We obtained LD > 0.6 SNPs from each GWAS by SNiPA (setting: GRCh37.1000 genomes, Phase 3 v5, European, Ensembl 87) (https://snipa.org/snipa3/) ().
2.11 Splicing type analysis
We used the “splicing events classification” pipeline from sQTLseekeR (v2.1) (Monlong et al., 2014).
2.12 Pathway analysis
We used EnrichR (Xie et al., 2021) with KEGG (Kanehisa, 2019), NCI-Nature (Schaefer et al., 2009), and Gene ontology (Gene Ontology et al., 2023). We required FDR < 0.1 for a pathway to be considered significant.
3 Results
3.1 Overview of the single-cell sQTL analysis
We collected eight full-length scRNA-seq pancreatic islet datasets (Segerstolpe et al., 2016; Wang et al., 2016; Xin et al., 2016; ; ; Lawlor et al., 2017; ; Marquez-Curtis et al., 2022; Supplementary Table S1), comprising of 4,841 cells from 54 healthy individuals and 2,574 cells from 21 type 2 diabetes (T2D) patients divided amongst six cell types (acinar, alpha, beta, gamma, delta, and ductal cells) (Supplementary Table S1). We used the CTAT pipeline (), which is based on the STAR aligner and the GATK-best practice variant calling pipeline. Importantly, it has been shown to perform well in a benchmark for genotype inference from scRNA-seq (Liu et al., 2019). To infer the genotype we combined all the mutation profiles from each cell within each individual. Since these profiles may contain somatic mutations, we only keep the germline variants predicted by Strelka2 (Kim et al., 2018). We further filtered out the SNPs with minor allele frequency (MAF) < 0.1. As a result, we could obtain ∼24,000 SNPs per individual, similar to the ∼35,000 SNPs that would be expected from exome sequencing (Marian, 2014) (Supplementary Tables S1, S2) The lower number of SNPs is expected since we are restricted to the genes that are expressed in our sample. For transcript quantification we used RSEM (Li and Dewey, 2011) since it has been reported as the best performing tool for isoform level expression quantification with full-length scRNA-seq data (Westoby et al., 2018). Prior to clustering we ran Harmony (Korsunsky et al., 2019) to correct the batch effect from healthy donors or T2D patients, respectively. Then, we removed 416 cells from the healthy and 283 cells from the T2D group that had been aligned to different cell types from the original paper, assuming they were batch outliers (Supplementary Figure S2; Supplementary Table S3). For cell type level sQTL analysis, we created pseudobulks for each cell type following the same strategy used in single-cell eQTL analysis (; Yazar et al., 2022). However, due to the small number of cells, we only analyzed alpha, beta, delta, and ductal cells (Figure 1A). We also created a pseudobulk of pancreatic islets called “collapsed” to mimic bulk pancreatic islet data (Supplementary Table S4).
FIGURE 1
After preprocessing, we used sQTLseeker2 (Garrido-Martin et al., 2021) to obtain cell type level sQTLs (Figure 1B). We followed the recommendations by the authors of sQTLseeker2 and filtered out the lowly expressed transcripts. We applied a linear model to adjust the expression levels, using age and gender as covariates. A 50 kbp window centered around each splice junction was used to identify candidate variants since previous studies have shown that most sQTLs occur within this distance, while larger windows reduce power due to multiple-hypothesis correction (Garrido-Martin et al., 2021; ; ). The number of variants tested was further reduced by merging those that had a linkage disequilibrium > 0.8. Not all cell types were annotated in all datasets, and for each cell type we had at least 15 individuals per genotype.
Our analysis yielded 136 sQTLs from alpha cells, 92 from beta cells, while only 1 or 2 sQTLs were significant in other cell types at an FDR of 0.1. For the collapsed dataset we identified 268 sQTLs (Supplementary Tables S5, S6). Previous analyses using bulk data in the GTEX database (Shan et al., 2019) have shown that the number of significant eQTLs is strongly correlated with the sample size (Pearson correlation coefficient: 0.86) (Supplementary Figure S1A). We also observed a similar pattern from our data by down-sampling (Supplementary Figure S1B). We assumed a similar relationship exists for sQTLs, and extrapolating from a brain study, the expected number of sQTLs from 54 samples is ∼300 (63 samples: 387, 124 samples: 849) (Zhang et al., 2020), similar to the 268 sQTLs that we observed in the collapsed data. Similarly, a study of sQTLs using bulk RNAseq from pancreatic islets (), reported 4,858 sQTLs from 400 samples using an FDR of 0.01. Assuming that the number of sQTLs scales linearly with the sample size, then 268 sQTLs for 54 samples is an expected result. Taken together, these comparisons show that the number of sQTLs from the “collapsed” data was within the expected range. For T2D patients, due to the sample size, we could not get any suitable number of significant results (Supplementary Table S5).
We also note that the effect sizes are often substantial. For both alpha and beta cells, > 92% of sQTLs result in a change in a maximum difference of relative abundance > 0.2 which is considered a large effect size in the sQTLseeker2 paper (Garrido-Martin et al., 2021). Comparing the distribution of sQTLs across cell types we found a substantial number of sQTLs and genes showing the alteration in the splicing affected by these variants (sgenes) that do not overlap with the “collapsed” result (Figure 2A). This result is consistent with single cell eQTL studies which have shown that bulk transcriptomes can mask the signal from individual cell types. For both alpha and beta cells, the overlap was only around 15%–20%, indicating that AS in different cell types is regulated by different SNPs (Figure 2B). At the same time, the splicing of a given gene was also affected differently between different cell types with only 15.8% of sgenes found in both alpha and beta cells, suggesting that the set of sQTLs affecting a given pathway will differ between cell types.
FIGURE 2
Interestingly, the expression levels of sgenes are typically not higher in the cell type where they are active. We performed differentially expressed gene (DEG) analysis between alpha and beta cells to obtain highly expressed genes for each cell type. The 32% for alpha cells and 20% for beta cells is the ratio of the DE genes among sgenes. However, we only found 1.7% of highly expressed genes in the alpha cell (27/1,626), and 0.32% of highly expressed genes in the beta cell (9/2,780) were cell type-specific sGenes. Hence, this implies that DE genes are distinct from sgenes.
3.2 Functional analysis of sQTL
Next, we characterized the sQTLs’ genomic locations (Figure 3A; Supplementary Table S7; ), and this revealed that most sQTLs were located outside of exonic regions, with only 26/137 alpha cell variants and 15/92 beta cell variants overlapping exons. Comparing each region to the total number of detected variants suggested that SNPs from intergenic regions were not enriched above chance level since the odd ratios were close to 1. A previous study of sQTLs in the pancreatic islet showed that 5′ and 3′ UTRs were the most enriched regions (). Although we could not find as many sQTLs from the 5′ UTR as in the previous report, we observed substantial numbers of sQTLs in 3′ UTR (odd ratios, alpha: 2.14, beta:1.77, collapsed: 1.57). We also observed a higher number of sQTLs in intronic regions which was not surprising given its key role in splicing (Jo and Choi, 2015), but the odd ratios were surprisingly low (alpha: 0.33, beta: 0.7, collapsed: 0.57). Although the splice site has been shown to be enriched for sQTLs (Garrido-Martin et al., 2021; ; Qi et al., 2022), we did not observe this in our study. One explanation could be that the splice junction can only be observed from unspliced transcripts. We also cannot rule out that the discrepancy from previous studies is due to variant calling from the transcripts. Second, we investigated whether sQTLs are associated with eQTLs (Supplementary Table S8). A previous bulk pancreatic islet study only reported a weak association between sQTLs and eQTLs (∼25%) (), suggesting that regulation of gene expression levels and splicing are largely independent. Similarly, when we compared cell type level sQTLs and eQTLs, we found a low overlap (<27%) (Figure 3B), suggesting that gene expression and splicing regulation are largely independent also at the cell type level. Third, we hypothesize that sQTLs would be enriched at RNA binding protein (RBP) (Zhang et al., 2020; Garrido-Martin et al., 2021; Qi et al., 2022) sites. We analyzed whether sQTLs are enriched at RBP sites using the POSTAR3 database (Zhao et al., 2022), and as expected odd ratios were higher than 1 for alpha, beta, and “collapsed” (Figure 3C; Supplementary Table S9). Fourth, since many sQTLs have been associated with disease (Zhang et al., 2020; Garrido-Martin et al., 2021; ; ; Qi et al., 2022), we investigated overlap with GWAS variants related to type 1 diabetes (T1D), type 2 diabetes (T2D), and glycemic traits () within LD > 0.6 SNPs using SNiPA (). However, we could not observe any GWAS association for our sQTL results, and we conjecture that this is due to the small number of variants identified.
FIGURE 3
3.3 Cell type-specific splicing in sQTL
Next, we characterized the splicing events involved in sQTL using the definitions from sQTLseeker (Monlong et al., 2014). Most of the splicing events were complex 5′ events, complex 3′ events, and skipped exon events, except unknown events (Monlong et al., 2014; Figure 4A; Supplementary Table S10). This finding was consistent for alpha, beta, and “collapsed”. Compared to previous reports, we observed a similar number of complex 5′ and 3′ events, but we observed more skipped exons (Monlong et al., 2014). We also carried out functional analysis for sgenes using EnrichR (Xie et al., 2021) to quantify the overlap with KEGG (Kanehisa, 2019), NIC-Nature (2016) (Schaefer et al., 2009), and Gene ontology (2023) (Gene Ontology et al., 2023; Figure 4B; Supplementary Table S11). A previous bulk pancreatic islet sQTL study () showed enrichment for terms related to virus infection (herpes simplex virus 1 infection and viral myocarditis) from the KEGG pathway. Our “collapsed” result was also enriched in bacterial response (yersinia infection and shigellosis), indicating an association with the inflammation pathway. Previously, it has been reported that the association between inflammation and pancreatic islets such as microbial DNA mediated pancreatic islet inflammation in obesity (Gao et al., 2022) or bacterial impact on T1D (). Hence, this result indicates that one of the mechanisms through which the genotype impacts obesity or T1D is via changed splicing pattern in response to inflammation. However, only a few pathways were enriched in the “collapsed”, indicating a dilution effect by averaging the signals from different cell types.
FIGURE 4
On the other hand, there were several enriched pathways from both alpha and beta cells that were found in only one cell type (Figure 4B; Supplementary Figure S3). One interesting finding was that the alpha cells are enriched for inflammation and fatty acid pathways, while the beta cells were enriched for the mTORC1 pathway. Both the association between the alpha cell and fatty acid pathway (), and the association between the beta cell and mTOR signaling (; ; ) were reported previously from a bulk study (), indicating cell type specific regulation. The alpha cells were also enriched in TGF-beta signaling, which was surprising since a tight association between beta cells and TGF-beta signaling has been reported previously (Lee et al., 2021; Wang et al., 2022). This result implies an additional functional association of TGF-beta signaling in alpha cells.
Another interesting pathway was the cell morphology pathway mediated by CDC42 which was enriched in both alpha and beta cells. It has been reported that CDC42 is related to cytoskeleton and cellular structure (; Watson et al., 2017; Pichaud et al., 2019). This suggests that specific variants can alter pancreatic islet morphology via splicing alterations in a cell type specific manner. In addition, it has been reported that the morphology of the pancreatic islet is associated with diabetes (; Okajima et al., 2022). Especially, CDC42 was reported as a regulator for not only diabetes, but also insulin secretion and development (Martinelli et al., 2018; Veluthakal et al., 2018; Huang et al., 2019; He et al., 2020). This implies a further association between splicing alteration via different genotypes and pancreatic islet diseases like diabetes. One example of cell type specific splicing regulation of a pathway is given by BCAR1 which was associated with CDC42 only in alpha cells (Figures 4B, 5A; Supplementary Figure S4A). Reassuringly, it was previously reported to affect cell differentiation via CDC42 (). By contrast, RAC1 was associated with CDC42 only in beta cells (Figures 4B, 5B; Supplementary Figure S4B), consistent with previous reports linking it to pancreatic islet morphogenesis via beta cells (Greiner et al., 2009). However, none of these three genes were observed in previous bulk pancreatic islet sQTLs (). We investigated CDC42 more carefully and found four SNPs associated with this sQTL from alpha and beta cells. Interestingly, even though they influence the same gene, the perturbation in each cell type was different resulting in varying ratios between isoforms depending on the combination of genetic variants (Figure 5C). This result indicates an additional layer of heterogeneity for sQTLs, which will lead to more diversification in the transcriptome in a cell type specific manner.
FIGURE 5
4 Discussion
We carried out a single-cell sQTL analysis of the pancreatic islet, identifying cell type level variants regulating AS. As expected, the single-cell resolution allowed us to identify variants that were not revealed when analyzing bulk samples. Interestingly, we also captured cell type level regulation of splicing event between different genotypes. Reassuringly, we found a coherent regulation of splicing events by sgenes in each cell type which were often related to diabetes, or a phenotype relevant to pancreatic islets (Monlong et al., 2014; Garrido-Martin et al., 2021; ; ; Yamaguchi et al., 2022). This indicates that function or phenotype of the pancreatic islet can affect individuals by different splicing events across cell types. In particular, we found a cell type specific sQTL for CDC42, a gene which had previously been related to pancreatic islet morphology, which has not been reported elsewhere.
However, our study has a several limitations. The cohort only included 54 healthy donors, and not all cell types were represented in all donors. The poor representation of some cell types resulted in few significant sQTLs and poor overlap with pancreatic islet related GWAS. Furthermore, we used a relatively loose threshold compared to other more well powered studies. Also, usually nominal p value and empirical p value by permutation test were used together for eQTL or sQTL analysis. Due to lack of the sample size we could not obtain a suitable number of sQTLs by empirical p value. Another limitation is our restricted genotyping, and we can only call variants overlapping expressed transcripts. Due to the low number of immature mRNAs, the number of variants from intronic regions or splicing sites is limited. Moreover, we did not utilize the population structure to regress out such effects or gene expression distribution (Price et al., 2006; Garrido-Martin et al., 2021). However, these regression models have been reported to result in false positives (), so given our small sample size this may not be such a severe drawback. Although it is possible to impute SNPs (Marchini and Howie, 2010), we did not use this approach since we cannot guarantee the reliability when only using SNPs from the transcriptome.
Nevertheless, comparison to bulk sQTL and single cell eQTL analyses suggest that we are able to detect the expected number of variants. Hence, by increasing the sample size it will most likely be possible to identify more splice affecting variants. Even though the scRNAseq protocols used in this study are designed to capture mature mRNAs, they will frequently profile non-exonic regions derived from nascent transcripts (Hangauer et al., 2013; ). Here, this was an advantage, as it allowed us to discover additional sQTL variants located in these regions. By using WGS for the genotyping, one would most likely obtain even more cell type level sQTLs. Recently, long read technologies have been combined with single cell RNAseq technologies, and this is likely to be a better and more cost effective approach for obtaining information about isoform usage (Rhoads and Au, 2015; Uapinyoying et al., 2020; Huang et al., 2021). In particular, since we showed that variant calling via short full-length is feasible, this provides further evidence that long read sequencing will be useful for sQTL studies since these technologies generally perform better for variant calling (Shiau et al., 2023; Kim et al., 2024; Penter et al., 2024).
The main limiting factor of our study is not the number of cells, but the number of donors. When we randomly downsampled the cells for each donor, there was almost no difference in the number of sQTLs, indicating a marginal effect by cell count per donor (Supplementary Figure S5). Another benefit from a larger and more diverse set of donors is that it will be possible to reduce errors through population-structure level batch correction and batch correction at a donor level. In addition, it will allow us to use stricter thresholds for the analysis, similar to what has been used in single-cell eQTL or bulk sQTL studies to improve the results.
Larger single cell sQTL studies are likely to further demonstrate that regulation of expression and splicing are largely disconnected at the genetic level. We conjecture that this will hold true not just for the pancreatic islets, but for other conditions and tissues more generally. The increased number of variants associated with isoform choice will most likely facilitate the interpretation of the many GWAS variants where we lack a mechanistic understanding of how they influence the phenotype. Studies of disease conditions, such as T2D, will be particularly helpful for interpreting GWAS variants, and they are likely to further our understanding of the origins and role of cellular heterogeneity in health and disease (). Taken together, our results provide a proof of concept that it is possible to identify cell type level sQTLs, and we envision that better powered studies will allow us to further uncover the genetic regulation of splicing.
Statements
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.
Author contributions
J-WC: Writing – original draft, Software, Writing – review and editing, Resources, Investigation, Formal Analysis, Validation, Visualization, Methodology, Data curation. JC: Writing – review and editing, Investigation, Formal Analysis, Validation, Visualization, Methodology. MH: Writing – review and editing, Conceptualization, Funding acquisition, Supervision, Writing – original draft, Project administration.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. This work was funded by the Evergrande Center and the Helmsley Foundation (2008-04050).
Acknowledgments
We would like to thank Ruben Chazarra-Gil and Marta Melé for helpful feedback on the manuscript.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbinf.2025.1657895/full#supplementary-material
References
1
AbdellatifA. M.Jensen SmithH.HarmsR. Z.SarvetnickN. E. (2019). Human islet response to selected type 1 diabetes-associated bacteria: a transcriptome-based study. Front. Immunol.10, 2623. 10.3389/fimmu.2019.02623
2
AgostiniF.ZagalakJ.AttigJ.UleJ.LuscombeN. M. (2021). Intergenic RNA mainly derives from nascent transcripts of known genes. Genome Biol.22 (1), 136. 10.1186/s13059-021-02350-x
3
ArmourS. L.FruehA.ChibalinaM. V.DouH.Argemi-MuntadasL.HamiltonA.et al (2023). Glucose controls glucagon secretion by regulating fatty acid oxidation in pancreatic alpha-cells. Diabetes72 (10), 1446–1459. 10.2337/db23-0056
4
ArnoldM.RafflerJ.PfeuferA.SuhreK.KastenmullerG. (2015). SNiPA: an interactive, genetic variant-centered annotation browser. Bioinformatics31 (8), 1334–1336. 10.1093/bioinformatics/btu779
5
Arzalluz-LuqueA.ConesaA. (2018). Single-cell RNAseq for the study of isoforms-how is that possible?Genome Biol.19 (1), 110. 10.1186/s13059-018-1496-z
6
AsaharaS. I.InoueH.WatanabeH.KidoY. (2022). Roles of mTOR in the regulation of pancreatic beta-cell mass and insulin secretion. Biomolecules12 (5), 614. 10.3390/biom12050614
7
AtlaG.Bonas-GuarchS.Cuenca-ArduraM.BeucherA.CrouchD. J. M.Garcia-HurtadoJ.et al (2022). Genetic regulation of RNA splicing in human pancreatic islets. Genome Biol.23 (1), 196. 10.1186/s13059-022-02757-0
8
BarrieE. S.SmithR. M.SanfordJ. C.SadeeW. (2012). mRNA transcript diversity creates new opportunities for pharmacological intervention. Mol. Pharmacol.81 (5), 620–630. 10.1124/mol.111.076604
9
BauerJ.GomboczE.SchulzH.HauslageJ.GrimmD. (2021). Interaction network provides clues on the role of BCAR1 in cellular response to changes in gravity. Computation9 (8), 81. 10.3390/computation9080081
10
BenningerR. K. P.KravetsV. (2022). The physiological role of beta-cell heterogeneity in pancreatic islet function. Nat. Rev. Endocrinol.18 (1), 9–22. 10.1038/s41574-021-00568-0
11
BlackD. L. (2000). Protein diversity from alternative splicing: a challenge for bioinformatics and post-genome biology. Cell103 (3), 367–370. 10.1016/s0092-8674(00)00128-8
12
Blandino-RosanoM.ChenA. Y.ScheysJ. O.AlejandroE. U.GouldA. P.TaranukhaT.et al (2012). mTORC1 signaling and regulation of pancreatic beta-cell mass. Cell Cycle11 (10), 1892–1902. 10.4161/cc.20036
13
Blandino-RosanoM.BarbaressoR.Jimenez-PalomaresM.BozadjievaN.Werneck-de-CastroJ. P.HatanakaM.et al (2017). Loss of mTORC1 signalling impairs beta-cell homeostasis and insulin processing. Nat. Commun.8, 16014. 10.1038/ncomms16014
14
Bonner-WeirS.SullivanB. A.WeirG. C. (2015). Human islet morphology revisited: human and rodent islets are not so different after all. J. Histochem Cytochem63 (8), 604–612. 10.1369/0022155415570969
15
BrotmanS. M.RaulersonC. K.VadlamudiS.CurrinK. W.ShenQ.ParsonsV. A.et al (2022). Subcutaneous adipose tissue splice quantitative trait loci reveal differences in isoform usage associated with cardiometabolic traits. Am. J. Hum. Genet.109 (1), 66–80. 10.1016/j.ajhg.2021.11.019
16
BrouardJ. S.BissonnetteN. (2022). Variant calling from RNA-seq data using the GATK joint genotyping workflow. Methods Mol. Biol.2493, 205–233. 10.1007/978-1-0716-2293-3_13
17
ButlerA.HoffmanP.SmibertP.PapalexiE.SatijaR. (2018). Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat. Biotechnol.36 (5), 411–420. 10.1038/nbt.4096
18
Camunas-SolerJ.DaiX. Q.HangY.BautistaA.LyonJ.SuzukiK.et al (2020). Patch-seq links single-cell transcriptomes to human islet dysfunction in diabetes. Cell Metab.31 (5), 1017–1031.e4. 10.1016/j.cmet.2020.04.005
19
CarithersL. J.MooreH. M. (2015). The genotype-tissue expression (GTEx) Project. Biopreserv Biobank13 (5), 307–308. 10.1089/bio.2015.29031.hmm
20
ChenJ.SpracklenC. N.MarenneG.VarshneyA.CorbinL. J.LuanJ.et al (2021). The trans-ancestral genomic architecture of glycemic traits. Nat. Genet.53 (6), 840–860. 10.1038/s41588-021-00852-9
21
ClawsonH.LeeB. T.RaneyB. J.BarberG. P.CasperJ.DiekhansM.et al (2023). GenArk: towards a million UCSC genome browsers. Genome Biol.24 (1), 217. 10.1186/s13059-023-03057-x
22
GTEx Consortium (2017). Genetic effects on gene expression across human tissues. Nature550 (7675), 204–213. 10.1038/nature24277
23
CuomoA. S. E.AlvariG.AzodiC. B.CarthyD. J.BonderM. J. (2021). Optimizing expression quantitative trait locus mapping workflows for single-cell studies. Genome Biol.22 (1), 188. 10.1186/s13059-021-02407-x
24
DahlA.GuillemotV.MeffordJ.AschardH.ZaitlenN. (2019). Adjusting for principal components of molecular phenotypes induces replicating false positives. Genetics211 (4), 1179–1189. 10.1534/genetics.118.301768
25
DanecekP.BonfieldJ. K.LiddleJ.MarshallJ.OhanV.PollardM. O.et al (2021). Twelve years of SAMtools and BCFtools. Gigascience10 (2), giab008. 10.1093/gigascience/giab008
26
DorajooR.AliY.TayV. S. Y.KangJ.SamyduraiS.LiuJ.et al (2017). Single-cell transcriptomics of East-Asian pancreatic islets cells. Sci. Rep.7 (1), 5024. 10.1038/s41598-017-05266-4
27
EngeM.ArdaH. E.MignardiM.BeausangJ.BottinoR.KimS. K.et al (2017). Single-cell analysis of human pancreas reveals transcriptional signatures of aging and somatic mutation patterns. Cell171 (2), 321–330.e14. 10.1016/j.cell.2017.09.004
28
Escobar-HoyosL.KnorrK.Abdel-WahabO. (2019). Aberrant RNA splicing in cancer. Annu. Rev. Cancer Biol.3 (1), 167–185. 10.1146/annurev-cancerbio-030617-050407
29
FangalV. D. (2020). CTAT mutations: a machine learning based RNA-seq variant calling pipeline incorporating variant annotation, prioritization, and visualization. Master thesis. Harvard Extension School. Available online at: https://dash.harvard.edu/entities/publication/48734d69-35ea-415a-9775-95db9191091a.
30
FengD.XieJ. (2013). Aberrant splicing in neurological diseases. Wiley Interdiscip. Rev. RNA4 (6), 631–649. 10.1002/wrna.1184
31
FilipenkoN. R.AttwellS.RoskelleyC.DedharS. (2005). Integrin-linked kinase activity regulates Rac- and Cdc42-mediated actin cytoskeleton reorganization via alpha-PIX. Oncogene24 (38), 5837–5849. 10.1038/sj.onc.1208737
32
GaoH.LuoZ.JiY.TangK.JinZ.LyC.et al (2022). Accumulation of microbial DNAs promotes to islet inflammation and beta cell abnormalities in obesity in mice. Nat. Commun.13 (1), 565. 10.1038/s41467-022-28239-2
33
Garrido-MartinD.BorsariB.CalvoM.ReverterF.GuigoR. (2021). Identification and analysis of splicing quantitative trait loci across multiple tissues in the human genome. Nat. Commun.12 (1), 727. 10.1038/s41467-020-20578-2
34
Gene OntologyC.AleksanderS. A.BalhoffJ.CarbonS.CherryJ. M.DrabkinH. J.et al (2023). The gene ontology knowledgebase in 2023. Genetics224 (1), iyad031. 10.1093/genetics/iyad031
35
GreinerT. U.KesavanG.StahlbergA.SembH. (2009). Rac1 regulates pancreatic islet morphogenesis. BMC Dev. Biol.9, 2. 10.1186/1471-213X-9-2
36
HangauerM. J.VaughnI. W.McManusM. T. (2013). Pervasive transcription of the human genome produces thousands of previously unidentified long intergenic noncoding RNAs. PLoS Genet.9 (6), e1003569. 10.1371/journal.pgen.1003569
37
HeX. Q.WangN.ZhaoJ. J.WangD.WangC. J.XieL.et al (2020). Specific deletion of CDC42 in pancreatic beta cells attenuates glucose-induced insulin expression and secretion in mice. Mol. Cell Endocrinol.518, 111004. 10.1016/j.mce.2020.111004
38
HinrichsA. S.KarolchikD.BaertschR.BarberG. P.BejeranoG.ClawsonH.et al (2006). The UCSC genome browser database: update 2006. Nucleic Acids Res.34 (Database issue), D590–D598. 10.1093/nar/gkj144
39
HuangQ. Y.LaiX. N.QianX. L.LvL. C.LiJ.DuanJ.et al (2019). Cdc42: a novel regulator of insulin secretion and diabetes-associated diseases. Int. J. Mol. Sci.20 (1), 179. 10.3390/ijms20010179
40
HuangK. K.HuangJ.WuJ. K. L.LeeM.TayS. T.KumarV.et al (2021). Long-read transcriptome sequencing reveals abundant promoter diversity in distinct molecular subtypes of gastric cancer. Genome Biol.22 (1), 44. 10.1186/s13059-021-02261-x
41
JoB. S.ChoiS. S. (2015). Introns: the functional benefits of introns in genomes. Genomics Inf.13 (4), 112–118. 10.5808/GI.2015.13.4.112
42
KanehisaM. (2019). Toward understanding the origin and evolution of cellular organisms. Protein Sci.28 (11), 1947–1951. 10.1002/pro.3715
43
KerimovN.TambetsR.HayhurstJ. D.RahuI.KolbergP.RaudvereU.et al (2023). eQTL Catalogue 2023: new datasets, X chromosome QTLs, and improved detection and visualisation of transcript-level QTLs. PLoS Genet.19 (9), e1010932. 10.1371/journal.pgen.1010932
44
KimS.SchefflerK.HalpernA. L.BekritskyM. A.NohE.KallbergM.et al (2018). Strelka2: fast and accurate calling of germline and somatic variants. Nat. Methods15 (8), 591–594. 10.1038/s41592-018-0051-x
45
KimH. S.GrimesS. M.ChenT.SatheA.LauB. T.HwangG. H.et al (2024). Direct measurement of engineered cancer mutations and their transcriptional phenotypes in single cells. Nat. Biotechnol.42 (8), 1254–1262. 10.1038/s41587-023-01949-8
46
KorsunskyI.MillardN.FanJ.SlowikowskiK.ZhangF.WeiK.et al (2019). Fast, sensitive and accurate integration of single-cell data with Harmony. Nat. Methods16 (12), 1289–1296. 10.1038/s41592-019-0619-0
47
LawlorN.GeorgeJ.BolisettyM.KursaweR.SunL.SivakamasundariV.et al (2017). Single-cell transcriptomes identify human islet cell signatures and reveal cell-type-specific expression changes in type 2 diabetes. Genome Res.27 (2), 208–222. 10.1101/gr.212720.116
48
LawrenceM.HuberW.PagesH.AboyounP.CarlsonM.GentlemanR.et al (2013). Software for computing and annotating genomic ranges. PLoS Comput. Biol.9 (8), e1003118. 10.1371/journal.pcbi.1003118
49
LeeJ. H.LeeJ. H.RaneS. G. (2021). TGF-beta signaling in pancreatic islet beta cell development and function. Endocrinology162 (3), bqaa233. 10.1210/endocr/bqaa233
50
LiB.DeweyC. N. (2011). RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinforma.12, 323. 10.1186/1471-2105-12-323
51
LiuF.ZhangY.ZhangL.LiZ.FangQ.GaoR.et al (2019). Systematic comparative analysis of single-nucleotide variant detection methods from single-cell RNA sequencing data. Genome Biol.20 (1), 242. 10.1186/s13059-019-1863-4
52
Lopez-BigasN.AuditB.OuzounisC.ParraG.GuigoR. (2005). Are splicing mutations the most frequent cause of hereditary disease?FEBS Lett.579 (9), 1900–1903. 10.1016/j.febslet.2005.02.047
53
ManolioT. A. (2010). Genomewide association studies and assessment of the risk of disease. N. Engl. J. Med.363 (2), 166–176. 10.1056/NEJMra0905980
54
MarchiniJ.HowieB. (2010). Genotype imputation for genome-wide association studies. Nat. Rev. Genet.11 (7), 499–511. 10.1038/nrg2796
55
MarianA. J. (2014). Sequencing your genome: what does it mean?Methodist Debakey Cardiovasc J.10 (1), 3–6. 10.14797/mdcj-10-1-3
56
Marquez-CurtisL. A.DaiX. Q.HangY.LamJ. Y.LyonJ.Manning FoxJ. E.et al (2022). Cryopreservation and post-thaw characterization of dissociated human islet cells. PLoS One17 (1), e0263005. 10.1371/journal.pone.0263005
57
MartinelliS.KrumbachO. H. F.PantaleoniF.CoppolaS.AminE.PannoneL.et al (2018). Functional dysregulation of CDC42 causes diverse developmental phenotypes. Am. J. Hum. Genet.102 (2), 309–320. 10.1016/j.ajhg.2017.12.015
58
MonlongJ.CalvoM.FerreiraP. G.GuigoR. (2014). Identification of genetic variants associated with alternative splicing using sQTLseekeR. Nat. Commun.5, 4698. 10.1038/ncomms5698
59
NathanA.AsgariS.IshigakiK.ValenciaC.AmariutaT.LuoY.et al (2022). Single-cell eQTL models reveal dynamic T cell state dependence of disease loci. Nature606 (7912), 120–128. 10.1038/s41586-022-04713-1
60
NeavinD.NguyenQ.DaniszewskiM. S.LiangH. H.ChiuH. S.WeeY. K.et al (2021). Single cell eQTL analysis identifies cell type-specific genetic control of gene expression in fibroblasts and reprogrammed induced pluripotent stem cells. Genome Biol.22 (1), 76. 10.1186/s13059-021-02293-3
61
OkajimaY.MatsuzakaT.MiyazakiS.MotomuraK.OhnoH.SharmaR.et al (2022). Morphological and functional adaptation of pancreatic islet blood vessels to insulin resistance is impaired in diabetic db/db mice. Biochim. Biophys. Acta Mol. Basis Dis.1868 (4), 166339. 10.1016/j.bbadis.2022.166339
62
OngenH.BuilA.BrownA. A.DermitzakisE. T.DelaneauO. (2016). Fast and efficient QTL mapper for thousands of molecular phenotypes. Bioinformatics32 (10), 1479–1485. 10.1093/bioinformatics/btv722
63
PenterL.BorjiM.NaglerA.LyuH.LuW. S.CieriN.et al (2024). Integrative genotyping of cancer and immune phenotypes by long-read sequencing. Nat. Commun.15 (1), 32. 10.1038/s41467-023-44137-7
64
PicardiE.D'ErchiaA. M.Lo GiudiceC.PesoleG. (2017). REDIportal: a comprehensive database of A-to-I RNA editing events in humans. Nucleic Acids Res.45 (D1), D750–D757. 10.1093/nar/gkw767
65
PichaudF.WaltherR. F.Nunes de AlmeidaF. (2019). Regulation of Cdc42 and its effectors in epithelial morphogenesis. J. Cell Sci.132 (10), jcs217869. 10.1242/jcs.217869
66
PriceA. L.PattersonN. J.PlengeR. M.WeinblattM. E.ShadickN. A.ReichD. (2006). Principal components analysis corrects for stratification in genome-wide association studies. Nat. Genet.38 (8), 904–909. 10.1038/ng1847
67
QiT.WuY.FangH.ZhangF.LiuS.ZengJ.et al (2022). Genetic control of RNA splicing and its distinct role in complex trait variation. Nat. Genet.54 (9), 1355–1363. 10.1038/s41588-022-01154-4
68
RamaswamiG.LiJ. B. (2014). RADAR: a rigorously annotated database of A-to-I RNA editing. Nucleic Acids Res.42 (Database issue), D109–D113. 10.1093/nar/gkt996
69
RhoadsA.AuK. F. (2015). PacBio sequencing and its applications. Genomics Proteomics Bioinforma.13 (5), 278–289. 10.1016/j.gpb.2015.08.002
70
SchaeferC. F.AnthonyK.KrupaS.BuchoffJ.DayM.HannayT.et al (2009). PID: the pathway interaction database. Nucleic Acids Res.37 (Database issue), D674–D679. 10.1093/nar/gkn653
71
SegerstolpeA.PalasantzaA.EliassonP.AnderssonE. M.AndreassonA. C.SunX.et al (2016). Single-cell transcriptome profiling of human pancreatic islets in health and type 2 diabetes. Cell Metab.24 (4), 593–607. 10.1016/j.cmet.2016.08.020
72
ShanN.WangZ.HouL. (2019). Identification of trans-eQTLs using mediation analysis with multiple mediators. BMC Bioinforma.20 (Suppl. 3), 126. 10.1186/s12859-019-2651-6
73
ShiauC. K.LuL.KieserR.FukumuraK.PanT.LinH. Y.et al (2023). High throughput single cell long-read sequencing analyses of same-cell genotypes and phenotypes in human tumors. Nat. Commun.14 (1), 4124. 10.1038/s41467-023-39813-7
74
SollisE.MosakuA.AbidA.BunielloA.CerezoM.GilL.et al (2023). The NHGRI-EBI GWAS Catalog: knowledgebase and deposition resource. Nucleic Acids Res.51 (D1), D977–D985. 10.1093/nar/gkac1010
75
TianC.ZhangY.TongY.KockK. H.SimD. Y.LiuF.et al (2024). Single-cell RNA sequencing of peripheral blood links cell-type-specific regulation of splicing to autoimmune and inflammatory diseases. Nat. Genet.56 (12), 2739–2752. 10.1038/s41588-024-02019-8
76
UapinyoyingP.GoecksJ.KnoblachS. M.PanchapakesanK.BonnemannC. G.PartridgeT. A.et al (2020). A long-read RNA-seq approach to identify novel transcripts of very large genes. Genome Res.30 (6), 885–897. 10.1101/gr.259903.119
77
van der WijstM.de VriesD. H.GrootH. E.TrynkaG.HonC. C.BonderM. J.et al (2020). The single-cell eQTLGen consortium. Elife9, e52155. 10.7554/eLife.52155
78
VeluthakalR.ChepurnyO. G.LeechC. A.SchwedeF.HolzG. G.ThurmondD. C. (2018). Restoration of glucose-stimulated Cdc42-Pak1 activation and insulin secretion by a selective epac activator in type 2 diabetic human islets. Diabetes67 (10), 1999–2011. 10.2337/db17-1174
79
WangB. D.LeeN. H. (2018). Aberrant RNA splicing in cancer and drug resistance. Cancers (Basel)10 (11), 458. 10.3390/cancers10110458
80
WangY. J.SchugJ.WonK. J.LiuC.NajiA.AvrahamiD.et al (2016). Single-cell transcriptomics of the human endocrine pancreas. Diabetes65 (10), 3028–3038. 10.2337/db16-0405
81
WangH. L.WangL.ZhaoC. Y.LanH. Y. (2022). Role of TGF-beta signaling in beta cell proliferation and function in diabetes. Biomolecules12 (3), 373. 10.3390/biom12030373
82
WatsonJ. R.OwenD.MottH. R. (2017). Cdc42 in actin dynamics: an ordered pathway governed by complex equilibria and directional effector handover. Small GTPases8 (4), 237–244. 10.1080/21541248.2016.1215657
83
WestobyJ.HerreraM. S.Ferguson-SmithA. C.HembergM. (2018). Simulation-based benchmarking of isoform quantification in single-cell RNA-seq. Genome Biol.19 (1), 191. 10.1186/s13059-018-1571-5
84
WolfeM. (2022). genomicMateSelectR: genomic mate selection.
85
XieZ.BaileyA.KuleshovM. V.ClarkeD. J. B.EvangelistaJ. E.JenkinsS. L.et al (2021). Gene set knowledge discovery with enrichr. Curr. Protoc.1 (3), e90. 10.1002/cpz1.90
86
XinY.KimJ.OkamotoH.NiM.WeiY.AdlerC.et al (2016). RNA sequencing of single human islet cells reveals type 2 diabetes genes. Cell Metab.24 (4), 608–615. 10.1016/j.cmet.2016.08.018
87
YamaguchiK.IshigakiK.SuzukiA.TsuchidaY.TsuchiyaH.SumitomoS.et al (2022). Splicing QTL analysis focusing on coding sequences reveals mechanisms for disease susceptibility loci. Nat. Commun.13 (1), 4659. 10.1038/s41467-022-32358-1
88
YazarS.Alquicira-HernandezJ.WingK.SenabouthA.GordonM. G.AndersenS.et al (2022). Single-cell eQTL mapping identifies cell type-specific genetic control of autoimmune disease. Science376 (6589), eabf3041. 10.1126/science.abf3041
89
ZhangY.YangH. T.Kadash-EdmondsonK.PanY.PanZ.DavidsonB. L.et al (2020). Regional variation of splicing QTLs in human brain. Am. J. Hum. Genet.107 (2), 196–210. 10.1016/j.ajhg.2020.06.002
90
ZhaoW.ZhangS.ZhuY.XiX.BaoP.MaZ.et al (2022). POSTAR3: an updated platform for exploring post-transcriptional regulation coordinated by RNA-binding proteins. Nucleic Acids Res.50 (D1), D287–D294. 10.1093/nar/gkab702
Summary
Keywords
alternative splicing, single-cell RNA-sequencing, pancreatic islet, splicing quantitative trait loci, CDC42
Citation
Cho J-W, Cao J and Hemberg M (2025) Single-cell splicing QTL analysis in pancreatic islets. Front. Bioinform. 5:1657895. doi: 10.3389/fbinf.2025.1657895
Received
02 July 2025
Accepted
28 August 2025
Published
10 September 2025
Volume
5 - 2025
Edited by
Bruno Andrade, INRA UMR1253 Science and Technologie du Lait and de l’œuf, France
Reviewed by
Youtao Lu, University of Pennsylvania, United States
Guanjue Xiang, Dana–Farber Cancer Institute, United States
Updates
Copyright
© 2025 Cho, Cao and Hemberg.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Martin Hemberg, mhemberg@bwh.harvard.edu
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.