Exhaustive capture of bladder cancer transcriptional and genomic variants integrating canonical and mapping-free protocols

Data extraction

The Binary Alignment Map (BAM) files of all BLCA patients were retrieved from the TCGA database with permission from dbGAP [12]. BAM files were converted to paired-ends fastq files using Picard tools with default parameters. Both of the RNAseq and WES sequencing data were involved in this study. All patients' clinical information was also obtained, including the survival time, stage, relapse, or metastasis. Patients were divided into two groups: complete remission and progressive group. Adapters sequences were trimmed using the cutadapt software with default parameters [13]. Duplicated reads were removed since these reads are majority generated from the PCR process instead of natural status in cells.

Canonical protocol based on genome reference

The latest version of the genome and annotation files were downloaded from the Gencode [14]. Reads from RNAseq data were mapped to the human genome of the hg38 version using the STAR software with default parameters [15]. Reads from the WES data were mapped to the human genome using the BWA algorithm with default parameters [16]. The BAM file processing is performed using the Cufflinks software [17]. Given genome annotation file (GTF) and BAM files, we performed differential expression analysis of genes using the Cuffdiff algorithm [18]. Cuffdiff calculates the FPKM of each gene in each sample by summing the FPKMs of all transcripts belonging to the host gene. Differentially expressed genes were selected by Cuffdiff with default parameters from the FPKM matrices. For WES data analysis, GATK was applied to call somatic variants using the BAM files from WES dataset, and all somatic variants were stored in the Variant Call Format (VCF) format files [19]. The single-nucleotide variants (SNVs) from all BLCA patients were aggregated as a matrix. SNVs observed in at least two patients were selected as recurrent SNVs.

Mapping-free protocol without genome reference

As the human genome is not completed and gaps are present, all variations within the 'unannotated' regions are not captured by canonical methods. Herein, we introduced a mapping-free protocol named DEkupl [9]. The reads are decomposed to k-mers, and all the k-mers different to reference are retained. In this way, we can exhaustively capture all the variations in the “dark genome regions”. The DEkupl estimates the gene expression using the Kallisto, which is also a reference-free software [20]. The expression level of both genes and k-mers were compared between the complete remission and progressive groups. Genes and k-mers with P values less than 0.05 and log2FC beyond 95% confidence interval were considered differentially expressed genes (DEGs). Consequently, the DEkupl software provides both differentially expressed genes detected by limma and novel events with variations exhibited as contigs.

Gene-level candidates through transcriptome analysis

As the gene expression estimators in the two protocols are based on different rationales, we compared the differential genes detected by both protocols between the complete remission and progressive groups. Cufflink and Kallisto were applied to detect DEGs in the mapping-based and mapping-free protocols, respectively. Only the consistent differential genes between the two protocols were considered to be bladder cancer-related candidates. Eventually, the consistent DEGs between responsive and progressive patients were selected as gene-level candidates.

Contig-level candidates through transcriptome analysis

Since we do not align reads to the reference in the mapping-free protocol, we retained all the contigs harboring differences to the reference. We can benefit from this because variations within the repeat or low complexity regions can also be kept, which would be ignored or discarded by the aligners such as STAR or BWA. Only the differentially expressed contigs (DECs) were retained and further mapped to the genome for annotation. The annotation process was done using GSNAP of version 2020-06-04 [21]. The DECs contain multiple events, including SNV, splice, split, lincRNA, polyA, repeat, and unmapped contigs. The unmapped contigs may either come from the exogenous microorganisms or unannotated human genome. All these events represented by contigs are absent in the transcriptome, thus the mapping-free protocol products were regarded as 'novel' events.

Convincing SNVs integrating WES and RNAseq analysis

In the WES data analysis, we screened recurrent SNVs present in at least two bladder cancer patients. However, many such SNVs do not make sense if their host genes do not show differential expressions or even stay silent. Only SNVs embedded in the DEGs are more likely biological factors to cancer progression. Therefore, we compared the host genes of all recurrent SNVs to the consistent DEGs (gene-level candidates) and only kept the intersection genes.

Survival analysis of genes and contigs

The convincing candidates were composed of consistent DEGs between two protocols and significant DECs detected by the mapping-free protocol. The candidates were compared to the patients' overall survival time using survival analysis. The P values were calculated using the log-rank test [22] with univariate cox regression. The Kaplan–Meier curves [23] were drawn for the top prognosis related candidates. Besides the univariate cox regression, we also applied the cox proportional-hazards model (coxPh), a multivariate regression model [24]. Genes and contigs significant in both univariate and multivariate cox analysis were selected as prognostic candidates.

Deeply investigate the novel events

The novel events detected by mapping-free protocol cannot be mapped to the reference by aligners. We cannot locate these novel events to specific genes or positions on the genome either. However, the contigs corresponding to these novel events show differential expression as cancer progresses. We compared signature contigs' expression profile between the early stage (stage I/II) and late-stage (stage III/IV) patients to deeply investigate the relationship between the novel events and cancer progression. Survival analysis was performed to select the potential signature contigs as prognostic indicators.

Diagnostic model construction integrating the signatures

The significant gene candidates were screened using the log-rank test. To take the collinearity between variables into account, we further applied the lasso regression [25] to select diagnostic signatures. The samples were randomly separated into two groups, in which 75% were training set and the rest 25% were test set. A support vector machine (SVM) model was constructed using the selected signatures [26]. Five-fold cross-validation [27] was used to assess the performance of the diagnostic model on the trainset. The samples were randomly split into five folds. We trained the model using four folds of samples in each cycle and tested it using the rest one fold. This procedure was repeated until all samples were predicted once as test sample. The ROC curve was shown with the sensitivity and specificity on the test set [28].

Validation on the independent GEO dataset

To further validate our model's performance of predicting progressive bladder cancer patients, we applied another independent dataset from the GEO database. The validation data accession is GSE31189, which comprises 52 bladder cancer patients and 40 normal samples. The performance of our diagnostic model was shown using the ROC curve.

Differentially expressed genes from RNAseq data

Differentially expressed genes (DEGs) between BLCA patients and normal samples were extracted from the RNAseq data integrating mapping-based and mapping-free protocols. Up/down-regulated genes were extracted according to the P values and fold change values. The selected genes present diverse expression levels between responsive and progressive patients. The volcano graphs of DEGs consistent between canonical RNAseq pipeline and DEkupl were drawn, as shown in Figure 1A.

Figure 1B shows the consistent DEGs between the two protocols. The consistent DEGs were used as stable prognostic risk genes. Eventually, we obtained 1240 consistent DEGs in total, including 727 up-regulated and 513 down-regulated genes, respectively.

To check if two groups of patients present differences on DEGs level, we performed the principal component analysis (PCA) [29]. As seen in Figure 1C, the two groups of patients can be distinguished using the first three components.

To further clarify the functions regulated by the DEGs, we enriched the DEGs to both gene ontology (GO) and the KEGG pathway [30] using the clusterProfiler package [31]. The top 10 GO biology process terms are shown in Figure 2A. The significant functions include axonogenesis and morphogenesis. Figure 2B shows the four hub functions and corresponding genes. These functions are supposed to be activated as most of the involved genes were up-regulated in progressive patients. Figure 2C suggests that some enriched functions share cross-talk genes. These genes achieve the communication between functions as bridges. Figure 2D is the most significant KEGG pathway that was enriched by 23 DEGs. It implies the genetic association between Bladder cancer and melanoma.

Recurrent SNVs from WES data

Large numbers of mutated genes do not express in the Bladder tissues. The significant number of SNVs is attributed to the gene length, for instance, TTN [32]. Therefore, the frequently mutated genes have little impact on the disease progression if they are not expressed. Given this hypothesis, we only investigated the SNVs from the differentially expressed genes. The most recurrently mutated genes were identified using the GATK across all patients. The frequency and patients harboring each mutated gene can be seen from Figure 3A.

As seen in Figure 3A, the top 20 differentially expressed genes with the highest mutational frequency can be observed in at least 5% of patients, in which the highest frequency is 7%. Interestingly, the most well-known cancer suppressor gene, TP53 [33], is absent in our top 20 gene list, which has a mutational frequency of 88% in bladder cancer patients. The reason is TP53 is not differentially expressed in progressive bladder cancer patients despite a high mutational frequency. In contrast, the genes in our list were recurrently mutated in bladder cancer patients and presented significant differences as the tumor progressed (Table 1).

To further investigate the internal relationship between the most mutated genes, we applied the maftools [34] to estimate the interactions between genes. There are two types of interactions between genes that are the co-occurrence and mutually exclusive relationship. Co-occurrence interaction indicates that two genes tend to be mutated simultaneously [35]. The mutually exclusive interaction suggests that two genes are barely observed to be mutated simultaneously [36]. It's worth mentioning that mutually exclusive interaction does not mean two genes could not be mutated simultaneously. Instead, it means the cells carrying both of the two mutated genes tend to trigger the apoptosis program and turn out to die. This phenomenon is called synthetic lethality [37], which is a significant avenue of cancer therapy.

As seen in Figure 3B, there are many mutually exclusive gene pairs whose P values are less than 0.1, marked as dots. The predominant interaction among the top 20 genes is co-occurrence. The gene pairs with co-occurrence relationships are marked with stars, suggesting significant P values less than 0.05.

Novel events identified by DEkupl

Since we also applied a mapping-free protocol, we could capture novel events unsolvable to the mapping-based methods. DEkupl captured all the k-mers absent in the human reference and then merged the k-mers into contigs. These contigs were therefore considered to be the context of local events. We eventually identified contigs composed of 139 SNVs, 412 intron events, 6 repeats, and 24 unmapped events. As DEkupl masked all the k-mers present in the reference, these contigs extended from the retained k-mers can be considered novel events. The detailed results can be seen from the Supplementary Table 1. All these novel events present significant differences between the two groups (Figure 4A).

In Figure 4A, most of the novel contigs were highly expressed in the progressive group. Meanwhile, we observed that some complete remission patients exhibit similar expression patterns as the progressive patients. It suggests that patients sensitive to therapy at the initial phase might carry a risk of relapse in the following phases.

Candidates SNVs screening

Besides the shared DEGs, we also screened the SNVs detected by both of the two protocols from WES data as convincing SNVs. The contigs generated by DEkupl were mapped to the human genome using GSNAP software. Then the genomic coordinates were compared with the GATK results. Finally, we obtained 450 convincing SNVs belonging to the DECs (Figure 4B). Figure 4B shows the locations of SNVs and the correlation between the host genes. Some SNVs are found to locate at closed loci forming a local cluster. And genes from the same chromosome tend to present a positive correlation. In contrast, genes from distant genomic loci tend to have a negative correlation.

Prognosis related indicators

Since the SNVs whose host genes belong to the DEGs were considered potential indicators, we intended to extract prognosis related indicators from the whole DEGs list. We initially used the log-rank test to estimate each gene. The patients were divided into two groups according to the median expression value of the queried gene. Then the log-rank test P value was computed. After ranking the genes P values in ascending order, we selected the top 6 genes and drew the KM curves, as shown in Figure 5.

As seen from Figure 5, patients with different expression levels of the six genes tend to present various surviving time. The high expression levels of the gene FABP4 and PCDHGA7 always imply a poor surviving prognosis. In contrast, the low expression levels of the gene PRSS30P, LOC100272228, KLRC2, and SLC10A5 indicate a poor surviving prognosis.

We also applied the coxPh, a multivariable regression method, to select prognostic indicators from the DEGs. The top 20 genes can be seen in Table 2.

The top indicator lists in the coxPh and log-rank test are not the same. This is because coxPh considers the internal interactions among genes, while the log-rank test treats each gene independently. However, among the 1240 DEGs, 4 of the top 6 genes detected by the log-rank test rank in the top25%. The ranks of PRSS30P, LOC100272228, PCDHGA7, and SLC10A5 are 79, 124, 310, and 43, respectively.

Survival analysis using the contigs of novel events

Besides DEGs, we also captured some novel events from the RNAseq data. These novel events were composed of SNVs, intron, splice, repeat, and unmapped contigs. It's worth mentioning that the unmapped contigs may come from transcripts produced by rearranged genes or result from exogenous viral genomes and could thus be highly relevant biologically. We therefore investigated the correlation between surviving time and the novel events, including unmapped contigs.

We compared the expression of contigs corresponding to the novel events, as shown in Figure 6A. The violin graph indicates the transcriptional difference between the early and late-stage patients. Survival analysis also proves the correlation between these novel events and survival prognosis. Interestingly, we showed the top 6 most significant contigs in Figure 6B. All the patients with high expressions of these 6 contigs tend to have a worse prognosis than the others. Our findings imply that the overexpression of these unannotated events may promote cancer progression.

Diagnostic model construction

To screen the diagnostic signatures, we herein applied a machine learning feature selection algorithm called lasso regression. We combined all the DEGs and convincing SNVs as features. Lasso regression determined the best combination of features according to the Log gamma criteria.

In Figure 7A, 17 signatures were selected by lasso regression via the cutoff of log gamma. The best log gamma cutoff was determined automatically by lasso regression. The coefficient of each signature can be seen in the Figure 7B. Positive and negative signatures were marked in blue and red, respectively. Finally, the SVM model was trained using the 17 signatures, and the AUC was 90% on the test set (Figure 8A). Our results demonstrate that the signatures in our findings can distinguish progressive patients from responsive patients. Many of the 17 genes are already verified by previous studies to be cancer related. Thorsen K found that ACSL5 shows tumor-specific alternative TSS usage in cancer samples [38]. The gene of ACSS3 is proved to promote BLCA cell growth [39]. ART3 and ASAH2B are two well-known cancer driver genes, especially in breast cancer patients [40–43]. AZGP1 is involved in colorectal cancer and androgen receptor-induced prostate cancer [44, 45]. DENND2D, PCDHB7, FAM113B, and GBX2 have prognostic impact on various cancers, including lung cancer, gastric cancer, ovarian and esophageal cancer [46–50]. TLR4 and NR1H3 are the prognostic factors of bladder cancer patients [51–53].

Validation and comparison

We trained a diagnostic model using the signatures and proved the performance of predicting progressive bladder cancer patients. To further investigate if our model can also predict bladder cancer patients from the normal samples, we applied another independent data from the GEO database. The validation data comprises 52 bladder cancer patients and 40 normal samples. The ROC curve on the validation data is shown in Figure 8B. The accuracy of predicting bladder cancer patients is 0.87, which is slightly lower than the discovery data.

We also compared our 17-signatures model to the other models published previously. The first model is an immune prognostic model (IPM) developed by Hongyan Li [54]. The IPM is based on immune-associated genes derived from differential analysis between wild-type TP53 and mutated TP53 bladder cancer patients. The IPM was validated using an independent cohort from the GEO database. IPM utilizes the expression of four essential genes (CTSG, TREML4, KRT1, and PPBP) to calculate the risk score for each patient. The other model we compared to is an autophagy prognostic model (APM) developed by WANG et al. [55]. The APM was constructed with autophagy-related genes (ARGs). The author screened three key prognostic ARGs (JUN, MYC, and ITGA3) to calculate each patient's risk score. The APM was validated by two GEO datasets (GSE13507 and GSE31684). The comparison between our model and the other two models is shown in Figure 9.

As shown in Figure 9, our 17-signature model outperforms the other two models in both TCGA and GEO cohort. APM and IPM apply four and three essential genes as signatures, respectively. From this point, our 17-signature model takes more genes into account and thus achieves better performance expectedly. On the other hand, APM and IPM focus on either immune related genes or autophagy related genes. In contrast, we considered the whole transcriptome and exome. In addition to biological interpretation, our model also has an advantage in terms of accuracy.