RNA-Seq accurately identifies cancer biomarker signatures to distinguish tissue of origin.

Metastatic cancer of unknown primary (CUP) accounts for up to 5% of all new cancer cases, with a 5-year survival rate of only 10%. Accurate identification of tissue of origin would allow for directed, personalized therapies to improve clinical outcomes. Our objective was to use transcriptome sequencing (RNA-Seq) to identify lineage-specific biomarker signatures for the cancer types that most commonly metastasize as CUP (colorectum, kidney, liver, lung, ovary, pancreas, prostate, and stomach). RNA-Seq data of 17,471 transcripts from a total of 3,244 cancer samples across 26 different tissue types were compiled from in-house sequencing data and publically available International Cancer Genome Consortium and The Cancer Genome Atlas datasets. Robust cancer biomarker signatures were extracted using a 10-fold cross-validation method of log transformation, quantile normalization, transcript ranking by area under the receiver operating characteristic curve, and stepwise logistic regression. The entire algorithm was then repeated with a new set of randomly generated training and test sets, yielding highly concordant biomarker signatures. External validation of the cancer-specific signatures yielded high sensitivity (92.0% ± 3.15%; mean ± standard deviation) and specificity (97.7% ± 2.99%) for each cancer biomarker signature. The overall performance of this RNA-Seq biomarker-generating algorithm yielded an accuracy of 90.5%. In conclusion, we demonstrate a computational model for producing highly sensitive and specific cancer biomarker signatures from RNA-Seq data, generating signatures for the top eight cancer types responsible for CUP to accurately identify tumor origin.

Introduction

Metastatic cancer of unknown primary origin (CUP) is an important clinical dilemma, comprising 3% to 5% of all new cancer cases [1], [2]. Without a firm histologic diagnosis, the clinical management of these patients varies widely [3], and despite protocol-driven guidelines, outcomes remain poor. Median survival is 6 to 9 months [4], with a 5-year survival rate of only 10% [5], [6]. The role of chemotherapy in the treatment of occult primary tumors is primarily palliative and does not improve long-term survival; the National Comprehensive Cancer Network (NCCN) panel encourages CUP patient enrollment in clinical trials when possible [7]. An accurate method for distinguishing tumor origin to tailor personalized therapies is therefore critical to the successful management of these malignancies.

Thus far, there have been a number of studies focused on identifying unique signatures that distinguish among different cancer types, using immunohistochemistry [8], [9], [10], [11], cytogenetic studies [12], [13], [14], comparative microarray analysis [15], [16], [17], [18], [19], [20], [21], [22], combined microarray and quantitative polymerase chain reaction techniques [23], [24], bead-based miRNA profiling [25], and more recently, limited, high-throughput sequencing data combined with microarray [26]. Currently, only qRT-PCR [23], [24] and microarray-based [19], [22] assays are commercially available for use, with diagnostic accuracies ranging from 74 – 85%.

Compared to traditional microarray technology, transcriptome sequencing (RNA-Seq) possesses a number of advantages, including unlimited genome coverage and discovery potential, a greater than 8000-fold dynamic range for quantifying gene expression levels, and the ability to identify splice variants, unmapped genes, and unrecognized non-coding RNAs [27], [28]. The rapidly decreasing cost of high-throughput sequencing methods has also improved the accessibility of these techniques for clinical application and allowed for the generation of large-scale datasets to robustly interrogate such clinical problems as CUP.

In a review of all published autopsies performed on CUP patients who died from cancer progression from 1944 to 2000, a primary tumor was successfully identified post-mortem in 73% of the 884 cases [29]. The most common tissues of origin were lung (27%), pancreas (24%), kidney (6%), colorectum (6%), stomach (5%), liver (5%), ovary (3%), and prostate (3%) [4], [30], [31], [32], [33], [34], [35], [36], [37], [38]. Our objective was therefore to identify lineage-specific biomarker signatures for each of these cancers, using a large, multi-cancer RNA-Seq database to distinguish tissues of origin from among different cancer types.

Material and Methods

Multi-Cancer RNA-Seq Database

Paired-end RNA-Seq data for 364 cancer samples from 22 different tissue types were used to compile a multi-cancer gene expression dataset as previously described [39]. This dataset was then supplemented with publically available RNA-Seq cancer data accessed from the International Cancer Genome Consortium [40] and The Cancer Genome Atlas [41]. This included four additional cancer types (acute myeloid leukemia, endometrial cancer, head and neck squamous cancer, and lung cancer). The dataset was restricted to those transcripts commonly annotated across all three datasets. The final composite data matrix was comprised of gene expression readouts for 17,471 transcripts from 3,244 cancer samples (139 cancer cell lines and 3,105 tissues) across 26 cancer types.

Model for Deriving an Optimal Biomarker Signature

R language [42] was used to program an algorithm to derive an optimal biomarker signature for a cancer type of interest (Figure 1). The 3,244 samples were randomly allocated to the training or test sets (Table 1). A maximum of 50 samples per tissue type were assigned to the training set (688 samples) and the remainder to the test set (2,556 samples). The training set was used to generate the optimal biomarker signatures for each cancer type, while the test set was reserved for final, external validation. Biomarker signatures were generated for the eight tissue types that account for approximately 80% of CUP cases (colorectum, kidney, liver, lung, ovary, pancreas, prostate, and stomach).

Figure 1

Algorithm for extracting optimal cancer biomarker signatures from RNA-Seq dataset. AUC, area under the receiver operating characteristic curve.

Table 1

Allocation of Cancer Samples to RNA-Seq Training and Test Sets

Cancer Type	All	Training Set	Test Set
Adrenal gland	3	3	0
Acute myeloid leukemia	174	50	124
Bladder	70	50	20
Breast	864	50	814
Cervix	8	8	0
Colorectum	244	50	194
Endometrium	333	50	283
Germ cell	1	1	0
Kidney	24	24	0
Liver	15	15	0
Head and neck	263	50	213
Lung	348	50	298
Lymphoma	11	11	0
Medulloblastoma	1	1	0
Melanoma	136	50	86
Merkel cell	3	3	0
Myeloproliferative neoplasm	9	9	0
Neuroblastoma	2	2	0
Neuroepithelioma	1	1	0
Oropharynx	4	4	0
Ovary	418	50	368
Pancreas	76	50	26
Prostate	154	50	104
Rhabdomyosarcoma	1	1	0
Salivary gland	4	4	0
Stomach	77	50	27
Total	3244	688	2556

Transcript normalization

Transcript reads were normalized with log transformation followed by quantile normalization to account for variations between and within datasets, such as differences in the amount of starting material and reported transcript units. The entire 688-sample training set was then divided into 10 randomly generated subsets, each with an equal proportion of samples of the cancer type of interest. A 10-fold cross-validation method was used to train the model on 9-fold and test each signature on the remaining 1-fold.

Univariate transcript analysis

Within each 9-fold training subset, area under the receiver operating characteristic (ROC) curve values (AUC) [43] were calculated for each of the 17,471 transcripts. The transcripts were then sorted by decreasing AUC.

Stepwise logistic regression

To generate an optimal signature, an iterative approach was used, rather than inputting all 17,471 transcripts into the model, which would increase computational burden and the likelihood of overfitting. The top 100 transcripts, based on univariate AUC rank, were introduced into a stepwise logistic regression model, to determine the optimal signature at each iteration for n, between 1 and 100 input transcripts. Logistic regression was performed in both directions to optimize the Akaike information criterion (AIC) [44], [45], [46] so that at each step, it was calculated whether the current signature would be improved not only by adding the next variable but also by discarding any of the variables present within the currently optimized signature. The final signatures were used to calculate the predicted likelihood of each sample in the remaining 1-fold being of that cancer type, given n input transcripts.

Biomarker signature selection

A final “cross-validated AUC” was determined for each signature generated from n transcripts, based on the calculated predictions for each sample compared to their true identities. The optimal biomarker signature was determined to be the one generated from the top n* number of transcripts that yielded the highest cross-validated AUC. The entire 688-sample training set was then used as the input training set to generate a final, optimal biomarker signature based on the top n* transcripts.

Internal validation

Each cancer biomarker signature was internally validated by using the entire 688-sample training set as the input. Each sample received a predicted value, m, between 0 and 1, indicating likelihood of the sample being the cancer type of interest. The predicted values were then used to generate ROC curves for each signature. Optimal score thresholds, k, (above which was defined as “positive” for that cancer type and below which was “negative”) were calculated by selecting the point on the ROC curve with the minimum distance from (0,1), which represents a perfect test of 100% sensitivity and specificity [47].

External validation

Each cancer biomarker signature was then externally validated against the reserved 2,556-sample test set using the optimal score thresholds. Overall sensitivity and specificity were calculated for each cancer signature.

Duplicate cancer predictions

Each of the 2,556 cancer samples in the reserved test set was tested using each of the eight cancer biomarker signatures. Those samples that predicted positive for more than one cancer type were assigned the cancer type that had the highest relative predicted value, defined as [m – k]/[1 – k].

Additional analysis

Graphs were plotted using GraphPad Prism. The heat map was generated with Cluster 3.0 [48] and visualized using TreeView [49]. Statistical analysis was performed using R and GraphPad Prism.

Results

The results of our biomarker-generating model are shown in Figure 2. For all eight cancer types, the maximum, cross-validated AUC was obtained within the first 100 input transcripts. Interestingly, cross-validated AUC plots of the colorectal, lung, pancreas, and stomach cancer samples yielded prominent peaks, beyond which the inclusion of additional transcripts worsened the biomarker signature’s accuracy. Conversely, liver and ovarian cancer samples yielded relatively flat curves of near-perfect cross-validated AUC’s, suggesting that these cancer types have such unique biomarker profiles that many highly accurate signatures may be generated. Optimal signatures for each tissue type were objectively determined by selecting the number of input transcripts, n*, that corresponded with the maximum cross-validated AUC (Figure 2, red points).

Figure 2

Scatter plots of 10-fold cross-validation method to determine optimal biomarker signatures for 8 different cancer types. Points highlighted in red indicate the highest, cross-validated AUC for each cancer type. AUC, area under the receiver operating characteristic curve.

Next, using the entire 688-sample training set as the input test set, the final list of transcripts was generated for each cancer biomarker signature, by performing stepwise logistic regression of the top n* transcripts. The entire model was then run again using a new, random allocation of training and test samples. The final biomarker signatures for each cancer type were concordant with the signatures generated from the first randomization (Table 2), with an overall cosine similarity measurement of 0.53 [50].

Table 2

Cancer Biomarker Signatures Generated from Two Separate Randomizations

Transcripts highlighted in red are common between the two signatures.

ROC curves were then generated for each biomarker signature, yielding high AUC’s (Figure 3); for comparison, lines of identity are shown, representing a random test that has no prognostic value. From each cancer signature ROC curve, threshold cut-offs minimizing the distance to (0,1) were calculated to use in subsequent external validation testing (Figure 3, red points).

Figure 3

Internal validation of eight cancer-specific biomarker signatures yields high area under the ROC curve values. The entire 688-sample RNA-Seq training set was used as the test set for each cancer signature. Dotted lines indicate lines of identity. Points of minimum distance to (0,1) are highlighted in red. ROC, receiver operating characteristic; AUC, area under the ROC curve.

The cancer biomarker signatures were then externally validated using the reserved 2,556-sample RNA-Seq test set. Each sample was tested against each of the 8 biomarker signatures and predicted to be positive or negative for that cancer type based on the threshold cut-offs calculated previously. The true identities of each of the 2,556 samples were then compared to the predicted identities, and high sensitivity (92.0% ± 3.15%; mean ± standard deviation) and specificity (97.7% ± 2.99%) were demonstrated for each biomarker signature (Table 3). A heat map representation of the 43 transcripts comprising the 8 biomarker signatures for all 3,244 samples illustrated the strength with which each cancer signature successfully distinguished samples from among the different cancer types (Figure 4). Despite similarities in expression patterns across different cancers, the transcript signatures were nonetheless highly specific, requiring elevated expression across all the biomarker transcripts for a given cancer type (Table 3). Of the 3,244 samples, 2.9% had positive predictions for more than one cancer. These samples with duplicate predictions were assigned the final identity of the cancer type with the highest relative predicted value, yielding an overall accuracy of 90.5% (Table 4).

Table 3

External Validation of Eight Cancer-Specific Biomarker Signatures Using 2,556-Sample RNA-Seq Test Set

	Colorectum	Kidney	Liver	Lung	Ovary	Pancreas	Prostate	Stomach
TP	174	0	0	274	360	24	95	24
TN	2355	2533	2552	2222	2188	2389	2424	2305
FP	7	23	4	36	0	141	28	224
FN	20	0	0	24	8	2	9	3
Sensitivity (%)	89.7	-	-	91.9	97.8	92.3	91.3	88.9
Specificity (%)	99.7	99.1	99.8	98.4	100	94.4	98.9	91.1

TP, true positive; TN, true negative; FP, false positive; FN, false negative.

Figure 4

RNA-Seq heat map of 8 cancer-specific biomarker signatures (rows) across all 3,244 cancer samples (columns).

Table 4

Overall Performance of RNA-Seq Biomarker Generating Algorithm for Predicting Tissue of Origin in 2,556 Cancer Samples

	Samples		%
TP	946	Sensitivity	95.0
TN	1366	Specificity	87.6
FP	194	PPV	83.0
FN	50	NPV	96.5
		Accuracy	90.5

Samples with duplicate cancer predictions were assigned the identity with the highest predicted value. TP, true positive; TN, true negative; FP, false positive; FN, false negative; PPV, positive predictive value; NPV, negative predictive value.

Discussion

In this study, we demonstrate an effective and efficient model for extracting highly sensitive and specific cancer biomarker signatures from a large RNA-Seq dataset. This technique yielded transcript signatures for the top eight cancer types that cause metastatic CUP. The robustness of the final signatures was demonstrated by external validation, through testing of a large test set randomly allocated a priori. By reserving these 2,556 samples solely for external validation, this test set represents a large dataset of “clinical unknowns” (i.e., identities unknown to the training model). The high sensitivities and specificities achieved with each of the cancer biomarker signatures therefore represent a realistic approximation of the accuracy with which the signatures may be used to predict tissue origin of a CUP sample in the clinical setting.

Many of the cancer biomarkers identified in this study have previously been well-characterized in their respective cancer types. The prostate cancer signature includes KLK3, which is responsible for encoding prostate-specific antigen, the serine protease used as a serum marker in prostate cancer screening and disease monitoring [51], as well as PRAC, which is known to be highly expressed in prostate cancers [52], [53]. Similarly, NKX3-1 is an androgen-regulated homeobox gene, which transcriptionally regulates oxidative damage response and is required for prostate stem cell maintenance; aberrant expression has been found to correlate strongly with prostate cancer progression [54], [55], [56], [57], [58], [59]. In addition, a recent study demonstrated that immunohistochemical staining with the kidney biomarker FXYD2, a Na-K-ATPase regulator, is highly sensitive and specific for renal cell carcinoma [60]. Similarly, NOX1 was highly expressed in our colorectal cancer samples, as confirmed in prior studies, which stimulates mitogenesis and angiogenesis though a ROS-mediated mechanism; NOX1 expression has also been found to correlate strongly with activating KRAS mutations, which are present in approximately 50% of colorectal tumors [61], [62]. IGFBP1, is a hepatocyte-derived secreted protein required for normal liver regeneration by inhibiting proapoptotic signals [63], with overexpression previously identified in hepatocellular cancers [64], as well as in our study. The lung biomarker NAPSA is a well characterized proteinase expressed in type II pneumocytes [65], [66], [67], whose expression has high sensitivity and specificity for distinguishing primary lung adenocarcinoma from metastatic pulmonary lesions from other primaries [67]. The pancreas biomarker NKX6-1 is a transcriptional regulator that has been shown to play an important role in beta cell differentiation during pancreatic development [68], [69], [70], [71]. Similarly, GKN1 is highly expressed in the gastric epithelium, providing protection to the antral mucosa and promoting healing after injury; it also acts as a tumor suppressor and is down-regulated compared to normal gastric tissues [72], [73], [74], [75], [76] but in our model was still significantly overexpressed as compared to other cancer types.

Our cancer signatures also identified transcripts that have not previously been associated with the cancer types of interest. Although not yet characterized in ovarian cancer, BEST1, which forms calcium-activated chloride channels across epithelial cells to promote cell proliferation [77], has been shown to be up-regulated in colon cancers [78]. Similarly, hypermethylation of DPYS, a gene important in pyrimidine metabolism, has been identified in prostate, colon, and breast cancers, as well as melanomas [79], [80], [81], but in our study, high expression was most sensitive and specific for kidney cancer. Interestingly, the transmembrane glycoprotein SI was highly expressed in both prostate and gastric cancer samples. Mutations in SI have previously been identified in head and neck, colorectal, and ovarian cancers, and in a recent study, SI mutations resulted in significant gene enrichment in oxidative phosphorylation, glycolysis/gluconeogenesis, and B-cell receptor signaling pathways, for promoting malignant progression in chronic lymphocytic leukemia [82].

In our pancreatic cancer signature, our model also identified several gene fusions not previously associated with this disease, including ANKHD1-EIF4EBP3, a readthrough transcript of the neighboring cell survival scaffolding gene ANKHD1 and the downstream translational repressor EIF4EBP3, both of which are effectors of the RAS/MAPK pathway [83], [84], which is known to play a critical role in the development and progression of pancreatic cancer [85], [86], [87], [88]. The prostate biomarker NKX3-1 and related family member NKX3-2 also comprised the pancreas signature. While the role of NKX3-2 in pancreatic cancer has not yet been characterized, its role in chondreogenesis and skeletal development has been well studied, acting as a transcriptional repressor downstream of SHH through interactions with the signal transduction protein SMAD4 [89], [90], [91], [92]. SHH and its related hedgehog-signaling pathways are well-known mediators of pancreatic carcinogenesis and are the targets of many new therapeutics [88], [93], [94], [95]. Similarly, inactivation of the tumor suppressor SMAD4 plays a critical role in the development of pancreatic cancer and correlates with increased tumor aggressiveness and poor prognosis [96], [97], [98], [99]. It is important to note that a common difficulty encountered in the analysis of pancreatic adenocarcinoma tissues is frequent contamination by a dense, desmoplastic stroma that characteristically surround these tumor cells and can occupy up to 90% of a tumor sample’s content [100]. However, in our study of 76 pancreatic cancer samples, we were nonetheless able to extract an 8-transcript signature to distinguish pancreatic samples from other cancer types with high sensitivity and specificity.

As compared to other studies focused on distinguishing tissue of origin for CUP, our study has multiple strengths. We analyzed a large number of cancer samples from 26 different tissue types (3,244 samples as compared to the previous studies analyzing fewer than 800 samples) [15], [16], [17], [18], [19], [23], [24], [25], [26]. In addition, we used multiple validation methods to strengthen our biomarker signatures, specifically reserving 2,556 samples for external validation testing, to yield an overall accuracy of 90.5%. This is as compared to previously reported classification accuracies of 76% to 89% [15], [16], [17], [19], [23], [24]. Finally, the use of RNA-Seq expression data has a number of potential advantages over microarray techniques, as previously outlined, including wide genome coverage, which allowed us to identify several new biomarkers, such as BEST1 in ovarian cancer and the gene fusion ANKHD1-EIF4EBP3 in pancreatic cancer. To our knowledge, this is the first CUP study using large-scale RNA-Seq data for both training and validation to demonstrate a highly accurate model for cancer prediction.

One of the limitations of our study is that, although RNA-Seq allows for the capability to detect unmapped genes, in this proof-of-concept study, we limited our analysis to only annotated transcripts. While computationally more intensive, a dataset comprised of chromosomal positions rather than annotated genes would allow for additional discovery of novel biomarkers and could potentially improve the accuracy. In addition, in our dataset, there were insufficient kidney and liver samples to allocate to our test set for external validation; however, the kidney and liver biomarker signatures nonetheless yielded strong specificities of > 99%.

We have demonstrated the strength of this model in its ability to accurately and efficiently distinguish samples of one type (i.e., cancer type of interest) from another (i.e., heterogeneous group of other cancer types). While this study focused specifically on deriving lineage-specific cancer signatures by RNA-Seq, this model may be applied to any large dataset to query other clinical questions, such as deriving a signature that distinguishes premalignant and cancerous lesions from inflammatory and other benign conditions. Similarly, given a dataset of patients with known clinical outcomes, our model could be used to derive biomarker signatures that identify patients who would respond to a given therapy or identify patients with worse outcomes compared to clinically and histologically-matched cohorts who should be targeted for aggressive treatment and surveillance.

Conclusions

In this study, we introduced a computational model that successfully extracted accurate, lineage-specific cancer signatures for the top eight tissue types that contribute to CUP using RNA-Seq. Through external validation of a large dataset, we have shown how these signatures may be used to accurately identify tumors of unknown origin, demonstrating the translational potential of not only our cancer biomarker signatures but also the model itself, which may be applied to other clinical queries.