Sequence data analysis methods and systems

Abstract

Disclosed herein is a sequence data analysis method, comprising: receiving single-cell RNA sequence data comprising RNA sequence reads; aligning the RNA sequence reads to a reference genome; extracting intron junctions from the aligned RNA sequence reads; determining read counts for the extracted intron junctions; and identifying splicing events using the read counts for the extracted intron junctions. In particular, the method comprises identifying splicing events from intron excision ratios across the extracted intron junctions and further performing splicing quantitative trait loci mapping.

Description

[0001] SEQUENCE DATA ANALYSIS METHODS AND SYSTEMS

[0002] TECHNICAL FIELD

[0003] The present disclosure relates bioinformatics and the analysis of sequence data. In particular, the present disclosure relates to splicing analysis.

[0004] BACKGROUND

[0005] Ribonucleic acid (RNA) alternative splicing is a crucial mediator of complex diseases. Genetic variants can affect splice site usage or splicing quantitative trait loci (sQTLs), thereby affecting the relative ratio of isoform expression. Studies have identified tens of thousands of tissue-specific sQTLs across a wide range of human tissues, yet a tissue consists of many cell types and tissue-level investigation may overlook cell typespecific genetic regulation. To address this gap, recent studies have mapped cell typespecific effects with bulk RNA-seq on cell cultures and Fluorescence-activated cell sorting (FACS) sorted primary cells, as well as with single-cell RNA-seq data (scRNA-seq). In particular, the increasing scale and comprehensive cell type coverage offered by single-cell RNA-seq has made it a popular choice for recent cell type-specific quantitative trait loci (QTL) studies. However, current single-cell sQTL studies must rely on pseudobulk analysis because existing sQTL mapping tools are designed for bulk RNA-seq data.

[0006] However, pseudobulk analyses suffers from several important aspects. First, pseudobulk aggregation miss important information about cell-to-cell variability. It is widely accepted that, within the same cell type, cell state is a key determinant of genetic regulation, and neither cell culture, FACS-sorted, or pseudobulk analysis provide information about cell state-dependent regulation. Second, single-cell expression quantitative trait locus (eQTL) studies have shown that each single cell can be treated as a data point, and proper statistical modeling of repeated single-cell measurements can boost statistical power eQTL calling. However, current pseudobulk sQTL procedure does not properly account for the repeated single-cell measurements within the same donor, which leads to a reduction in power. Third, current isoform-, exon-, and intron-based bulk level sQTL mapping assumes even RNA-seq read coverage typical for bulk RNA-seq. In contrast, scRNA-seq data often violates this assumption with 573’-biased read coverage, and it is unclear whether isoform-, exon-, and intron-based splicing phenotypes are appropriate for end-biased read coverage. In a recent scRNA-seq sQTL study, pseudobulk analysis missed intron retention events as a key component of splicing regulation.

[0007] SUMMARY

[0008] According to a first aspect of the present disclosure a sequence data analysis method is provided. The sequence data analysis method comprises: receiving single-cell RNA sequence data comprising RNA sequence reads; aligning the RNA sequence reads to a reference genome; extracting intron junctions from the aligned RNA sequence reads; determining read counts for the extracted intron junctions; and identifying splicing events using the read counts for the extracted intron junctions.

[0009] In an embodiment, identifying splicing events using the read counts for the extracted intron junctions comprises identifying splicing events from site usage ratios of detected able splice sites.

[0010] in an embodiment, the method further comprises performing splicing quantitative trait loci mapping using the splice site usage ratios.

[0011] In an embodiment, identifying splicing events using the read counts for the extracted intron junctions comprises identifying splicing events from intron excision ratios across the extracted intron junctions.

[0012] In an embodiment, the method further comprises performing splicing quantitative trait loci mapping using the intron excision ratios across the extracted intron junctions.

[0013] In an embodiment, the method further comprises identifying disease-associated splicing events by assessing colocalization between the splicing quantitative trait loci and genome-wide association study loci. In an embodiment, the disease associated splicing event is associated with blood traits, anthropometric traits and / or immune-related traits.

[0014] In an embodiment, aligning the RNA sequence reads to a reference genome is carried out in two-pass mode.

[0015] In an embodiment, the method further comprises detecting fusion genes from an identified splicing event.

[0016] In an embodiment, the method further comprises detecting a mis-splicing event associated with a disease.

[0017] In an embodiment, the disease include rare diseases that affects less than 65 out of 100 000 individuals such as spinal muscular dystrophy, Duchenne muscular dystrophy, and familial amyloid polyneuropathy.

[0018] In an embodiment, the method further comprises identifying a target gene having a mutation which leads to an increased risk of a complex (polygenic) disease, such as Graves’ disease, rheumatoid arthritis, asthma, atopic dermatitis, and systemic lupus erythematosus.

[0019] According to a second aspect of the present disclosure a sequence data processing system configured to carry out a method as set out above is provided.

[0020] According to a third aspect of the present disclosure, a computer readable medium carrying computer executable instructions which when executed on a processor cause the processor to carry out a method configured to carry out a method as set out above is provided

[0021] BRIEF DESCRIPTION OF THE DRAWINGS

[0022] In the following, embodiments of the present invention will be described as non-limiting examples with reference to the accompanying drawings in which: FIG.1 is a block diagram showing a sequence data analysis system according to an embodiment of the present invention;

[0023] FIG.2 is a flow chart showing a sequence data analysis method according to an embodiment of the present invention;

[0024] FIG.3Ato FIG.3I illustrate a method of identifying alternative splicing events according to an embodiment of the present invention;

[0025] FIG.4Ato FIG.4H illustrate investigation of disease associated alternative splicing with single-cell resolution using an embodiment of the present invention;

[0026] FIG.5A to FIG.5E show a schematic of an integrative single-cell splicing analysis and QTL caller according to an embodiment of the present invention;

[0027] FIG.6A to FIG.6F show a comparison of an integrative single-cell splicing analysis pipeline according an embodiment of the present invention and a LeafCutter pipeline;

[0028] FIG.7A to FIG.7H show an application of an integrative single-cell splicing analysis pipeline according to an embodiment of the present invention to identify cell typespecific genetic regulation of splicing in the dorsolateral prefrontal cortex;

[0029] FIG.8A to FIG.8I show an application of an integrative single-cell splicing analysis pipeline according to an embodiment of the present invention to study Intron retention sQTLs and context-biased sQTLs;

[0030] FIG.9A to FIG.9F show an application of an integrative single-cell splicing analysis pipeline according to an embodiment of the present invention to study cell-state dependent sQTLs identified in inhibitory neurons; and

[0031] FIG.10A to FIG.10H show an application of an integrative single-cell splicing analysis pipeline according to an embodiment of the present invention to study mechanisms contributing to neuronal diseases. DETAILED DESCRIPTION

[0032] The present disclosure relates to an integrative single-cell splicing analysis system and method.

[0033] FIG.1 is a block diagram showing a sequence data analysis system according to an embodiment of the present invention. The sequence data analysis system 100 is a computer system with memory that stores computer program modules which implement sequence data analysis methods according to embodiments of the present invention.

[0034] The sequence data analysis system 100 comprises a processor 110, a working memory 112, an input interface 114, an output interface 116, and program storage 120. The processor 110 may be implemented as one or more central processing unit (CPU) chips. The program storage 120 is a non-volatile storage device such as a hard disk drive which stores computer program modules. The computer program modules are loaded into the working memory 112 for execution by the processor 110. The input interface 114 is an interface which allows data to be received by the sequence data analysis system 100, for example input RNA sequence data comprising sequence reads. The input interface 114 may be a wireless network interface such as a Wi-Fi or Bluetooth interface, alternatively it may be a wired interface. The output interface 116 is an interface which allows the sequence data analysis system 100 to output results of the sequence analysis processing.

[0035] The program storage 120 stores a sequence alignment module 122, an intron junction extraction module 124, a splicing event identification module 126, a model construction module 128 and a downstream analysis module 130. The computer program modules cause the processor 110 to execute various sequence data analysis processing which is described in more detail below. The program storage 120 may be referred to in some contexts as computer readable storage media and / or non-transitory computer readable media. As depicted in FIG.1, the computer program modules are distinct modules which perform respective functions implemented by the sequence data analysis system 100. It will be appreciated that the boundaries between these modules are exemplary only, and that alternative embodiments may merge modules or impose an alternative decomposition of functionality of modules. For example, the modules discussed herein may be decomposed into sub-modules to be executed as multiple computer processes, and, optionally, on multiple computers. Moreover, alternative embodiments may combine multiple instances of a particular module or sub-module. It will also be appreciated that, while a software implementation of the computer program modules is described herein, these may alternatively be implemented as one or more hardware modules (such as field-programmable gate array(s) or application-specific integrated circuit(s)) comprising circuitry which implements equivalent functionality to that implemented in software.

[0036] Although sequence data analysis system 100 is described with reference to a computer, it should be appreciated that the sequence data analysis system 100 may be formed by two or more computers in communication with each other that collaborate to perform a task. For example, but not by way of limitation, an application may be partitioned in such away as to permit concurrent and / or parallel processing of the instructions of the application. Alternatively, the data processed by the application may be partitioned in such away as to permit concurrent and / or parallel processing of different portions of a data set by the two or more computers. In an embodiment, virtualization software may be employed by the sequence data analysis system 100 to provide the functionality of a number of servers that is not directly bound to the number of computers in the sequence data analysis system 100. In an embodiment, the functionality disclosed above may be provided by executing the application and / or applications in a cloud computing environment. Cloud computing may comprise providing computing services via a network connection using dynamically scalable computing resources. A cloud computing environment may be established by an enterprise and / or may be hired on an as-needed basis from a third party provider.

[0037] FIG.2 is a flow chart showing a sequence data analysis method according to an embodiment of the present invention. The method 200 shown in FIG.2 is carried out by the sequence data analysis system 100 shown in FIG.1

[0038] In step 202, the sequence data analysis system 100 receives single-cell RNA sequence data. The single-cell RNA sequence data comprises RNA sequence reads. The RNA sequence reads may be a single cell or a single nucleus RNA-seq dataset. In step 204, the sequence alignment module 122 is executed by the processor 110 to align the sequence reads with a reference genome.

[0039] In step 206, the intron junction extraction module 124 is executed by the processor 110 to extract intron junctions from the aligned RNA sequence reads.

[0040] In step 208, the intron junction extraction module 124 is executed by the processor 110 to determine read counts for the extracted intron junctions.

[0041] In step 210, the splicing event identification module 126 is executed by the processor 110 to identify splicing events. In some embodiments, the splicing events are identified from site usage ratios of detected able splice sites. In other embodiments, the splicing events are identified from intron excision ratios across the extracted intron junctions.

[0042] Following step 210, the model construction module 128 may be executed by the processor 110 to construct a binomial mixed model and then the downstream analysis module 130 is executed by the processor 110 to perform cell type specific sQTLs detection, cell-state dependent sQTLs identification, context-biased sQTLs discovery and / or colocalization analysis.

[0043] FIG.3A to FIG.3I illustrate a method of identifying alternative splicing events according to an embodiment of the present invention.

[0044] As shown in FIG.3A, the input to the processing was based on the Asian Immune Diversity Atlas (AIDA) cohort comprised 503 healthy individuals from East Asian 302, Southeast Asian 304, and South Asian populations 306. The scRNA-seq data generated from this cohort consisted of a total of around 1 million peripheral blood mononuclear cells (PBMCs) 308. Alternative splicing events were quantified in singlecell 310 and pseudobulk 312 resolutions for downstream analysis, including differential splicing 314, cis- and trans-QTL 316, dynamic sQTL 318, and colocalization 320.

[0045] RNA-seq data were aligned to the human reference genome GRCh38 primary assembly and GENCODE v32 using STARsolo v2.7.10a77 with options - soloCBmatchWLtype 1MM --soloUMIdedup 1MM_Directional_UMItools. Two-pass mode was used to enable novel splice junction discovery. The -waspOutputMode option was used to reduce allelic mapping bias. For each sample, the corresponding post-QC VCF file was used for WASP filtering. Deduplication of reads was performed based on cell barcodes (CB) and unique molecular identifier (UMI) tags from the bam files using MarkDuplicates from Picard. Only uniquely mapped reads in proper pairs that passed the WASP filter were retained for downstream analysis.

[0046] The methods next proceed in one of two ways. The first implementation uses samtools, LeafCutter, and QTLtools as described immediately below. The second implementation uses integrative single-cell splicing analysis and sQTL caller (ISSAC), below with reference to FIG.5A to FIG.5E.

[0047] Reads from the same cell type of each donor were extracted using custom scripts to make pseudobulk bam files, and reads from same individuals were pooled using samtools (v1.16.1) merge. Regtools were used to extract intron junctions and LeafCutter to quantify intron usage levels. The prepare_phenotype_table.py script from LeafCutter was used to generate phenotype files in sQTL mapping. Introns with zero read counts in more than 40% of the samples or with insufficient variation (standard deviation < 0.005) were removed.

[0048] cis-sQTL mapping was performed using QTLtools using intron excision ratios and a cis-window of 1 Mb on both sides of the junction. Eight PCs derived from splicing ratios, five genotype PCs, sex, and age information were used as covariates in the linear model. The number of phenotypes and genotype PCs are chosen to maximize sQTL discovery. Grouped permutations (-grp option) were used to jointly compute an empirical p-value over all intron clusters of a gene. QTLtools was run using the permutation mode (1000 permutations), and beta-approximated permutation p-values were adjusted for multiple test correlation using the qvalue package. The significance threshold was set at FDR < 0.05. Conditional sQTL analysis was done by forward stepwise regression followed by a backwards selection step. The gene-level significance threshold was set to be the maximum beta-adjusted p-value over all sGenes in a given cell type. A scan for cis-sQTLs using QTLtools were performed to correct for all previously discovered variants and all covariates. If the beta-adjusted P- value for the lead variant was insignificant at the gene-level threshold, the forward stage was complete, and the procedure moved on to the backward stage. If this p-value was significant, the lead variant was added to the list of discovered cis-QTLs as an independent signal and the forward step moved on to the next iteration. The backwards stage consisted of testing each variant separately, controlling for all other discovered variants.

[0049] FIG.3B shows a comparison of locations on a nuclide transcript of Read 1 and Read 2 5’ scRNA-seq. Read 1 of 5’ scRNA-seq was biased towards the transcription start site (TSS), and read 2 spread more evenly across the gene body.

[0050] FIG.3C shows base coverage rate per gene against read count as shows in FIG.3C, base coverage increased with read count. Left: fraction of base coverage across different read count bins (fraction of base coverage = covered bases / all bases). Right: a median of 85.3% of exonic bases are covered across all genes.

[0051] FIG.3D shows replication of LeafCutter intron discoveries in GENCODE, PacBio MAS-seq, and Snaptron. Upper panel: 59.3% of LeafCutter discoveries were annotated in GENCODE, and 85.9% replicated in PacBio long-read sequencing from four individuals. Lower panel: Close to 93% of detected splice junctions appeared in more than 1,000 samples, 98.8% in more than 100, and 99.5% in more than 10.

[0052] FIG.3E shows 21 distinct PBMC subtypes with sufficient cell counts. Cell types were colored by their hematopoietic lineage. Numbers below cell type labels indicated the sample size for sQTL calling.

[0053] FIG.3F shows the number of alternatively spliced (AS) genes detected per cell across 21 cell types at the single-cell level. Diamonds indicate the average number of detected genes (NODG) per cell. The dashed line indicated the number of AS genes detected by OneKIK. As shown in FIG.3F, a median of 1,146 (range of medians per cell type: 1,013 - 2,081 ) alternatively spliced (AS) genes per cell is observed.

[0054] FIG.3G shows NODG positively correlated with the number of AS genes. For the same increment in NODG, AIDA (5’ assay) detected a significantly higher number of AS genes than OneKIK (3’ assay). Each point represented a single cell. OneKIK data was analyzed with the identical pipeline, and the estimate was 4.3 times (1,146 vs.

[0055] 267; P < 2.2 x10-308) as much as OneK1K7. Because AIDA libraries have higher average sequencing depth than OneK1K (53,846 vs 33,656 reads / cell), it was considered whether this could explain their difference in AS gene detection. The number of detected genes (NODG) was used as a proxy for sequencing depth and observed that the number of AS and NODG per cell were highly correlated (Pearson’s r = 0.95). An average of 66% (ordinary least squares (OLS); 95% CI: 65.8% – 66.3%) of expressed genes had detectable AS events in AIDA, and only 12.1% (OLS; 95% CI: 12.0% - 12.2%) in OneK1K suggesting that difference in AS gene detection was not simply due to sequencing depth and should be attributed to exon painting and other factors.

[0056] FIG.3H shows number of detected AS genes per pseudobulk cell type. Intron junction usage was quantified on the pseudobulk level by grouping cells into pseudobulk.bam files per cell type per donor. Junction usage was quantified for each pseudobulk using LeafCutter, which grouped introns sharing a splice site into clusters and quantified junction usage as proportions between intron and cluster counts. LeafCutter detected a median of 7,721 (range of medians per cell type: 5,341 – 9,683) AS genes per pseudobulk.

[0057] FIG.3I shows the number of detected AS genes scaled with the number of cells in a pseudobulk. The number of AS genes is scaled as a sigmoid function with the number of cells per pseudobulk, saturating at ~11,500 genes (coefficient of determination R2 = 0.92). Because pseudobulk is less noisy and lower in sparsity than single-cell quantifications, we used LeafCutter for differential splicing and sQTL analysis and SpliZ quantifications only when single-cell analysis was required.

[0058] To identify alternative splicing mechanisms that mediated polygenic disease risk, GWAS summary statistics were compiled for 20 traits focused on Asian populations. Among these were autoimmune diseases: Graves’ disease (GD), rheumatoid arthritis (RA), and systemic lupus erythematosus (SLE); inflammatory diseases: atopic dermatitis (AD) and asthma; blood-related traits: white blood cell count (WBC) lymphocyte count (LYM), hemoglobin (Hb), platelet count (PLT), monocyte count (MON), basophil count (BAS), hematocrit (Ht), mean corpuscular volume (MCV), neutrophil count (NEU), mean corpuscular hemoglobin concentration (MCHC), mean corpuscular hemoglobin (MCH), red blood cell (RBC), eosinophil (EOS); and anthropometric traits: height and BMI.

[0059] FIG.4A to FIG.4H illustrate investigation of disease associated alternative splicing with single-cell resolution using an embodiment of the present invention.

[0060] To identify alternative splicing mechanisms that mediated polygenic disease risk, GWAS summary statistics were compiled for 20 traits focused on Asian populations. Among these were autoimmune diseases: Graves’ disease (GD), rheumatoid arthritis (RA), and systemic lupus erythematosus (SLE); inflammatory diseases: atopic dermatitis (AD) and asthma; blood-related traits: white blood cell count (WBC) lymphocyte count (LYM), hemoglobin (Hb), platelet count (PLT), monocyte count (MON), basophil count (BAS), hematocrit (Ht), mean corpuscular volume (MCV), neutrophil count (NEU), mean corpuscular hemoglobin concentration (MCHC), mean corpuscular hemoglobin (MCH), red blood cell (RBC), eosinophil (EOS); and anthropometric traits: height and BMI.

[0061] Colocalization was conducted between cis-sQTLs from 19 cell types and GWAS for 20 complex traits with COLOC65 and determined the proportion of GWAS loci that colocalized with sQTLs in any cell type.

[0062] FIG.4A shows cell-type-specific colocalizations between cis-sQTLs. The cis-sQTLs are from 19 cell types and 20 complex traits (AD, atopic dermatitis; SLE, systemic lupus erythematosus; Hb, hemoglobin; RA, rheumatoid arthritis; PLT, platelet count; WBC, white blood cell count; MON, monocyte count; BAS, basophil count; Ht, hematocrit; MCV, mean corpuscular volume; LYM, lymphocyte; NEU, neutrophil; GD, Graves’ disease; MCHC, mean corpuscular hemoglobin concentration; MCH, mean corpuscular hemoglobin; RBC, red blood cell; EOS, eosinophil; BMI, body-mass index).

[0063] FIG.4B shows Heritability enrichment (prop, h2 / prop. SNP) for 20 traits mediated by cis-sQTLs from 19 cell types. Autoimmune diseases and inflammatory diseases were highlighted in bold. The proportion varied across different GWAS traits, with immune- related diseases like atopic dermatitis (50%) and systemic lupus erythematosus (47.22%) having the highest proportion, while non-immunological traits like height (25.14%) and BMI (23.85%) having the lowest proportion of colocalization. The mean proportion for each trait category was highest for autoimmune diseases (41.30%), followed by inflammatory diseases and immune-related blood traits (33.33% and 34.19%, respectively), and lowest for non-immune-related blood traits and anthropometric traits (32.94% and 24.92%, respectively). We did not find a significant correlation between the proportion of colocalization events and GWAS sample size (Pearson’s r = -0.17). As an orthogonal validation, we applied stratified linkage disequilibrium score regression (S-LDSC) to estimate heritability enrichment from GWAS summary statistics by each cell type66. The heritability of SLE, AD, GD, and RA showed higher enrichment in PBMC sQTLs than other traits, whereas the heritability of anthropometric traits, especially BMI, was not enriched in PBMC sQTLs.

[0064] FIG.4C shows colocalization for 28 example sGenes across 19 cell types in the five disease traits. The shading of each circle indicated the associated disease trait. The inset showed the total number of colocalized loci across the five diseases. A total of 53 colocalized loci were identified among the five autoimmune and inflammatory diseases (COLOC H4 > 0.75). As a positive control, the results identified cell-type-specific colocalization between SLE GWAS and an sQTL that regulates the splicing of IRF5 exon one. IRF5 is a well-known risk gene for SLE67 and encodes a transcription factor involved in the production of type I interferons. The putative causal SNP rs2004640 disrupted the 5’ splice site of exon 1B, leading to nonsense-mediated decay (NMD) and downregulation of IRF5 expression.

[0065] FIG.4D shows Gene expression, eQTL, junction reads, sQTL, and H4 posterior probability (sQTL-GWAS colocalization) for TCHP across 19 cell types. High junction usage between exons four and five led to sQTLs and sQTL-GWAS colocalization. sQTL effects led to intron retention between exons four and five, which subsequently resulted in nonsense-medicated decay of TCHP transcript and manifested as eQTLs. No strong correlation was observed between TCHP expression and eQTLs.

[0066] Graves’ disease (GD), an autoimmune hyperthyroidism with a significantly higher incidence rate in Asian than in European population was studied. Out of the 16 GWAS loci for GD, 5 (31.25%) colocalized c / s-sQTLs in at least one cell type. Notably, out of the 5 GD GWAS loci that colocalized with c / s-sQTL, 4 (80%) did not colocalize with c / s-eQTLs (H4 < 0.75), indicating that c / s-sQTLs mediated additional GWAS risk loci independent of c / s-eQTLs. Colocalized loci showed varying degrees of cell-type specificity. Using H4 = 0.75 as a cutoff, two sGenes colocalized in a single cell type, one colocalized in three cell types, one colocalized in four cell types, and one colocalized in six cell types. Colocalization analysis captured TCHP as a cell-type-specific risk gene for GD, although rare cells might be underpowered to detect shared effects. It is hypothesized that cell-type-specific colocalization was driven by differential usage of the affected intron across cell types. Indeed, cell types with high usage of the intron junction between TCHP exons four and five harbored c / s-sQTLs, which led to cell-type-specific sQTL-GWAS colocalization. Furthermore, c / s-sQTL effects led to intron retention between exons four and five and NMD, which manifests as c / s-eQTLs. Trichoplein, the protein encoded by TCHP, has been shown to negatively regulate ciliogenesis and inhibit tumor growth. However, the role of trichoplein in Graves’ disease has not been elucidated. A recent study in the Japanese population identified a GWAS locus associated with Graves’ disease (P = 8.6x10-14).

[0067] FIG.4E shows cell-type-specific colocalization of TCHP between Graves’ disease GWAS and sQTLs in seven cell types. The GWAS lead variant rs74416240 had a significant genetic effect on TCHP exon four usages in monocytes, NK, CD4+, and CD8+T cells (FDR < 0.05) and strongly colocalized with TCHP c / s-sQTL.

[0068] FIG.4F shows Minor allele frequency of rs74416240 of five AIDA populations and five major populations in 1000 Genomes. This shows an East Asian bias of rs74416240 minor allele. Notably, the MAF of rs74416240 was high in East Asian (Japanese = 0.18; Chinese = 0.15, Korean = 0.13), modest in Southeast Asian (Malay = 0.04), and absent in South Asian populations within AIDA. We observed the same pattern of allele frequencies in 1000 Genomes populations, in which the minor allele of rs74416240 only appeared in EAS but no other populations.

[0069] FIG.4G shows Gene model of TCHP with three isoforms. rs74416240 was located in the 5’ splice site of the intron junction between exons four and five. The lead variant rs74416240 resided in the last nucleotide of TCHP exon four and was predicted to disrupt the 5’ splice site.

[0070] FIG.4H shows a minigene experiment to validate the effect of rs74416240 on TCHP exon four splicing in K562 cells. The Universal Minigene Vector (UMV) backbone alone corresponded to the band with the smallest molecular weight on the gel image. The test region, containing the 57-nt long exon four plus the 200-bp flanking sequences, was cloned into UMV. Two identical minigene constructs with one nucleotide difference at rs74416240 (ref = G; alt = A) were transfected into K562 cells. The reference allele (G) predominantly led to the normal isoform, and the alternate allele (A) led to intron retention. SpliceAl prediction suggested a high likelihood of donor loss at TCHP exon four 5’ splice site (Prob[donor loss] = 0.92). To validate this effect, we conducted a minigene experiment to test the effect of rs74416240 on TCHP exon four splicing in K562 cells. We used a universal minigene vector (UMV) constituting a CMV promotor and two constitutive MCAD exons. The test region contained 57-nt long exon four plus the 200-bp flanking sequences.

[0071] Two identical minigene constructs were transfected with one nucleotide difference at rs74416240 (ref = G; alt = A) and observed two distinct mRNA products. TCHP exon four inclusion and a low level of intron retention were observed for the reference allele (G) construct. The alternate allele (A) construct, where the 5’ splice site was disrupted, revealed a nearly complete intron retention isoform. These results suggest that rs74416240 is the causal variant for TCHP sQTL and possibly contributes to Graves’ disease.

[0072] FIG.5A to FIG.5E show a schematic of an integrative single-cell splicing analysis and QTL caller according to an embodiment of the present invention. The caller comprises three stages: The first stage integrates neighboring cells with similar gene expression profiles into metacells to mitigate sparsity in scRNA-seq data while allowing cell statedependent splicing analysis and QTL calling. The second stage constructs splicing phenotypes as the ratio of splice site usage, and the third stage uses binomial mixed effect models and score tests to map sQTLs. Because scRNA-seq data are often challenged by 5’ or 3’ sequencing bias, ISSAC quantifies splice site usage ratios to mitigate the influence of end-biased read distribution. Compared to isoform-based methods, ISSAC does not require any transcript reference and enables the discovery of novel cell-type-specific splicing events. Compared to intron-based methods such as LeafCutter, the caller does not rely on intron cluster generation, which is often challenged by uneven read coverage.

[0073] FIG.5A shows scRNA-seq preprocessing: After single cell clustering based on gene expression, single cells of one individual are divided into several metacells per cell type based on principal components embeddings obtained from scRNA-seq gene expression matrix.

[0074] FIG.5B shows Junction reads extraction: junctools and juncstats modules enable extracting UMI-based junction read counts from aligned RNA-seq files. To estimate splice site usage ratio, ISSAC first quantifies junction read counts for each sample using its junctools module. Different from bulk RNA-seq focused methods, the proposed method first collapses PCR-duplicate reads using cell barcode and UMI to mitigate high PCR amplification biases in single-cell data.

[0075] The splice ratio is estimated at the splice site level. First junction read counts are quantified for each sample using the junctools module. PCR-duplicate reads are collapsed using cell barcode and UMI to mitigate high PCR amplification biases in single-cell data. Reads mapped to the same junction with the same cell barcode and UMI will be retained only once.

[0076] For each splice site, the usage ratio is quantified as the proportion of junction reads that included this splice site, relative to all junction reads that could include it. To determine which reads could include the splice site, its splice partners

[0077] - nearby splice sites that joins with the target site are identified to define introns supported by junction reads - and then all junction reads involving these splice partners are counted.

[0078] Site Usage Ratio is defined as the junction reads used by the site divided by all the junction reads related to the site which equal to the sum of the junction reads used by the site and junction reads competing the usage of the site. For a splice site without any competing junction reads - such as intron retention or unprocessed nascent transcript - the splice site usage ratio is calculated as the number of junction reads supporting the splice site divided by the total number of reads associated with the site, which is the sum of junction reads supporting the site and reads crossing the site without splicing. Then the splice read count of the site and total read count will be used in to fit the null binomial mixed model.

[0079] FIG.5C shows Splice site partner definition: junction reads supporting the usage of the splice site c and junction reads competing the usage of the splice site c.

[0080] Splice site usage ratio is calculated using the following equation.

[0081] SUCc=

[0082] (fi

[0083]

[0084] The numerator is defined as the sum of junction reads utilizing site c (I1 + I2) and the denominator is defined as the sum of supporting junction reads (h + I2) and junction reads competing the usage of site c (E1).

[0085] Reads mapped to the same junction with the same cell barcode and UMI will be retained only once. For each splice site, the usage ratio is quantified as the proportion of junction reads that spanned this splice site, relative to all junction reads that could span the splice site. To determine which reads could span the splice site, splicing partners are identified. The splicing partners are nearby splice sites that join with the target site to form introns supported by junction reads and then counts junction reads originating from these splice partners.

[0086] FIG.5D shows Binomial mixed model construction: a null binomial mixed model is constructed by estimating parameters for fixed effect such as sex, age, ancestry PCs et al and variance components for random effect represented by genetic relationship matrix; associations between genotype and splice site usage ratio are estimated through single-variant score tests after null model construction. To map sQTL for each splice site, first fits a null generalized mixed-effect model is fitted to estimate coefficients of fixed effects and variance components and then performs score tests to test for the association between each genetic variant and splice site usage ratio The details of fitting null binomial mixed model and estimating parameters for both fixed effects and variance components are introduced.

[0087] The model can be written as

[0088] logit(πi) = χi+ βi+ bi

[0089] Here Xi denotes fixed covariate such as sex, age and ancestry PCs for donors and bias the random effect is used to account for intraindividual sample-to-sample variability. Penalized quasi-likelihood (PQL) and REML methods are used to iteratively estimate the model parameters. Below is the parameter estimation procedure:

[0090] Under the null hypothesis that genotype has no impact on the model, η = χβ + b is the linear predictor in the model fitting process with χβ as fixed effects and b as random effects.

[0091] Xi denotes one splice site’s splice read count in sample i;

[0092] rii denotes one splice site’s total read count in sample i;

[0093] b denotes random effects and b corresponds to ZV(0, TG / ? M);

[0094] Link function is 77 = log( ^j for binomial regression.

[0095] Here yt~ Binomial (n^ni ) and variance for ytVar yP) — nin^l — TTP).

[0096] The quasi-likelihood of the observations conditional on the random effects can be denoted as ql y\b~) and the log-likelihood of random effects distribution could be denoted as Z(Z?)Thus the joint log quasi likelihood is gZ(y|h) + Z(b) and marginal quasi likelihood is: | [ql(y|b) + l(b)]db — ql(y) =

[0097] l?b"

[0098] log

[0099]

[0100] And here ql(yL\b) = quasi-likelihood for the ith individual when

[0101]

[0102] the random effect b is given and cq is a constant which could be omitted. Here we need to iteratively maximize the ql(y\b~) and Z(h).

[0103] Then working response variable and working weights are obtained as following at each iteration:

[0104]

[0105] Var(πi) = diag(πik(1 − πik))

[0106] Here, the working response helps linearize the nonlinear relationship in the binomial model and working weight ensures observations with higher variance contribute less to parameter updates. The variance for the working response variable could be written as

[0107]

[0108] The initial number of τ was set to 0.

[0109] To estimate β and u, the REML pseudo-log-likelihood can be written as

[0110]

[0111] Note: p = rank(X)

[0112] The detailed estimation procedures are as follows:

[0113] 1) Fit null model: Here, iteratively reweighted least squares (IRLS) are used to obtain initial estimates of fixed effect β0and set initial random effect u0= 0 (τ = 0); The IRLS algorithm iterates over the following steps until convergence: (i) initialize the fixed effect coefficients p as 0 (ii) update the linear predictor η = Xβ; update the mean response with inverse logit function πik= 1 / (1+exp(-ηik)); update working response ηik=

[0114]

[0115] χβ + ... update working weights wik= πik(1 - πik); (iii) update β using

[0116]

[0117] weighted least squares equation pi+1= XTWX)~1XTWr]l (iv) check for convergence until |?i+1- Pt | < threshold (default: 0.001).

[0118] 2) After obtaining initial estimates for fixed effects, we started estimating for variance components.

[0119] At iteration k, update working response, working weights and variance, πik= 1 / (1+exp(-ηik))

[0120] Working response: ηi* = χβ + u + (yik- πiknik) / (nikπik(1-πik)) working weights: wik= πik(1 - πik) Variance of working response: Var(ηi*) = wik-1+ τGRM

[0121] 3) At each iteration time, estimate variance components T given p, u >'n - px

[0122] pMJ y*. = ~ I— y~~) log(. logt log |X* {VarC f “iX|

[0123]

[0124]

[0125]

[0126] Here in order to maximize above REML pseudo-log-likelihood, we used “opt. optimize” function based on nlopt library in C++ to optimize the function to obtain T which could maximize the likelihood (https: / / nlopt.readthedocs.io / ).

[0127] 4) Then we obtained the converged τ to update β, u, Var(η*);

[0128]

[0129] Iterate step (2) ~ (4) until < tol. Default tol = 10'6

[0130] Step3: After obtaining ηnullfrom above process, πnull= 1 / (1+exp(-ηnull)) and residualsnull= y - n × πnullwere obtained. Then test statistics is computed as following:

[0131]

[0132] Here, G = Goriginal- X(XTWX)-1XTGoriginalis the covariate-adjusted genotype. To estimate the normalized parameter r to control Type I error rate, a permutation method is performed by randomly simulating a certain number of genotypes with MAF > 0.05 (default: n = 1000) and obtained corresponding test statistics. The test statistics are assumed following normal distribution under null hypothesis with mean 0. Then the variance for the test statistics could be obtained and set as the parameter to normalize the distribution of test statistics to a normal distribution N(0,1). Then, corresponding P values could be obtained from variance-adjusted test statistics and effect size was defined as the normalized regression slope between covariate-adjusted genotype and adjusted phenotype

[0133]

[0134] To improve statistical power while maintaining low sparsity intrinsic to single-cell data, the phenotype is prepared for downstream cis-sQTL based on the concept of metacells. For each donor, initial estimates for each cell cluster were obtained based on PC embedding obtained from snRNA-seq gene expression matrix through k nearest neighbors’ graph construction and Louvain clustering (with m as the default parameters of the size of each cluster). Then clusters are iteratively decomposed with less than m cells to adjacent larger clusters to ensure each cluster having at least m cells. Finally, each cluster with at least m cells was treated as one metacell for downstream sQTL mapping tasks.

[0135] FIG.5E shows Downstream analysis: cell-type specific sQTLs detection, cell-state dependent sQTLs identification, context-biased sQTLs discovery and colocalization analysis are enabled, used a binomial random variable to model the splice site usage ratios, and fixed effects to capture the genetic regulatory effects and potential confounders such as sex, age, ancestry, and PEER factors or splicing PCs. Because each individual has multiple metacells per cell type, a genetic relationship matrix (GRM) is used to model repeated measurements within the same individual. As an additional benefit, the GRM explicitly models genetic relatedness across donors and obviates the need to remove genetically related donor from population-scale cohorts. Penalized quasi-likelihood (PQL) and REML are used to iteratively estimate fixed effects and the variance component, respectively. To mitigate the computational burden associated with matrix inversion and determinant calculation, the preconditioned conjugate gradient (PCG) method is employed to approximate the product of the inverse matrix and a vector, thereby circumventing the need for explicit matrix inversion during variance component estimation under the REML framework.

[0136] Next, perform score tests are performed to obtain corresponding P-values for each genetic variant. Effect size was defined as the regression slope between genotype and adjusted phenotype:

[0137]

[0138] FIG.6A to 6F show a comparison of an integrative single-cell splicing analysis pipeline according an embodiment of the present invention and a LeafCutter pipeline. The performance of the integrative single-cell splicing analysis pipeline in null and causal settings was evaluated to assess calibration and power. The integrative single-cell splicing analysis pipeline, comparing ISSAC to LeafCutter, the only method applied to map single-cell pseudobulk-level sQTL so far. Because the integrative single-cell splicing analysis pipeline and LeafCutter use different phenotype quantification, single-cell isoform expression was simulated with scIsoSim to ensure a fair comparison and minimize bias towards either phenotypic setup. First, both methods were assessed with null simulations, where no genetic variant is associated with splice site usage. To ensure robustness, we simulated lower (N=20) and higher (N=100) cell counts per individual.

[0139] FIG.6A shows null simulations of the integrative single-cell splicing analysis pipeline and the LeafCutter pipeline. The dark points 602 indicate the an integrative single-cell splicing analysis pipeline and the light points 604 indicate the LeafCutter pipeline. The P values shown in FIG.6A where obtained when the pipelines were applied to sitebased phenotype and intron-based phenotype in the simulated scRNAseq datasets; 1000 random SNPs with MAF > 0.05 were generated to measure the calibration of both pipelines.

[0140] FIG.6B shows a comparison of the power of the integrative single-cell splicing analysis pipeline and the LeafCutter pipeline when applied to simulated scRNA-seq datasets with various effect size. The dark points and line 606 indicate the an integrative singlecell splicing analysis pipeline and the light points and line 608 indicate the LeafCutter pipeline. The integrative single-cell splicing analysis pipeline and the LeafCutter pipeline were evaluated in causal simulations where a subset of genes has genetically regulated isoform expressions. The same individuals and cells from the null simulation were used and 334 genes were randomly selected as causal genes, whose isoform expressions are associated with genotypes. It was observed that the integrative singlecell splicing analysis pipeline obtained 16.0-40.3% higher power than LeafCutter across a wide range of effect sizes (β=1.1 to 1.9) at a false discovery rate (FDR) < 0.05. FIG.6C shows a comparison of the power of the integrative single-cell splicing analysis pipeline and the LeafCutter pipeline when applied to simulated scRNA-seq datasets with different sample size. The dark points and line 610 indicate the an integrative single-cell splicing analysis pipeline and the light points and line 612 indicate the LeafCutter pipeline. Sample sizes of N = 50, 100, 200, 400 were used. For common cell type, each individual has 100 cells and for rare cell type, each individual has 20 cells. As shown in FIG.6C, the integrative single-cell splicing analysis pipeline obtained 15.0-72.6% higher power than LeafCutter across a range of sample sizes (N=50 to 400). The power improvement is more pronounced in low-abundance cell types (N=20) than highly abundant cell types (N=100).

[0141] FIG.6D shows overall performance in terms of sGenes & ratios of significant phenotype when applying the integrative single-cell splicing analysis pipeline and the LeafCutter pipeline to three cell types with varied sample size in the 10X 5’ scRNA-seq AIDA phase I freeze II cohort and the same genes’ sites & introns were considered here to ensure comparison fairness.

[0142] FIG.6E shows Type I error rate of the integrative single-cell splicing analysis pipeline and the LeafCutter pipeline when applied to gdT GZMBhi, Naive B and CD16+ NK coming from 10X 5’ scRNA-seq AIDA cohort; 50 splice sites and 50 splice introns were randomly selected from each cell type and 250 times’ permutations were applied to the phenotype ID to disrupt the potential associations between variants +-1 MB around the splice site(intron) and splice sites(intron). The integrative single-cell splicing analysis pipeline’s power improvement is larger for rarer cell types, agreeing with the simulation results above. To assess the type I error rates of ISSAC and LeafCutter, 50 sGenes were randomly selected from each of the three cell types and the association between genetic variants and the phenotype was disrupted by permuted each phenotype 250 times. Using an FDR cutoff of 0.05, our results suggested that both ISSAC and LeafCutter has well-controlled type I error rates as false-positive associations remained below 5%.

[0143] FIG.6F shows time consumption of the integrative single-cell splicing analysis pipeline and R package Ime4 glmer function when applied to site-based phenotype generated by three cell types with varied sample size (gdT GZMBhi, Naive B & CD16+ NK) in the 10X 5’ scRNA-seq AIDA cohort. To assess performance on a real-world dataset, the integrative single-cell splicing analysis pipeline was applied to the Asian Immune Diversity Atlas datasets (phase I freeze II), for which single-cell sQTL was mapped by LeafCutter on the pseudobulk level. Three cell types with low (gdT GZMBhi; N = 11860; 39.14 cells / donor), medium (Naive B; N = 36778; 65.91 cells / donor) and high (CD16+ NK; N = 151,655; 259.68 cells / donor) cell counts were selected, and obtained 1093, 3287, and 15310 metacells for the three cell types, respectively. To ensure fairness, we used the same set of splice junctions for the integrative single-cell splicing analysis pipeline and LeafCutter phenotypic construction, as well as the same set of tested genes for gdT GZMBhi (N = 2103), Naive B (N = 1957), and CD16+ NK (N = 2081). ISSAC identified 2.4-fold (267 vs. 111), 1.57-fold (371 vs. 236), and 1.55-fold (771 vs. 496) as many sGenes as LeafCutter.

[0144] FIG.7A to FIG.7H show an application of an integrative single-cell splicing analysis pipeline according to an embodiment of the present invention to identify cell typespecific genetic regulation of splicing in the dorsolateral prefrontal cortex.

[0145] The human brain exhibits exceptionally high transcriptomic diversity, characterized by extensive alternative splicing and a rich repertoire of isoforms. Although prior studies have explored the genetic regulation of splicing at the bulk tissue level, and cell typespecific genetic regulation of gene expression, cell type specific genetic regulation of splicing in the brain remains unexplored. The integrative single-cell splicing analysis pipeline was applied to identify cell type-specific genetic regulation of splicing in the dorsolateral prefrontal cortex (DLPFC) using snRNA-seq data from two aging cohorts.

[0146] FIG.7A shows schematics of cell type specific sQTL analysis in DLPFC(BA9) regions.

[0147] 3 million cells from De Jager cohort and Tsai cohort were integrated and divided into 7 major cell types & 67 subcell types for downstream sQTL analysis. As shown in FIG.7A, after preprocessing and quality control, 3,177,748 cells of 722 specimens from 530 donors (192 donors shared across two cohorts) with both snRNAseq and wholegenome sequencing (WGS) data were retained for analysis. Following previous singlenucleus studies of the DLPFC, the combined dataset was classified into seven major cell types, including excitatory neurons, inhibitory neurons, astrocytes, oligodendrocytes, oligodendrocyte progenitor cells, and endothelial cells. FIG.7B shows numbers of sGenes and splice sites (Bonferroni adjusted P value < 0.05) detected in the seven major cell types (top right); Numbers of shared significant splice sites between each major cell type (bottom left), each line connecting black dot representing one shared set. Among all the 33,512 sSites (42,998 cis-sQTLs) that we detected across all cell types, 27,391 (81.7%) were specific to one of the seven major cell types, and 6,121 (18.3%) were shared between at least two cell types

[0148] FIG.7C shows the relationship between number of sGenes and number of cells in the 7 major cell types. The shaded area on either side of the linear regression line represents the 95% Cl.

[0149] FIG.7D shows the Relationship between number of sGenes and number of cells in the 67 subcell types. The shaded area on either side of the linear regression line represents the 95% Cl.

[0150] FIG.7E shows the numbers of sGenes and splice sites (Bonferroni adjusted P value < 0.05) detected in the 67 subcell types (bottom left); The proportion of sGenes shared by both major cell types and subcell types are indicated in the top right plot.

[0151] FIG.7F shows fold enrichment of sSNPs in functional annotation category; the dashed line indicates the threshold of enrichment.

[0152] FIG.7G shows sharing between sGenes and eGenes detected from seven major cell types. The distribution plot indicates the PP3 and PP4 of colocalized results between shared sGenes and eGenes. The proportion of eQTL and sQTL effects that share the same underlying genetic regulatory mechanism was estimated. By comparing the sGenes to eGenes from an existing study using the same cohort, it was found that an average of 70.3% (47.7% - 94.9%) sGenes were not eGenes in their corresponding cell type. Among the overlapping sGene-eGenes, 21.7% have H4 >0.75, indicating shared causal variants, while the majority of 69.2% have H3 > 0.75, indicating distinct regulatory variants. FIG.7H shows an example of celltype specific sQTLs; the lead SNP rs72786713 (C> G) modulated the increase usage of the splice site chr16:+:56683886 only significantly in astrocytes but not in another six major cell types. Excitatory neurons drove the largest cell type-specific sQTL discovery due to its high abundances in the DLPFC, and numerous cell type-specific sQTLs were discovered within other cell types. One example of astrocytes-specific sQTLs is MT1X. The mutation of lead SNP rs72786713 (C> G) led to increased usage of chr16:+:56683886 (FDR = 1.014x10-16, p = 0.09) only in astrocytes but not in other cell types.

[0153] Overall, the results indicated that the majority of cell type-specific sQTLs and eQTLs were regulated by distinct genetic mechanisms.

[0154] FIG.8A to FIG.8I show an application of an integrative single-cell splicing analysis pipeline according to an embodiment of the present invention to study Intron retention sQTLs and context-biased sQTLs.

[0155] Among the sQTL-eQTL colocalization events, we further investigated the shared genetic regulatory mechanism. Previous studies observed that alternative splicing, such as poison exon or intron retention, could lead to premature stop codon and nonsense mediated decay.

[0156] FIG.8A shows Numbers of sGenes discovered based on single intron clusters colocalizing with corresponding eGenes in major cell types (H4 > 0.75). A total of 298 colocalizing sQTLs observed in 230 sGenes regulated single-intron clusters, indicating potential intron retention events.

[0157] FIG.8B shows a violin plot showing the mutation of rs6713 (C> T) modulated the decrease usage of splice site chr2:+:71078317 within NAGK. In excitatory neurons, the 5’ splice site of NAGK exon 10 (chr2:+:71078317) was associated with the lead SNP rs7613 (G> A) (P =4.24x10-71, / 3 = -0.174).

[0158] FIG.8C shows Colocalization event between NAGK eQTLs and sQTLs near the splice site chr2:+:71078317; rs6713 was the lead SNP for both eQTLs and sQTLs located 9 bp downstream of the splice site, and NAGK sQTLs-eQTLs colocalized with an H4 of 0.99.

[0159] FIG.8D shows an Intron retention model of NAGK indicates when lead SNP rs6713 mutates from C to T, Intron9 between exon9 and exonlO could be retained and that led to nonsense mediated decay, further contributing to downregulation of NAGK and protein aggregates such as mHtt or a-syn not cleared. The intron retention event between exons 9 and 10 within NAGK likely led to nonsense-mediated decay and downregulation of NAGK expression. NAGK was known for interacting with the dynein light chain roadblock type 1 to clear protein aggregates. The downregulation of NAGK may impair its ability to clear protein aggregates and contribute to neurodegeneration.

[0160] FIG.8E shows Manhattan plots summarizing AD-biased sGenes in the major cell types. The FDR in the y axis - log10(FDR) represents the interaction term FDR between splice phenotype and the interaction term GenotypexAD. AD pathology may alter the cellular and tissue microenvironment, leading to AD-biased genetic regulation. Conversely, certain inherited splice-regulatory variants may alter the risk of AD. For each significant lead c / s-sQTL, we performed AD-biased sQTL analysis with a binomial mixed model to test for genotype -by- AD interaction. We identified a total of 372 AD-biased sQTLs in 183 sGenes across seven major cell types (Interaction FDR < 0.01)

[0161] FIG.8F shows a violin plot showing the associations between the lead SNP rs72786713 and splice site chr16:+:56683386 of MT1 X in astrocytes is ADbiased. One example of AD-biased splice sGenes is MT1X in astrocytes (Fig. 4f). The lead SNP rs72786713 (C> G) modulated the splice sites (chr16:+:56683386) only in AD patients but not in healthy controls (control: P = 0.362, / 3 = 0.001; AD: P = 9.53 X 10-41, / 3 = 0.142). MT1X expression was previously reported to increase along the temporal progression of AD in astrocytes, which may explain AD-biased sQTL effects on MT1X. Notably, MT / Xwas also a cell type-specific sGene in astrocytes but not in other cell types. FIG.8G shows Numbers of sex-biased splice sites in the major cell types (GenotypexSex interaction FDR < 0.01). Similar to AD-biased sQTLs, genome-wide sex-biased sQTLs detection was performed with a binomial mixed model to test for genotype-by-sex (GxS) interaction. 463 sex-biased sQTLs were identified in 207 sGenes across seven major cell types.

[0162] FIG.8H shows a violin plot showing the association between the lead SNP rs12551220 and splice site chr9:+:136980855 of PTGDS in oligodendrocyte is female specific. One example is that, in oligodendrocyte, the lead SNP rs12551220 (G> A) modulated the splice site chr9:+: 136980855 within PTGDS only in females but not in males (Fig. 4h; male: P = 0.992, [3 = 0.038; female: P = 1.60x10-39, p = -0.248). PTGDS functions by converting prostaglandin H2 to prostaglandin D2 and plays a role in preventing the aggregation of amyloid-beta peptides, thereby conferring protection against Alzheimer’s disease.

[0163] FIG.8I shows a violin plot showing the association between the lead SNP rs17426032 and splice site chr5:+:94649597 of SLF1 in inhibitory neurons is male-specific. SLF1 is an example of male-biased sQTLs. Within inhibitory neurons, the lead SNP rs17426032 (T> C) modulated the usage of splice site chr5:+:94649597 only in males but not in females (male: P = 3.14×10-12, β = -0.120; female: P = 0.178, β = -0.019). SLF1 has been reported to participate in the DNA damage response pathway58, and has been implicated in two neurological disorders, form agnosia and Alcardi syndrome.

[0164] FIG.9A to FIG.9F show an application of an integrative single-cell splicing analysis pipeline according to an embodiment of the present invention to study cell-state dependent sQTLs identified in inhibitory neurons. Regulatory effects of non-coding variants are often modulated by continuous cell states. Previous studies have shown that up to one-third of cis-eQTLs were influenced by multimodally defined cell states. However, there has been minimal direct evidence of cell state-dependent genetic regulation of alternative splicing. Pseudobulk splicing quantifications can only detect mean cell-type effects and would mask cell state-dependent sQTLs. In contrast, metacells from the integrative single-cell splicing analysis pipeline allow the capture of continuous cell states using low-dimensional embedding from principal components analysis (PCA). Here, inhibitory neurons are focused on due to their relative abundance and elevated cellular diversity compared to other cell types. Inhibitory neurons regulate brain activity by modulating the activity of excitatory neurons and play a critical role in synaptic plasticity.

[0165] FIG.9A shows Gene Ontology BP enrichment results of PC1 marker genes (P< 0.05). PC1 is enriched with dendrite growth genes in Gene Ontology (GO) enrichment analysis.

[0166] FIG.9B shows Gene Ontology BP enrichment results of PC3 marker genes (P< 0.05). PC3 is enriched with synaptic organization genes in GO enrichment analysis.

[0167] FIG.9C shows examples of cell-state dependent sQTLs; the heatmap means effect size of the cell-state dependent sQTLs within PC-low (low 1 / 3 of PC), PC-mid (mid 1 / 3 of PC), PC-high (high 1 / 3 of PC) groups, respectively with the color bar corresponding to the cell-state represented by PC1~8. To map cell state-dependent sQTLs, genotype-by-cell-state (GxS) interaction effect was added into the binomial mixed-effect model. Here, 369 c / s-sQTLs were identified to be cell state-dependent along the eight PCs in the major cell types (Interaction FDR < 0.05).

[0168] FIG.9D shows a heatmap coloured by Pearson’s correlations between PCX(1 ~8) and selected marker genes’ normalized expression. Each metacell is scored along the top eight PCs and observed that each PC correlated with well-defined inhibitory neuron functions. ADARB2 and SOX6 expression contributed most to PC1 (Pearson’s r = 0.85 and -0.81, respectively). ADARB2 is key marker gene in central ganglionic eminence (CGE)-derived inhibitory neurons, and its mutations has been linked several neurological disorders. SOX6 is a gene that ensures proper inhibitory neurons formation and function. ERBB4 and RALYL expression were highly correlated with PC3 (Pearson’s correlation = -0.71 and 0.65, respectively). ERBB4 has been reported to encode a receptor tyrosine kinase and promote inhibitory synapse formation, and RALYL encodes an RNA binding protein and its decreased expression was correlated with the gradual progression of Alzheimer’s disease. FIG.9E shows a Uniform Manifold Approximation Projection (UMAP) plot of PC1 and one example of PC1-dependent sQTL PREX2 and Violin plots show the varied sQTL effect for cells in PC1-low (right), PC1-mid (left top) and PC1-high groups (left bottom). PREX2 sQTL was one example showing interaction with PC1 (interaction P = 1.46 x 10-7). The genetic effect between lead SNP rs7817030 (A> G) and the splice site (chr8:+:68053096) in PREX2 showed gradual increase along PC1 (PC1 -low P = 0.591, P = -0.002; PC1-mid P = 8.88 x 10-5, (3 = 0.074; PC1-high P = 2.27 x 10-10, β = -0.153).

[0169] FIG.9F shows a UMAP plot of PC3 and one example of PC3-dependent sQTL LSAMP; Violin plots show the varied sQTL effect for cells in PC3-low (right top), PC3-mid (left) and PC3-high groups (right bottom). PREX2 has been reported to be involved in dendrite morphology and synaptic plasticity in Purkinje cells, congruent with PC1's representation of the dynamic changes of dendrite development. Another example is LSAMP; the effect size of its lead SNP rs9944185 (C> A) gradually decreased as PC3 increased (Fig. 5f; PC3-low P = 1.04 x 10-35, β = -0.217; PC3-mid P = 2.72 x 10-18, β =-0.198; PC3-high P = 3.14 x 10-6, β = -0.112). LSAMP was known to participate in the generation of a neuronal surface glycoprotein which acts as an adhesion molecule during axon guidance and synaptic formation. Cell state-dependent splicing observed in LSAMP aligns with PC3’s role in the dynamic remodeling of synaptic organization.

[0170] These results highlight the critical role that cell states play in the genetic regulation of splicing.

[0171] FIG.10A to FIG.10H show an application of an integrative single-cell splicing analysis pipeline according to an embodiment of the present invention to study mechanisms contributing to neuronal diseases. To prioritize functional genes for complex neurological disorders, colocalization was performed to assess whether genetic association with sGene and trait is driven by the same underlying causal variants. Here, COLOC analysis was performed between the seven major cell types and six well-powered neurological disorder GWAS including Alzheimer’s disease (AD), Amyotrophic lateral sclerosis (ALS), Parkinson’s disease (PD), Schizophrenia (SCZ), Lewy body dementia (LBD), and Neuroticism. FIG.10A shows colocalization events between cell-type specific sGenes with risk SNPs of AD, ALS, PD, SCZ, LBD and Neuroticism GWAS. The piechart reports H4 over 0.75 of colocalization events detected by COLOC method. A total of 349 colocalization events across 142 sGenes with the six neurological diseases (H4 > 0.75) were identified.

[0172] FIG.10B shows Heritability enrichment (proportion h2 / proportion variant) for five neuron-related traits and one non-neuronal trait BMI mediated by cis-sQTLs from six major cell types including Exc, Inh, Oli, OPC, Ast and Mic. The dashed line represents the threshold for enrichment tendency within the trait GWAS for the cell type-specific sQTLs. To understand implications of cell types in each disorder, stratified LD score regression (S-LDSC) was applied to estimate heritability enrichment of cell typespecific sGenes within GWAS summary statistics. Numerous cell type-specific enrichments across six GWAS (Enrichment Score > 1) were observed. For instance, microglia sGenes were highly enriched in AD and ALS GWAS but not for the other five traits. It has been observed that most AD risk genes were highly or exclusively in microglia, and activated microglia has been reported as a neuropathological hallmark of ALS80. Furthermore, oligodendrocyte sGenes were highly enriched in Parkinson’s disease. Oligodendrocyte-specific gene expression was reported to be distinctly associated with PD risk, but neuroinflammation-related cells, such as microglia, play a less causal role in PD81. As a quality control, BMI82 was used as negative controls and observed a lack of overall enrichment.

[0173] FIG.10C shows colocalization plots between SCZ GWAS and Exc, Inh, Oli sQTL around GNL3 loci. The colocalization between the GNL3 sQTL and SCZ GWAS was cell type-specific (inhibitory neurons: H4 = 0.98; excitatory neurons: H4 = 0.77, oligodendrocyte: H4 = 0.48) and was therefore not identified by previous bulk-level sQTL studies.

[0174] FIG.10D shows a Gene model demonstrating the mechanism how sSNP rs1108842 regulates the alternative 3’ splice site of exonl within GNL3. The putative causal SNP rs1108842 (A> C) mediated the splice site usage of chr3:+:52686105 (excitatory neurons: P = 2.09 x 10-22; β = 0.084; inhibitory neurons: P = 0.0001; β = 0.061; oligodendrocyte: P= 0.64; β = -0.004) and chr3:+:52686202 (excitatory neurons: P = 2.44 x 10-16, β = -0.077; inhibitory neurons: P = 2.29 x 10-7; β = -0.056; oligodendrocyte: P = 0.0431; β = -0.055) in a cell type-specific manner. The shift of GNL3 exonl 3’site from chr3:+:52686202 to chr3:+:52686105 increases the risk of Schizophrenia.

[0175] FIG.10E shows colocalization plots between Neuroticism GWAS and Exc, Inh, Oli, OPC, Ast sQTL around TRPT1 loci. To further functionally validate the mechanism of putative regulatory variants, a colocalization event was identified between neuroticism and TRPT1 sQTL, appearing in ectodermal cells including excitatory neurons, inhibitory neurons, OPC, oligodendrocyte and astrocytes.

[0176] FIG.10F shows a sashimi plot demonstrates the exon skipping event mediated by the lead SNP rs11549690 mutating from G to A. The highlighted numbers represent the UMI-based junction reads of the intron connecting exon6 and exon8. TRPT1 was predicted to enable tRNA 2’ - phosphotransferase activity and catalyze the last step of tRNA splicing. The lead SNP rs 11549690 (G> A) modulated the usage of the splice site chr11:-:64224099 of exon 7 (excitatory neurons: P = 7.37 x 10-89, β = -0.180; inhibitory neurons: P = 3.89 x 10-15, β = -0.101; oligodendrocyte: P = 2.31 x 10-20, β = -0.139; OPC: P = 6.58 x 10-18, β = -0.167; astrocytes: P = 5.98 x 10-22, β = -0.177). The lead SNP rs11549690 is located 10 bp upstream of the splice site chr11:-:64224099. Switching from the reference (G) to alternative allele (A) increased exon-7 skipping read counts from 50 to 322 (6.4-fold).

[0177] To narrow down on the most likely causal variant, SpliceAl was used to predict the delta score induced by each variant. The lead variant was predicted to have the highest delta-donor score (ΔDonor = 0.21) on the splice site chr11:-:64224099, and we proceeded with functional validation to test the effect of rs11549690 on TRPT1 exon 7 (114-nt) skipping in the SH-SY5Y neuroblastoma cell line. A minigene construct whose test region consisted of the full sequences of TRPT1 exons 6, 7 and 8 and intervening introns was designed. Two minigene constructs that are identical in sequence, except for a one nucleotide difference at rs11549690 (G or A), into SH-SY5Y cells were transfected. Both minigene constructs showed the presence of bands corresponding to exons 6 and 8 with either exon 7 inclusion (around 400-bp) or exon 7 exclusion (around 300-bp), but with the reference (G) allele showing a higher percentage of inclusion, and the alternative allele (A) showing a shift towards a higher percentage of exclusion. These functional validation experiments further support the causal effect of this SNP to the exon skipping events within TRPT1.

[0178] FIG.10G shows a gene model demonstrating the mechanism how sSNP rs11549690 regulates the exon skipping events between exon6 and exon8.

[0179] FIG.10H shows a minigene experiment to validate the effect of rs11549690 on TRPT1 exon 7 splicing in SH-SY5Y neuroblastoma cell line. The empty pcDNA3.1 (+) backbone alone corresponded to the band with the smallest molecular weight on the gel image. The test region containing TRPT1 exon 6, exon 7 and exon 8 was cloned into the pcDNA3.1(+) plasmid. Two identical minigene constructs with one nucleotide difference at rs11549690 (reference = G; alternative = A) were transfected into SH-SY5Y neuroblastoma cell line. The reference allele (G) predominantly led to the normal isoform with exon7 inclusion; the alternative allele (A) led to exon 7 exclusion.

[0180] These results highlight rs11549690 is a causal variant for TRPT1 sQTL and contributes to Neuroticism risk via this exon skipping event. In conclusion, these findings provided potential target genes whose alternative splicing mediates complex neurological diseases at the single-cell level.

[0181] As described above, the present disclosure provides an integrative single-cell splicing analysis pipeline, a method to map cell type-specific and cell state dependent sQTLs using sc / snRNA-seq data. Different from contemporary pseudobulk analysis, the method provided in the present disclosure directly models splice site usage ratios on the metacell-level using a binomial mixed-effect model.

[0182] We have shown via extensive simulation and real-world data that the integrative single-cell splicing analysis pipeline was well-calibrated and detected 1.5- to 2.4-fold as much sQTLs in diverse cell types compared to pseudobulk methods, with a higher increase in power for low-abundance or rare cell types. Ablation studies demonstrated that the integrative single-cell splicing analysis pipeline derived its power boost from a combination of splice site-based phenotypic construction, binomial distributional assumption, and metacell quantification with mixed-effect models. To accommodate for ever-increasing sample sizes, the integrative single-cell splicing analysis pipeline was implemented in C++ with several modern techniques to improve runtime speed and achieved 11.2- to 115.1 -fold faster speed than glmer function in the Ime4 R package.

[0183] The integrative single-cell splicing analysis pipeline has been applied to address the critical gap in identifying cell type-specific and cell state-dependent sQTLs in the human DLPFC, a brain region characterized by substantial isoform diversity and functional heterogeneity. More than 80% of sSites are unique to one of the seven major cell types, and thousands of sSites were identified only in subcell types and not in major cell types, highlighting the heterogeneity in splicing regulation. To address a key limitation of pseudobulk splicing quantification, the integrative single-cell splicing analysis pipeline may be used to uncover hundreds of cell state-dependent sQTLs in inhibitory neurons, revealing genetic effects that are modulated by dendritic development and trans-synaptic signaling. We prioritized 142 risk genes for neurological disorders, including MTCH2, GNL3, and TRPT1. Many colocalization events are specific to one single cell type, highlighting the functional heterogeneity of disease pathology. We used minigene assays to functionally validate the genetic regulation mechanism of TRPT1. These results provide guidance for future investigations of drug targets for treating patients with neurological disorders with celltype resolution.

[0184] Embodiments of the present disclosure can be used as diagnostic tools, including being used as a companion diagnostic for a drug, for a range of diseases. For example, embodiments of this disclosure can be used to detect fusion genes from single-cell / single-nucleus RNA-seq data from cancer patients. For a second example, the invention can be used to detect mis-splicing that causes the spinal muscular dystrophy (SMD).

[0185] Embodiments of the present disclosure can be used to identify new target genes for complex diseases. Complex disease means a disease that is not caused by the mutation in a single gene, but whose risk can be influenced by many genes and environmental factors. Target gene means a gene whose mutation could lead to increased risk of complex diseases. For example, this disclosure has identified IRF5 as a key risk gene for rheumatoid arthritis and systemic lupus erythematosus. For a second example, this disclosure has identified TCHP as a key risk gene for Graves’ disease (autoimmune hyperthyroidism).

[0186] Whilst the foregoing description has described exemplary embodiments, it will be understood by those skilled in the art that many variations of the embodiments can be made within the scope and spirit of the present invention.

Claims

CLAIMS1. A sequence data analysis method comprising:receiving single-cell RNA sequence data comprising RNA sequence reads; aligning the RNA sequence reads to a reference genome;extracting intron junctions from the aligned RNA sequence reads; determining read counts for the extracted intron junctions; and identifying splicing events using the read counts for the extracted intron junctions.

2. The sequence data analysis method according to claim 1, wherein identifying splicing events using the read counts for the extracted intron junctions comprises identifying splicing events from site usage ratios of detectable splice sites.

3. The sequence data analysis method according to claim 2, further comprising performing splicing quantitative trait loci mapping using the splice site usage ratios.

4. The sequence data analysis method according to claim 1, wherein identifying splicing events using the read counts for the extracted intron junctions comprises identifying splicing events from intron excision ratios across the extracted intron junctions.

5. The sequence data analysis method according to claim 4, further comprising performing splicing quantitative trait loci mapping using the intron excision ratios across the extracted intron junctions.

6. The sequence data analysis method according claim 3 or claim 5, further comprising identifying disease-associated splicing events by assessing colocalization between the splicing quantitative trait loci and genome-wide association study loci.

7. The sequence data analysis method according claim 6, wherein the disease associated splicing event is associated with blood traits, anthropometric traits and / or immune-related traits.

8. The sequence data analysis method according to any preceding claim, wherein aligning the RNA sequence reads to a reference genome is carried out in two-pass mode.

9. The sequence data analysis method according to any preceding claim, further comprising detecting fusion genes from an identified splicing event.

10. The sequence data analysis method according to any one of claims 1 to 7, further comprising detecting a mis-splicing event associated with a disease.

11. The sequence data analysis method according to claim 10, wherein the disease is a rare disease that affects less than 65 out of 100 000 individuals such as spinal muscular dystrophy, Duchenne muscular dystrophy, and familial amyloid polyneuropathy.

12. The sequence data analysis method according to any one of claims 1 to 7, further comprising identifying a target gene having a mutation which leads to an increased risk of a disease, such as Graves’ disease, rheumatoid arthritis, asthma, atopic dermatitis, and systemic lupus erythematosus.

13. A sequence data processing system configured to carry out the method according to any one of claims 1 to 12.

14. A computer readable medium carrying computer executable instructions which when executed on a processor cause the processor to carry out a method according to any one of claims 1 to 12.

Images