CancerSplicingQTL: a database for genome-wide identification of splicing QTLs in human cancer

Abstract Alternative splicing (AS) is a widespread process that increases structural transcript variation and proteome diversity. Aberrant splicing patterns are frequently observed in cancer initiation, progress, prognosis and therapy. Increasing evidence has demonstrated that AS events could undergo modulation by genetic variants. The identification of splicing quantitative trait loci (sQTLs), genetic variants that affect AS events, might represent an important step toward fully understanding the contribution of genetic variants in disease development. However, no database has yet been developed to systematically analyze sQTLs across multiple cancer types. Using genotype data from The Cancer Genome Atlas and corresponding AS values calculated by TCGASpliceSeq, we developed a computational pipeline to identify sQTLs from 9 026 tumor samples in 33 cancer types. We totally identified 4 599 598 sQTLs across all cancer types. We further performed survival analyses and identified 17 072 sQTLs associated with patient overall survival times. Furthermore, using genome-wide association study (GWAS) catalog data, we identified 1 180 132 sQTLs overlapping with known GWAS linkage disequilibrium regions. Finally, we constructed a user-friendly database, CancerSplicingQTL (http://www.cancersplicingqtl-hust.com/) for users to conveniently browse, search and download data of interest. This database provides an informative sQTL resource for further characterizing the potential functional roles of SNPs that control transcript isoforms in human cancer.


INTRODUCTION
Single nucleotide polymorphisms (SNPs) are the most frequent genetic variants in humans and represent a valuable resource for investigating the genetic basis of diseases (1). Genome-wide association studies (GWAS) have found abundant SNPs associated with various traits and diseases, but most of risk loci lack clear molecular mechanisms (2,3). Expression quantitative trait locus (eQTL) studies have been employed to identify SNPs that may influence the expression levels of genes, thereby contributing to the phenotype outcome (4)(5)(6). However, only a moderate proportion of GWAS-identified loci are strong eQTLs (7), which might be partly due to the small sample sizes, the tissues studied, and a focus on overall gene level expression measurements without consideration of transcript isoforms (8).
Alternative splicing (AS) is a molecular mechanism that produces multiple distinct transcript isoforms from a single gene. The invention of RNA sequencing greatly facilitated the identification of AS on a genomic scale (9). In human, AS can occur in ∼90% of genes in a cell type-, conditionor species-specific manner, which is thought to extensively increase the number of proteins over the number of genes in a genome (10,11). In cancer, aberrant splicing patterns are frequently observed and known to contribute to carcinogenesis, de-differentiation and metastasis (12). Many D910 Nucleic Acids Research, 2019, Vol. 47, Database issue cancer-specific transcript isoforms have been identified (13). For example, an alternatively spliced transcript isoform of the gene encoding spleen tyrosine kinase is frequently expressed in breast cancer cells but never in matched normal tissues (14). Available evidence reveals that at least 20% of disease-causing single base-pair mutations affect splicing (15). Common genetic variation that affects splicing regulation, referred to as splicing quantitative trait loci (sQTLs), can lead to differences in alternative splicing between individuals, consequently influence disease susceptibility and drug response (16). Thus, the identification of sQTLs, especially in cancer tissues, might represent an important step toward fully understanding the contribution of genetic variants in tumorigenesis and development.
Because of the significance of sQTLs, several studies have performed genome-wide sQTL identifications on different human tissues, such as whole blood and brain (8,(17)(18)(19). These large-scale transcriptome studies using highthroughput genotyping method and deep RNA sequencing have revealed widespread sQTLs throughout the genome. However, no database comprehensively provides sQTLs for a large number of cancer samples. To bridge this gap, we have developed a computational pipeline to systematically identify sQTLs in 33 cancer types incorporating 9026 tumor samples from The Cancer Genome Atlas (TCGA). We identified millions of sQTLs across cancer types, and constructed a user-friendly database, CancerSplicingQTL (http://www.cancersplicingqtl-hust.com/) for users to conveniently browse, search and download data of interest.

Values of splicing events collection and processing
Percent Spliced In (PSI) values of each AS event were downloaded from the TCGASpliceSeq (http://projects.insilico.us. com/TCGASpliceSeq/PSIdownload.jsp) with default parameter (20). PSI value is a common, intuitive ratio for quantifying splicing events (11). The value is calculated by the transcript element present divided by the total number of reads covering the splicing event ( Figure 1A). Six types of AS events were analyzed in CancerSplicingQTL, including skipped exon, retained intron, alternative donor sites, alternative acceptor sites, alternate terminator and alternate promoter ( Figure 1B). For each cancer type, probes were filtered using the following criteria: (i) the rate of missing PSI value >0.1, (ii) mapping to locations on sex chromosome ( Figure 1E). Finally, an average of 34 942 AS events per cancer type were used for analyses. To minimize the effects of outliers on the regression scores (21)(22)(23), the values for each probe across samples per cancer type were transformed into a standard normal distribution based on rank.

Genotype data collection, imputation and processing
We downloaded genotype data (level 2) of 10 944 tumor samples from the TCGA data portal (https://portal. gdc.cancer.gov/), which detected the genotypes using Affymetrix SNP Array 6.0 containing 898 620 SNPs for each sample. Of these samples, 9026 samples were available with PSI data. To increase the power for sQTL discovery, we imputed autosomal variants for all samples in each cancer type using IMPUTE2, with 1000 Genomes Phase 3 as the reference panel as described in our previous study (24). To improve computation efficiency, we used the two-step procedure of IMPUTE2, which includes pre-phasing and the imputation of the phased data. Following criteria were used to exclude SNPs: (i) imputation confidence score, INFO < 0.4, (ii) minor allele frequency (MAF) < 5%, (iii) SNP missing rate ≥ 5% for best-guessed genotypes at posterior probability ≥ 0.9 and (iv) Hardy-Weinberg Equilibrium P-value < 1 × 10 −6 estimated by Hardy-Weinberg R package (Figure 1C). After imputation and quality filtering, an average of 4 516 897 genotypes per cancer type were remained in the sQTL analyses.

Covariates
In QTL analyses, covariates are often included to correct for the known and unknown confounders and increase the sensitivity of analyses (25). The top five principal components (PC) calculated by smartPCA in the EIGENSOFT program (26) were included to control for ethnicity differences, as they account for 10% of the variation explained with diminishing returns (0.5% or smaller contribution) for subsequent PCs, which are sufficient to represent the major population structure found in the TCGA dataset. Furthermore, to remove the hidden batch effects and other confounders from the AS data, we used PEER software (27) to infer hidden determinants, and selected the first 15 PEER factors from the AS data as covariates. The hidden batch effects and ethnic differences respectively accounted for an average of 19.9% and 1.19% of contribution to PSI variance in all cancers, which were described in details at the Supplementary Table S1. Other common confounders, specifically age, sex, and tumor stage, were included as additional covariates (22,28,29) ( Figure 1D).

Identification of survival-associated sQTLs
As many AS are involved in cancer prognosis (31), sQTLs may alter gene splicing and thereby influence the prognosis. To prioritize promising sQTLs, we additionally identified sQTLs that might be associated with patient survival times. For each sQTL, we examined the associations between the sQTL and patient overall survival times. For each sQTL, samples were classified into three groups: homozygous genotype AA, heterozygous genotype Aa and homozygous genotype aa (A and a represent two alleles of one SNP). The log-rank test was used to examine the differences in survival time, and Kaplan-Meier (KM) curves were plotted to represent the survival times for each group. sQTLs with FDR < 0.05 were defined as survival-associated sQTLs.

Identification of GWAS-associated sQTLs
The identification of causal variants is a major challenge for post GWAS studies (32). Thus we integrated the sQTLs with known GWAS risk loci to facilitate interpretation of the function of genomic variants. We downloaded all the known risk tag SNPs identified in GWAS stud- ies from the National Human Genome Research Institute (NHGRI) GWAS catalog (http://www.ebi.ac.uk/gwas/, accessed by 1 March 2018) (2). Then we obtained GWAS linkage disequilibrium (LD) regions of these risk tag SNPs from SNAP (https://personal.broadinstitute.org/plin/snap/ ldsearch.php) (33) with parameters (SNP data set: 1000 Genomes; r 2 (the square of the Pearson correlation coefficient of linkage disequilibrium) threshold: 0.5; population panel: CEU (Utah Residents with Northern and Western European Ancestry), and distance limit: 100 kb). sQTLs that overlapped with GWAS tag SNPs and LD SNPs were defined as GWAS-related sQTLs.

Samples in CancerSplicingQTL
In total, CancerSplicingQTL included 9026 tumor samples with both genotype data and PSI data available for 33 cancer types (Figure 2A). The sample size for each cancer type ranged from 36 for cholangiocarcinoma (CHOL) to 1090 for invasive breast carcinoma (BRCA) ( Table 1). After imputation and quality control of the genotype data, an average of 4 516 897 SNPs per each cancer type were used for analyses, ranging from 2746 175 for BRCA to 5 120 270 for acute myeloid leukemia (LAML). After removing AS events with a rate of missing PSI beta value > 0.1 or mapping to sex chromosome, an average of 34 942 splicing events per cancer type were used for analyses, ranging from 24 707 for uterine corpus endometrial carcinoma (UCEC) to 43 937 for esophageal carcinoma (ESCA).

sQTLs in CancerSplicingQTL
The CancerSplicingQTL mainly contains three datasets that are sQTLs, survival-sQTLs and GWAS-sQTLs. In the sQTL analysis, the associations between each SNP and AS events within the ±100 kb window around the SNP were analyzed for sQTL mapping by linear regression. We totally identified 7 945 857 sQTL-AS pairs at a per-tissue FDR < 0.05 in 33 cancer types. In total, there are 4 599 598 sQTLs across cancer types, ranging from 64 in CHOL to 574 577 in thyroid carcinoma (THCA), with a median of 124 542 sQTLs per cancer type ( Table 1). The number of sQTLs was significantly correlated with the number of samples (Spearman correlation Rs = 0.96, P-value = 1.44 × 10 −18 , Figure  2B). These sQTLs affect a median of 4847 AS events of 2857 unique genes per cancer type. Most of sQTLs were centered on AS events and 50% of sQTLs located at ±31 kb region flanking the AS events ( Figure 2C). 42.8% of sQTLs were associated with multiple AS events, and of these affected AS events, 33.6%, 27.3% and 18.7% were alternate terminator (AT), exon skip (ES) and alternate promoter (AP), respectively ( Figure 2D).
To prioritize promising sQTLs, we linked sQTLs to patient survival times and known GWAS loci. We found 17 072 sQTLs associated with patient overall survival times across different cancer types at FDR < 0.05. The number of survival-sQTLs ranged from 0 in CHOL and lymphoid neoplasm diffuse large B-cell lymphoma (DLBC) to 1643 in prostate adenocarcinoma (PRAD). We also linked sQTL results to NHGRI GWAS Catalog data and found 1 180 132 sQTLs that overlapped with GWAS linkage disequilibrium (LD) regions of one or multiple traits.

DATABASE ORGANIZATION AND WEB INTERFACE
CancerSplicingQTL was built based on the NodeJS 8.10.0 (https://nodejs.org/en/) framework with MongoDB 3.6.5 (https://www.mongodb.com/) as its database engine. It runs on a Linux-based Nginx Web server, while ReactJS (https: //reactjs.org/), a modern JavaScript library, is used for building user interfaces. We have tested it on Google Chrome (preferred), Firefox or Apple Safari browsers. The SplicingQTL website is available online (http://www. cancersplicingqtl-hust.com/) and requires no registration.
We provided a user-friendly web interface that facilitates searching, browsing and downloading the three datasets. Users can enter the 'sQTL/survival-sQTL/GWAS-sQTL' pages by clicking on the corresponding button in the browser bar ( Figure 3A) or on hyperlinks embedded in the corresponding images in the 'Modules' section on the 'home' page ( Figure 3C). Two query sections 'Single Search' and 'Batch Search' are provided for comprehensive queries across all three datasets ( Figure 3B). In the 'Single Search' section, users can select a specific cancer type (e.g. BRCA) and input an SNP ID (e.g. rs936227), gene symbol (e.g. ULK3) or genomic region (e.g. chr15:75 100 000-75 200 000) to search sQTLs across all datasets. If users do not select cancer type, it will return results for all cancer types. The 'Batch Search' section allows users to input multiple cancer types, SNPs, genes or genomic regions of interest. For instance, inputting 'rs936227' and 'rs9989230' in the 'SNP ID' box, will return a complete list of matched entries across cancer types. In addition, a summary of the sample size, sQTL number and AS type distribution is also  shown on the 'home' page. Putting the cursor over a cancer name on the-hand side human anatomy diagram, the matched results will show on the right-hand side figures. All data in the database can be downloaded from the 'Download' page. A detailed tutorial showing how the data were collected and processed is available on the 'Help' page. Can-cerSplicingQTL welcomes any feedback by email via the 'Contact' page.

Query on the 'sQTLs' page
To query sQTLs, CancerSplicingQTL allows users to search by selecting a cancer type from a pull-down menu, or by entering a SNP ID or gene symbol. After users click the 'Search' button, the query results are displayed in a table showing SNP ID, SNP genomic position, SNP alleles, related gene symbol, splicing type, splicing exon, splicing ID (the same as TCGASpliceSeq annotation), splice position, beta value (effect size of SNP on PSI value), r value (correlation coefficient) and P-value of sQTL ( Figure 3D). By clicking the hyperlink 'Box Plot' on the right of each record, a vector diagram of a boxplot will display the association between SNP genotypes and normalized PSI values. For example, our analysis showed that at XPNPEP3 first exon, the PSI values of individuals carrying the homozygote rs5758123 aa is significantly higher than that of individuals carrying the homozygote rs5758123 AA and heterozygous rs5758123 Aa in stomach adenocarcinoma (P-value = 1.34 × 10 −53 , Figure 3G).

Query on the 'survival-sQTLs' page
A table with SNP ID, SNP genomic position, SNP alleles, Log-rank test P-value, and median survival time for each genotype group is displayed on the survival-sQTLs page ( Figure 3E). Search boxes are designed to retrieve specific cancer types and SNPs. If users select a specific cancer type or input a gene or SNP ID, the table will be reconstructed to display the results of the query. Each record embeds a hyperlink 'KM Plot', showing the association between SNP genotypes and overall survival times. For example, our analysis showed that patients with rs7597845 AA allele have a better prognosis than other patients with breast cancer (Pvalue = 9.09 × 10 −6 , Figure 3F).

Query on the 'GWAS-sQTLs' page
A complete list of the SNP information, regulated splice site, related gene information and related GWAS-traits are provided on the 'GWAS-sQTL' page. Search boxes are designed to retrieve a specific cancer type, phenotype or SNP. In addition, users can select a different LD threshold from the 'LD' dropdown box to prioritize SNPs. For example, the GWAS-catalog has collected 263 tag SNPs of breast cancer risk loci. We found that 83 tag SNPs have 1402 sQTLs in their LD regions (r 2 ≥ 0.5) affecting the splicing events of 100 genes. Causal variants of breast cancer could be existed among these sQTLs.

SUMMARY AND FUTURE DIRECTIONS
In summary, CancerSplicingQTL is a comprehensive sQTL resource that uses large cancer samples to evaluate the ef-fects of genetic variants on gene splicing. It provides a userfriendly interface for users to query, browse, and download sQTLs. To the best of our knowledge, CancerSplic-ingQTL is the first public database focusing on cancerspecific sQTLs. Millions of vector diagrams of sQTL box plots and KM plots are provided for scientific usage. We also identified numerous sQTLs associated with patient survival times or located in known GWAS loci that will be promising candidates for genetic research. Biologists can download entire datasets for further integrative studies.
Cancer genomics studies are developing rapidly (34,35), and we expect the number of cancer samples with genotype and splicing profiles to increase dramatically. In the future, we will continue to update CancerSplicingQTL to include more cancer samples and maintain it as a useful resource for the research community. We will add more genetic and splicing information into the database. We believe that Cancer-SplicingQTL will be an important resource for human cancer genetics, providing opportunities to bridge the knowledge gap from variants in sequence to phenotypes.