A computational framework to trace tumor tissue-of-origin of 19 cancer types based on RNA sequencing

doi:10.21203/rs.3.rs-1457167/v1

PDF

Research Article

A computational framework to trace tumor tissue-of-origin of 19 cancer types based on RNA sequencing

https://doi.org/10.21203/rs.3.rs-1457167/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Carcinoma of unknown primary (CUP) is a type of metastatic cancer with tissue-of-origin (TOO) unidentifiable by traditional methods. Most CUP patients have poor prognosis since no therapy targeting TOO is allowed. Thus, it’s critical to develop accurate computational methods to infer TOO. While qPCR or microarray-based methods are effective in predicting TOO for most cancer types, the overall prediction accuracy is yet to be improved. Here, we propose a computational framework to trace TOO of 19 cancer types based on RNA sequencing (RNA-seq). Specifically, we download the RNA-seq data of 7000+ tissue samples covering 19 cancer types with known TOO from TCGA. By feature selection, 90 genes are finally selected to train a random forest model for TOO inference; the 90 genes are enriched in both tissue-specific functions and tissue-general functions. The cross-validation accuracy of our framework reaches 97.55% across all cancer types. Furthermore, we collected an independent cohort of samples in GEO as testing samples. The accuracy on the independent data is 74% despite the differences in experiment procedures and pipelines. In conclusion, we develop an accurate yet robust computational framework for identifying TOO, which might be promising in clinical applications.

Cancers of unknown primary

tissue-of-origin

RNA sequencing

random forest

feature selection

Carcinoma of unknown primary (CUP) is a type of metastatic cancers with unknown cancer origin. CUP occupies 3–5% of all cancer incidences in the United States [1, 2]. Although there is no drug specifically approved for CUP, multiple guidelines have been approved for treating this disease using multiagent cytotoxic chemotherapy [3, 4]. However, the responses of CUP patients to non-targeted chemotherapies are poor with a 5-year survival rate around 11% [5, 6].

In order to solve this problem, there are multiple diagnostic methods proposed in the past decades. From 1980s to 2010s, immunohistochemistry is the mainstream method to identify cancer primary tissue [7–14]. However, this method is labor-intensive, requires highly skilled physicians, and has different predicted accuracies across different cancer types [15]. Imaging techniques such as PET/CT and ultrasound can also assist in clinical diagnosis on CUP [16–19]. However, they require experienced pathologists and their diagnostic accuracies vary from ~ 30% to ~ 90%, which is not high enough for safe clinical usage.

With the development of sequencing techniques, gene expression profiling-based methods are promising in auxiliary diagnosis for CUP patients [20–28]. The theory behind these methods is that the expression profile of metastatic tumor tissue is similar to that of the original tumor site instead of the metastatic tissue [29].

For example, Bloom et al. combined the cDNA and oligonucleotide platform with artificial neural network (ANN) to trace the primary tumor origin [27], and obtained an accuracy of 83–88% on different platforms, showing that the molecular profiling can be used to trace the primary origin of a metastatic tumor. Tissue of Origin (TOO) is a product based on the microarray data [30]. In the study, 100 samples with primary tumors were used as the initial training set. The Wilcoxon rank score was used to obtain the highly distinguishable genes for every cancer and merged to 1100 genes. One-versus-rest Support vector machine (ovr-SVM) was used to estimate the predictive accuracy of the genes. In the independent validation set, an accuracy of 85% was gained. In addition, Monzon et al. examined 2039 tumors using a dataset containing 1550-gene expression profile [26]. The model was trained on both primary and metastatic tumors as well as the well-differentiated and undifferentiated tumors. The final accuracy was 87.8%, which was similar to the results of Su’s study.

CancerTYPE ID is a product based on RT-PCR data, in which 578 labeled samples covering 39 tumor types were included in the dataset for the tissue-of-origin tracing. Among them, 75% were primary tumors and 25% were metastatic tumors. The dataset was split to a 466-sample dataset (frozen) and a 112-sample test set (FFPE) according to the sample type. K-nearest neighbor algorithm (KNN) (k = 5) was applied to the classification problem, which reached an accuracy of 84% in the leave-one-out cross validation. The result also showed there was no difference in the accuracies on predictions of primary or metastatic tumor. 30 genes per tumor were chosen, resulting in 979 distinct genes. A taxonomy tree was constructed and 30 genes at each tree node were selected, merging to 1001 distinct genes. 61–80 genes were chosen from the 1001 genes by genetic algorithm and KNN. The accuracy reached 86% by using the 74-gene list on the test set. And among the 100 genetic algorithm (GA) selected gene sets, 90 genes showed up in 20 sets. The union of the 74-gene set and the 90-gene set generated a 126-gene set. Furthermore, the criterion of PCR assay filtered another 39 genes, leaving 87 genes at last. By appending 5 reference genes, a 92-gene list was generated [25].

Xu et al. reported a multiple-platform 154-gene panel to detect the primary origin of metastatic tumors and the overall accuracy was 92% [31]. They chose the 154 genes by recursive feature selection followed by specific filtering standards. The array data was used in the training phase and the RNA-seq data was used in the validation phase.

Recently, Liang et al. contributed to the field of CUP by using copy number alteration (CNA) as the input and deep learning as the method for accurate CUP inferring for six cancer types [32]. However, since the CNA data is less available than other kinds of data such as expression profiles or mutation profiles, inferring CUP using CNA is not recommended currently.

Despite the relatively good performances shown from the above-mentioned attempts, there are a few limitations. First, microarray, RT-PCR and RNA-seq are quite different platforms. The model training in one platform might not work well on another platform. The RNA-seq data has been proven more reliable than microarray, and thus a model depending purely on RNA-seq might be preferable. Second, ensemble-based methods like random forest usually have better performance than those of genetic algorithm and KNN. Here, we proposed a computational framework to trace tumor tissue-of-origin of 19 cancer types based on RNA sequencing using the random forest algorithm. Specifically, we applied random forest in feature selection and model training on the cancer genome atlas (TCGA) data. We then used an independent RNA-seq dataset from Gene Expression Omnibus (GEO) to validate our model.

Dataset preparation

We used 2 sources of data in this work. The largest one is 8085 samples from TCGA were collected and further split into training dataset and test dataset. For independent validation, another dataset from GEO was obtained. The GEO datasets containing metastatic 42 samples were collected by searching the keywords for each cancer on the GEO website (see Table S1 for more detail). Samples were prepared as described in Materials and Methods and RNA-seq data were prepared for the predictions by the trained random forest model. All datasets described here were summarized in Fig. 1.

First, 8085 RNA-seq profiles of samples in TCGA [33] covering 21 main cancers (all cancer abbreviations are supplied in Table 1) were collected. The TCGA data were quantified by RSEM method. In TCGA, lung cancer was further classified as lung squamous cell carcinoma and lung adenocarcinoma since they are classified as two major tumors in TCGA. We also merged the two cancers, colon adenocarcinoma (COAD) and rectum adenocarcinoma (READ) to COADREAD, since they are similar to each other [34]. T-SNE plot for COAD and READ also showed that they can hardly be distinguished (Fig. 2a). Table 1 showed the distribution of the collected data of all cancers. The TCGA dataset contains 7713 primary tumor samples and 372 metastatic tumor samples. We checked the distribution of each cancer (shown in Table 1), and found that 352 (94.62%) metastatic tumor samples belong to Skin Cutaneous Melanoma (SKCM), which is much higher than that in the primary tumor samples (1.04%). Due to the distinct distribution of the primary and metastatic tumor, the SKCM samples were removed from the data in case it would affect the results, leaving 7633 samples covering 19 cancers in the primary tumor sample dataset and 20 metastatic tumor samples covering 5 cancers (BRCA, CESC, COAD, HNSC and THCA). The 20 metastatic tumor samples were used as an independent validation dataset.

Table 1

The number of patients per cancer from TCGA database.
Cancer	Amount	Percentage	Full name
BLCA	301	3.90%	Bladder Urothelial Carcinoma
BRCA	1056	13.69%	Breast invasive carcinoma
CESC	258	3.35%	Cervical squamous cell carcinoma and endocervical adenocarcinoma
COAD	451	5.85%	Colon adenocarcinoma
GBM	153	1.98%	Glioblastoma multiforme
HNSC	480	6.22%	Head and Neck squamous cell carcinoma
KIRC	526	6.82%	Kidney renal clear cell carcinoma
KIRP	222	2.88%	Kidney renal papillary cell carcinoma
LAML	173	2.24%	Acute Myeloid Leukemia
LGG	439	5.69%	Brain Lower Grade Glioma
LIHC	294	3.81%	Liver hepatocellular carcinoma
LUAD	486	6.30%	Lung adenocarcinoma
LUSC	428	5.55%	Lung squamous cell carcinoma
OV	261	3.38%	Ovarian serous cystadenocarcinoma
PAAD	142	1.84%	Pancreatic adenocarcinoma
PRAD	379	4.91%	Prostate adenocarcinoma
READ	153	1.98%	Rectum adenocarcinoma
SKCM	80	1.04%	Skin Cutaneous Melanoma
STAD	415	5.38%	Stomach adenocarcinoma
THCA	500	6.48%	Thyroid carcinoma
UCEC	516	6.69%	Uterine Corpus Endometrial Carcinoma

Among the 7633 primary tumor samples, 7579 were frozen samples and 49 were Formalin-Fixed Paraffin Paraffin-embedded (FFPE) samples and 5 samples were ambiguous since they were labeled different types in their sample or portion information. Since FFPE samples sustained from RNA degradation and chemical modifications [35, 36], the sample storage method was also taken as a factor that might influence the classification accuracy.

Since there might be outliers in the data, such as wrong labels or technical noise from the surgery or the experiment, we performed outlier removal using one-class support vector machine [37]. The outlier detection was performed as described in Materials and Methods. After outlier removal, 7174 samples were left, and the sample number per cancer were listed in Table 2.

Table 2

Sample numbers from the TCGA frozen dataset for each cancer after removing outliers by SVM.
Cancer	Original number	Number after removal
BLCA	298	283
BRCA	1041	985
CESC	257	244
COADREAD	591	562
GBM	153	145
HNSC	480	455
KIRC	522	493
KIRP	222	211
LAML	173	162
LGG	439	417
LIHC	294	278
LUAD	476	448
LUSC	428	404
OV	261	247
PAAD	142	130
PRAD	375	356
STAD	415	393
THCA	500	477
UCEC	512	484

Random forest algorithm performed well in frozen datasets during cross validation

We ran 10-fold cross validations on the frozen datasets to evaluate the algorithm and feature number to use. For most cases, it’s relatively less costly to use panels to determine the expression level of the specific genes. However, the coverage of the panel is also confined to the cost of the panel. Our goal is to lower the cost (gene number) while keeping the primary tracing accuracy as high as possible. We used the of random forest algorithm to choose the most discriminating genes (see Materials and Methods). We chose the best G genes from the training set in each fold of cross validation, where G is an integer between 10 and 200 with 10 as the step size.

The average accuracies of all 10-fold cross validations were plotted against the feature numbers (Fig. 2b). In the frozen dataset, the accuracy stabilized after 90 genes (accuracy = 97.55%). When using more genes, the accuracy increases for less than 0.1% for each increase of gene number. Thus, 90 were set as the optimal feature number for frozen datasets. The confusion matrix of the cross validation of the 90 genes from the frozen dataset was plot as Fig. 2c.

Selected genes were significantly enriched in both tissue-specific functions and tissue-general functions.

We selected 90 genes for further usage in frozen dataset. The log2 transformed average expression values for selected genes in different cancers were plotted as heatmaps (Fig. 3). In some extent, the heatmap explained how those genes could be used to distinguish different cancer types. Some of the genes can be very informative to make the judgment of the primary origin, as they have a much different expression in some cancers. For example, SCGB2A2 has an extremely low expression in LAML compared to other cancers while it has a very high expression in BRCA, showing the gene expression levels could be used to distinguish different cancer types.

We also performed enrichment analysis on the gene set. The enrichment results were shown in Fig. 4. The genes were enriched significantly in development of human organs, such as urogenital system development, endocrine system development, thyroid gland development, prostate gland morphogenesis, renal system development, lung cell differentiation and reproductive system development et al. There are other genes involving in hormone metabolic process, hormone-mediated signaling pathway. Remarkably, some genes involved in the organ or tissue-specific development are more likely to be differentially expressed in tumors and normal tissues. For example, HOXB13 gene, which encodes nuclear transcription factor which belongs to the HOX family, is known to involve in embryonic development and in the termination of the differentiated tissue [38, 39]. The previous studies revealed that HOXB13 was overexpressed in numerous cancers, such as prostate cancer [40], breast cancer [41], ovarian cancer [42], cervical cancer [43], which demonstrated that HOXB13 might be a promising biomarker of carcinogenesis in certain types of cancer. In addition, SOX17 gene encodes a member of the SOX family (SRY-related HMG-box) of transcription factors involved in multiple developmental processes [44], as well as in carcinogenesis [45]. It has been reported that the SOX17 expression was down-regulated in several types of cancer cells [46–48], and it was associated with a poor prognosis of cancers. [48–50]. Remarkably, we found ESR1 was a gene signature of our 90 genes panel. ESR1 encodes a nuclear estrogen receptor which is involved in the regulation of cell proliferation and differentiation in target tissues. The overexpressed ESR1 gene results in two-thirds of breast cancer and it has been defined the ER + subtype of breast cancer [51]. ESR1 is known as a biomarker to evaluate endocrine resistance, disease relapse and an overall poor prognosis [52–54].

The model trained from TCGA frozen dataset performed well on FFPE dataset in specific cancer types.

We first checked whether the 90-gene set could help to discriminate the samples from different cancers. t-SNE were used to visualize the samples from 19 cancer types of the frozen dataset (Fig. 5a) and the samples the FFPE dataset (Fig. 5b) using only 90 genes. The samples from different cancer types were clustered well using much fewer genes. However, the FFPE dataset was less clustered, indicating the FFPE samples might not be a good choice for predicting primary origin.

We trained a model using 90 genes on the 7194-sample frozen dataset. We then predicted the 49-sample FFPE dataset using the trained model, gaining the overall accuracy of 59.18%. We checked the precisions for each tumor type in the predictions of FFPE samples, and found that the STAD and COADREAD were not well classified since COADREAD samples were prone to be predicted as STAD, leading to a low precision on STAD and a low recall on COADREAD. This indicated that the two tumor types tended to be similar after being fixed by formalin. Thus, we need to take the sample storage method into account in further investigations. There are 6 samples having a low confidence in predicting, thus we set a threshold of 0.3 for the max probability of all cancers, leaving 43 predictable samples, which lifted the accuracy to 67.44%. We plotted the confusion matrix for the FFPE dataset in Fig. 5c.

The model trained from TCGA dataset performed well in independent datasets.

For testing our method, we prepared 2 independent datasets: (1) the metastatic dataset from TCGA containing 20 metastatic samples; (2) a metastatic dataset from GEO containing 42 samples.

The model trained from TCGA primary tumor dataset achieves 94.74% accuracy in metastatic tumor dataset for 19 predictable samples. The model was tested on the independent metastatic dataset in TCGA, outputting the result shown in Table S2 and Table S3.

To expand the validation dataset, we further validated the trained models on the GEO datasets. The STAR-RSEM pipeline was used to extract the expression value for all genes [55]. The GEO metastatic dataset contained 42 predictable samples and an accuracy of 74.29% was achieved. We plotted the confusion matrix for the FFPE dataset in Fig. 5d. Specifically, GSE83132 contains samples from cell lines, which might lead to inaccurate RNA expression profiling.

Here we demonstrated the random forest algorithm that could be used in the problem of primary-tracing using the RNA-seq data with minimal gene set. The overall accuracy using less than 100 genes is high in the cross-validation results, showing its fitting ability and robustness. The gene number is positively related to the overall accuracy. However, the best cross validation result came from the 90-gene or 60-gene set instead of the 100-gene set, showing that there might be a risk of over-fitting if too many genes were used as features. Thus, some further methods performed well in other fields should be explored [56, 57].

Figure 2c showed the confusion matrix of the classification using the 90 genes, in which the most samples lay on the diagonal of the matrix. It’s worthy to point out that the recall of the READ was relatively low and most of the erroneous READ samples were classified as the COAD, leading to a relatively low precision for COAD. It might be due to the connectivity between colon and rectum, as well as the similarity in expression profile. It is necessary to further study to accurately classify the two relatively adjacent cancers.

We also showed that the FFPE samples might affect the result of primary tracing prediction, which was not mentioned by Ma et al. who also used mixed datasets [25]. If possible, frozen and FFPE datasets should not be simply mixed for a larger dataset. Instead, STAD and COADREAD samples should be dealt with carefully.

It’s also worthy to point out that the probability filtration could increase the accuracy tremendously. Since the filtration is a procedure that drops some low-confidence samples, there are also some results which have correct labels with relatively low probability, the procedure might lower the recall rate. Therefore, it requires careful design to set the threshold at which the max probability and gap. Besides, further optimization should be made to raise the model’s generalization performance.

We found that during the test phase using the mRNA profiles of the samples from the GEO datasets, the accuracy remained high, showing the robustness of the random forest model and the selected features. However, the accuracy dropped a little bit in the validation datasets, which is possibly due to differences in experiment protocols that might also affect the final expression quantification, which might furthermore affect the prediction result. Further study could be focused on the effect of those differences.

Data preparation

The TCGA RNA-seq and array data were downloaded from ICGC Data Portal (https://dcc.icgc.org/releases/release_26/).

The normalized expression value of each sample and each gene from the downloaded table of each cancer type were extracted, generating a M * N matrix where M is the number of the sample number and N is the number of the gene number. All the samples were labeled by its cancer type.

Since the primary goal is to trace the primary tumor, we removed all the metastatic tumors in the TCGA dataset for test set by checking TCGA identifier, leaving the samples with primary tumors as training dataset.

For the GEO dataset, we further filtered samples by the following criteria: samples containing less than 12853 genes were removed for those samples might be problematic in sequencing. The number 12853 was inferred from TCGA dataset, since the sample TCGA-AA-A004-01A-01R-A00A-07 has the least non-zero genes, which contains only 12853 non-zero genes.

Outlier removal

The scikit-learn class OneClassSVM uses the support vector machine to fit the bound of the samples. It was used to detect the outliers and the nu parameter was set to 0.05 and the rbf kernel was used. All outliers will be removed before any training.

Feature selection

Considering the gene number for all genes is costly, the gene number could be shrunk to a small size to lower the cost.

The random forest algorithm [58] was applied for the feature selection. Random forest consists of a number of binary decision trees. Each decision tree could use up to $\sqrt{selected gene number}$ genes for classification. Every node in the decision tree is a judge on a single gene, which is designed to split the dataset into two so that similar expression values end up in the same set. Gini impurity was used as the criterion to further split each node. The Gini impurity was calculated by

$${I}_{g}\left(f\right)=1-\sum _{i=1}^{n}{{f}_{i}}^{2}$$

where $f$ is a vector of all frequencies of different cancers in the current node and ${f}_{i}$ is the frequency of the $i$th cancer.

For a single tree, the decrease of the Gini impurity splitting at each feature could be calculated. For a forest, the impurity decrease from each feature can be averaged and the features are ranked according to this measure. For every node in every decision tree, the Gini impurity was calculated. By traversing all possible values of all left features, the best feature-value combination will be decided according to the most dramatic drop of Gini impurity. The decrease of the Gini impurity, i.e., the importance score for the feature will be calculated by

$${\varDelta I}_{g}\left(F\right)={I}_{g}-{I}_{g0}-{I}_{g1}$$

where $F$ is a feature, ${I}_{g}$ is the Gini impurity before splitting, ${I}_{g0}$ and ${I}_{g1}$ is the Gini impurity of the left/right node after splitting. The value chosen for the feature is hidden in convenience.

After calculating the decrease for all features in all trees, the importance score for a feature in the random forest will be calculated by

$$S\left(F\right)=\frac{\sum _{t}^{T}\sum _{F}^{A}{\varDelta I}_{g}\left(F\right)}{\sum _{t}^{T}I\left(F in t\right)}$$

Where $F$ is a specific feature, $T$ is the set of all trees, $t$ is one decision tree from the set $T$, $A$ is all features available for tree $t$, and $I\left(F in t\right)$ is the indicator function for whether feature $F$ is available for tree $t$.

After sorting by importance scores of all features, the top N genes with the highest scores are selected accordingly to the requirement.

The selected genes within all the samples were extracted to the new matrix for further classification procedure.

In this work, 2000 base estimators were used to choose the best features.

Functional annotation

For the analysis of biological significance, the functions were annotated for the specific gene set. Gene ontology [59, 60] was used as the database for the enrichment analysis. Genes were clustered by R package clusterProfiler [61]. The visualization was done by R package ggplot2 [62].

Cross-validation

In every N-fold cross validation (where N is an integer), all the samples were stratified into N subsets by different random seeds. And the algorithm was repeated N times. During each repeat, one of the N subsets was used as the test set and the other N-1 subsets were consolidated to a training set. Features were selected within the training set and those features were from the training set was used to train a model. The test set was then used to evaluate the model. Then the average error across all N trials was computed.

Classification

The samples were labeled as their cancer type. After feature selection, the expression profile corresponding to the genes in the gene set of each sample was extracted and l1-normalized within each sample. The l1-normalized gene expression levels were calculated by

$${G}_{i}=\frac{{g}_{i}}{\sum _{i=1}^{n}{g}_{i}}$$

where ${G}_{i}$ is the $i$th normalized gene expression level for one sample, ${g}_{i}$ is corresponding to the raw gene expression level and $n$ is the total number of the selected genes.

Random forest was used for the classification of different cancers. The base estimator number was set to 2000 for the random forest classifier. For each decision tree, sub-samples are drawn with replacement by bootstrapping method. Each decision tree will use up to $\sqrt{selected gene number}$ genes to be compatible with the feature selection process. For each tree in the random forest, the Gini impurity was calculated for each possible feature before and after the next split and the feature that leads to the greatest decline will be used. The nodes will be split until all samples in the node were of the same cancer type.

The trained decision trees were used as voters for the cancer classification and the cancer with the largest number of trees was reported. The probability corresponding to each cancer were also given by the random forest classifier. The results will be further filtered if necessary.

The feature selection and classification were implemented using scikit-learn package [37].

Probability filtration

The largest probability that the classifier predicts one sample could be reckoned as the confidence of the classifier. However, even the largest probability could be small, indicating a low-confidence prediction. Thus, the results could be further filtered if necessary.

There were 2 criteria for the filtration considering both the value of the probability and the gap between the largest probability and the second largest probability.

The first criterion is applied by

$$\text{max} p>M$$

where $p$ is the probability vector containing the probabilities corresponding to all the cancers and $M$ is the minimum threshold for the $\text{max}p$. The second criterion is applied by

$${p\text{'}}_{1}-{p\text{'}}_{2}>G$$

where $p\text{'}$ is the sorted $p$ in descending order and $G$ is the minimum threshold for the gap between the largest probability and the second largest probability. Those predictions that fell out of either of the criteria will be dropped.

RNA sequencing data alignment and quantification

STAR (v2.7.2b) is used for the alignment of the raw sequencing data [63]. The hg19 genome was used as the reference genome, using default parameters except the outFilterType was set to BySJout.

RSEM (v1.3.1) is used for quantification of all genes after the alignment, which is based on expectation maximization algorithm [55]. The default parameters are used in this procedure.

Heatmap visualization

To plot the heatmaps, the R package pheatmap was used. Before plotting, the expression values were first added 1e-12 and transformed by log2.

t-SNE visualization

Since expression profile was in very high dimension. Thus, the visualization of such data requires dimension reduction. t-SNE is a popular way to reduce the dimension for many kinds of data. t-SNE could convert high-dimensional data into low-dimensional data but still keeping their proximity [64]. manifold.TSNE from Scikit-learn package [37] was used to reduce the dimension of the data and matplotlib package was used for visualization.

In conclusion, this study developed a method by RNA-seq and machine learning to assist in identification and validation of the primary origin of tumor tissue. Our study focused on 19 types tumors and applied random forest algorithm using only 90 genes as features, which resulted in a relatively high accuracy in terms of predicting the origin of the primary cancer site. This work indicated the RNA-seq and machine learning algorithms might be used in clinical practice when pathological study alone fails to indicate the primary origin site of certain cancer. In order to ensure the precision of the classification, the probability could also be a constraint, though there was a tradeoff between the precision and recall.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

All classification codes are available at https://github.com/wangbo00129/ClassifyBySklearn. Data used in this work were in the Result section Dataset preparation.

Competing interests

The authors declare no conflict of interest.

Funding

This study was partially supported by the Natural Science Foundation of Hunan, China (Grant No. 2018JJ2461).

Authors' contributions

JY and GT conceived the project; BW and HY implemented the experiments and analyzed the data; YZ prepared the data and performed literature search; BW, YZ and JY wrote the manuscript; all authors approved the final manuscript.

Acknowledgments

Not applicable.

Sokilde R, Vincent M, Moller AK, Hansen A, Hoiby PE, Blondal T, Nielsen BS, Daugaard G, Moller S, Litman T: Efficient identification of miRNAs for classification of tumor origin. The Journal of molecular diagnostics: JMD 2014, 16(1):106–115.
He B, Lang J, Wang B, Liu X, Lu Q, He J, Gao W, Bing P, Tian G, Yang J: TOOme: A Novel Computational Framework to Infer Cancer Tissue-of-Origin by Integrating Both Gene Mutation and Expression. Frontiers in bioengineering and biotechnology 2020, 8:394.
Natoli C, Ramazzotti V, Nappi O, Giacomini P, Palmeri S, Salvatore M, Landriscina M, Zilli M, Natali PG, Tinari N et al: Unknown primary tumors. Biochimica et biophysica acta 2011, 1816(1):13–24.
Agwa E, Ma PC: Overview of various techniques/platforms with critical evaluation of each. Current treatment options in oncology 2013, 14(4):623–633.
Varadhachary GR, Raber MN: Carcinoma of unknown primary site. N Engl J Med 2014, 371(21):2040–2040.
He B, Dai C, Lang J, Bing P, Tian G, Wang B, Yang J: A machine learning framework to trace tumor tissue-of-origin of 13 types of cancer based on DNA somatic mutation. Biochim Biophys Acta Mol Basis Dis 2020, 1866(11):165916.
Yam LT, Janckila AJ, Lam WK, Li CY: Immunohistochemistry of prostatic acid phosphatase. The Prostate 1981, 2(1):97–107.
de Almeida PC, Pestana CB: [Use of immunohistochemistry in detecting the primary site in neoplasm metastasis]. AMB: revista da Associacao Medica Brasileira 1989, 35(3):84–87.
de Almeida PC, Pestana CB: Immunohistochemical markers in the identification of metastatic breast cancer. Breast cancer research and treatment 1992, 21(3):201–210.
Brown RW, Campagna LB, Dunn JK, Cagle PT: Immunohistochemical identification of tumor markers in metastatic adenocarcinoma. A diagnostic adjunct in the determination of primary site. American journal of clinical pathology 1997, 107(1):12–19.
Nap M: Immunohistochemistry of CA 125. Unusual expression in normal tissues, distribution in the human fetus and questions around its application in diagnostic pathology. The International journal of biological markers 1998, 13(4):210–215.
Hameed O, Humphrey PA: Immunohistochemistry in diagnostic surgical pathology of the prostate. Seminars in diagnostic pathology 2005, 22(1):88–104.
Park SY, Kim BH, Kim JH, Lee S, Kang GH: Panels of immunohistochemical markers help determine primary sites of metastatic adenocarcinoma. Archives of pathology & laboratory medicine 2007, 131(10):1561–1567.
Idikio HA: Immunohistochemistry in diagnostic surgical pathology: contributions of protein life-cycle, use of evidence-based methods and data normalization on interpretation of immunohistochemical stains. International journal of clinical and experimental pathology 2009, 3(2):169–176.
Kulkarni A, Pillai R, Ezekiel AM, Henner WD, Handorf CR: Comparison of histopathology to gene expression profiling for the diagnosis of metastatic cancer. Diagnostic pathology 2012, 7:110.
Chiti A, Fanti S, Savelli G, Romeo A, Bellanova B, Rodari M, van Graafeiland BJ, Monetti N, Bombardieri E: Comparison of somatostatin receptor imaging, computed tomography and ultrasound in the clinical management of neuroendocrine gastro-entero-pancreatic tumours. European journal of nuclear medicine 1998, 25(10):1396–1403.
Guntinas-Lichius O, Peter Klussmann J, Dinh S, Dinh M, Schmidt M, Semrau R, Mueller RP: Diagnostic work-up and outcome of cervical metastases from an unknown primary. Acta oto-laryngologica 2006, 126(5):536–544.
Kroiss A, Putzer D, Decristoforo C, Uprimny C, Warwitz B, Nilica B, Gabriel M, Kendler D, Waitz D, Widmann G et al: 68Ga-DOTA-TOC uptake in neuroendocrine tumour and healthy tissue: differentiation of physiological uptake and pathological processes in PET/CT. European journal of nuclear medicine and molecular imaging 2013, 40(4):514–523.
Prowse SJ, Shaw R, Ganeshan D, Prowse PM, Hanlon R, Lewis-Jones H, Wieshmann H: The added value of 18F-fluorodeoxyglucose positron emission tomography computed tomography in patients with neck lymph node metastases from an unknown primary malignancy. The Journal of laryngology and otology 2013, 127(8):780–787.
Golub TR, Slonim DK, Tamayo P, Huard C, Gaasenbeek M, Mesirov JP, Coller H, Loh ML, Downing JR, Caligiuri MA et al: Molecular classification of cancer: class discovery and class prediction by gene expression monitoring. Science 1999, 286(5439):531–537.
Ramaswamy S, Tamayo P, Rifkin R, Mukherjee S, Yeang CH, Angelo M, Ladd C, Reich M, Latulippe E, Mesirov JP et al: Multiclass cancer diagnosis using tumor gene expression signatures. Proceedings of the National Academy of Sciences of the United States of America 2001, 98(26):15149–15154.
Greco FA, Erlander MG: Molecular classification of cancers of unknown primary site. Molecular diagnosis & therapy 2009, 13(6):367–373.
Monzon FA, Koen TJ: Diagnosis of metastatic neoplasms: molecular approaches for identification of tissue of origin. Archives of pathology & laboratory medicine 2010, 134(2):216–224.
Rosenwald S, Gilad S, Benjamin S, Lebanony D, Dromi N, Faerman A, Benjamin H, Tamir R, Ezagouri M, Goren E et al: Validation of a microRNA-based qRT-PCR test for accurate identification of tumor tissue origin. Modern pathology: an official journal of the United States and Canadian Academy of Pathology, Inc 2010, 23(6):814–823.
Ma XJ, Patel R Fau - Wang X, Wang X Fau - Salunga R, Salunga R Fau - Murage J, Murage J Fau - Desai R, Desai R Fau - Tuggle JT, Tuggle Jt Fau - Wang W, Wang W Fau - Chu S, Chu S Fau - Stecker K, Stecker K Fau - Raja R et al: Molecular classification of human cancers using a 92-gene real-time quantitative polymerase chain reaction assay. 2005(1543–2165 (Electronic)).
Monzon FA, Lyons-Weiler M, Buturovic LJ, Rigl CT, Henner WD, Sciulli C, Dumur CI, Medeiros F, Anderson GG: Multicenter validation of a 1,550-gene expression profile for identification of tumor tissue of origin. Journal of clinical oncology: official journal of the American Society of Clinical Oncology 2009, 27(15):2503–2508.
Bloom G, Yang IV, Boulware D, Kwong KY, Coppola D, Eschrich S, Quackenbush J, Yeatman TJ: Multi-platform, multi-site, microarray-based human tumor classification. The American journal of pathology 2004, 164(1):9–16.
Liu H, Qiu C, Wang B, Bing P, Tian G, Zhang X, Ma J, He B, Yang J: Evaluating DNA Methylation, Gene Expression, Somatic Mutation, and Their Combinations in Inferring Tumor Tissue-of-Origin. Front Cell Dev Biol 2021, 9:619330.
Erlander MG, Ma XJ, Kesty NC, Bao L, Salunga R, Schnabel CA: Performance and clinical evaluation of the 92-gene real-time PCR assay for tumor classification. The Journal of molecular diagnostics: JMD 2011, 13(5):493–503.
Su AI, Welsh Jb Fau - Sapinoso LM, Sapinoso Lm Fau - Kern SG, Kern Sg Fau - Dimitrov P, Dimitrov P Fau - Lapp H, Lapp H Fau - Schultz PG, Schultz Pg Fau - Powell SM, Powell Sm Fau - Moskaluk CA, Moskaluk Ca Fau - Frierson HF, Jr., Frierson Hf Jr Fau - Hampton GM et al: Molecular classification of human carcinomas by use of gene expression signatures. 2001(0008-5472 (Print)).
Xu Q, Chen J, Ni S, Tan C, Xu M, Dong L, Yuan L, Wang Q, Du X: Pan-cancer transcriptome analysis reveals a gene expression signature for the identification of tumor tissue origin. Modern pathology: an official journal of the United States and Canadian Academy of Pathology, Inc 2016, 29(6):546–556.
Liang Y, Wang H, Yang J, Li X, Dai C, Shao P, Tian G, Wang B, Wang Y: A Deep Learning Framework to Predict Tumor Tissue-of-Origin Based on Copy Number Alteration. Frontiers in bioengineering and biotechnology 2020, 8:701.
Cancer Genome Atlas Research N: Comprehensive genomic characterization defines human glioblastoma genes and core pathways. Nature 2008, 455(7216):1061–1068.
Zuo S, Dai G, Ren X: Identification of a 6-gene signature predicting prognosis for colorectal cancer. Cancer cell international 2019, 19:6.
Pennock ND, Jindal S, Horton W, Sun D, Narasimhan J, Carbone L, Fei SS, Searles R, Harrington CA, Burchard J et al: RNA-seq from archival FFPE breast cancer samples: molecular pathway fidelity and novel discovery. (1755–8794 (Electronic)).
Landolt L, Marti HP, Beisland C, Flatberg A, Eikrem OS: RNA extraction for RNA sequencing of archival renal tissues. (1502–7686 (Electronic)).
Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, Blondel M, Prettenhofer P, Weiss R, Dubourg V et al: Scikit-learn: Machine Learning in Python. J Mach Learn Res 2011, 12(null):2825–2830.
Decker B, Ostrander EA: Dysregulation of the homeobox transcription factor gene HOXB13: role in prostate cancer. Pharmacogenomics and personalized medicine 2014, 7:193–201.
Chung N, Jee BK, Chae SW, Jeon YW, Lee KH, Rha HK: HOX gene analysis of endothelial cell differentiation in human bone marrow-derived mesenchymal stem cells. Molecular biology reports 2009, 36(2):227–235.
Ewing CM, Ray AM, Lange EM, Zuhlke KA, Robbins CM, Tembe WD, Wiley KE, Isaacs SD, Johng D, Wang Y et al: Germline mutations in HOXB13 and prostate-cancer risk. N Engl J Med 2012, 366(2):141–149.
Ang MK, Ooi AS, Thike AA, Tan P, Zhang Z, Dykema K, Furge K, Teh BT, Tan PH: Molecular classification of breast phyllodes tumors: validation of the histologic grading scheme and insights into malignant progression. Breast cancer research and treatment 2011, 129(2):319–329.
Miao J, Wang Z, Provencher H, Muir B, Dahiya S, Carney E, Leong CO, Sgroi DC, Orsulic S: HOXB13 promotes ovarian cancer progression. Proceedings of the National Academy of Sciences of the United States of America 2007, 104(43):17093–17098.
Lopez R, Garrido E, Pina P, Hidalgo A, Lazos M, Ochoa R, Salcedo M: HOXB homeobox gene expression in cervical carcinoma. International journal of gynecological cancer: official journal of the International Gynecological Cancer Society 2006, 16(1):329–335.
Kamachi Y, Kondoh H: Sox proteins: regulators of cell fate specification and differentiation. Development 2013, 140(20):4129–4144.
Dong C, Wilhelm D, Koopman P: Sox genes and cancer. Cytogenetic and genome research 2004, 105(2–4):442–447.
Yang T, Li XN, Li L, Wu QM, Gao PZ, Wang HL, Zhao W: Sox17 inhibits hepatocellular carcinoma progression by downregulation of KIF14 expression. Tumour biology: the journal of the International Society for Oncodevelopmental Biology and Medicine 2014, 35(11):11199–11207.
Zhang W, Glockner SC, Guo M, Machida EO, Wang DH, Easwaran H, Van Neste L, Herman JG, Schuebel KE, Watkins DN et al: Epigenetic inactivation of the canonical Wnt antagonist SRY-box containing gene 17 in colorectal cancer. Cancer research 2008, 68(8):2764–2772.
Kuo IY, Wu CC, Chang JM, Huang YL, Lin CH, Yan JJ, Sheu BS, Lu PJ, Chang WL, Lai WW et al: Low SOX17 expression is a prognostic factor and drives transcriptional dysregulation and esophageal cancer progression. International journal of cancer 2014, 135(3):563–573.
Yin D, Jia Y, Yu Y, Brock MV, Herman JG, Han C, Su X, Liu Y, Guo M: SOX17 methylation inhibits its antagonism of Wnt signaling pathway in lung cancer. Discovery medicine 2012, 14(74):33–40.
Tang CY, Lin J, Qian W, Yang J, Ma JC, Deng ZQ, Yang L, An C, Wen XM, Zhang YY et al: Low SOX17 expression: prognostic significance in de novo acute myeloid leukemia with normal cytogenetics. Clinical chemistry and laboratory medicine 2014, 52(12):1843–1850.
Bailey SD, Desai K, Kron KJ, Mazrooei P, Sinnott-Armstrong NA, Treloar AE, Dowar M, Thu KL, Cescon DW, Silvester J et al: Noncoding somatic and inherited single-nucleotide variants converge to promote ESR1 expression in breast cancer. Nature genetics 2016, 48(10):1260–1266.
Angus L, Beije N, Jager A, Martens JW, Sleijfer S: ESR1 mutations: Moving towards guiding treatment decision-making in metastatic breast cancer patients. Cancer treatment reviews 2017, 52:33–40.
Musgrove EA, Sutherland RL: Biological determinants of endocrine resistance in breast cancer. Nature reviews Cancer 2009, 9(9):631–643.
Gururaj AE, Rayala SK, Vadlamudi RK, Kumar R: Novel mechanisms of resistance to endocrine therapy: genomic and nongenomic considerations. Clinical cancer research: an official journal of the American Association for Cancer Research 2006, 12(3 Pt 2):1001s-1007s.
Li B, Dewey CN: RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC bioinformatics 2011, 12:323.
Meng Y, Lu C, Jin M, Xu J, Zeng X, Yang J: A weighted bilinear neural collaborative filtering approach for drug repositioning. Brief Bioinform 2022:bbab581.
Yang J, Ju J, Guo L, Ji B, Shi S, Yang Z, Gao S, Yuan X, Tian G, Liang Y et al: Prediction of HER2-positive breast cancer recurrence and metastasis risk from histopathological images and clinical information via multimodal deep learning. Comput Struct Biotechnol J 2022, 20:333–342.
Breiman L: Random Forests. Machine Learning 2001, 45(1):5–32.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT et al: Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nature genetics 2000, 25(1):25–29.
The Gene Ontology C: The Gene Ontology Resource: 20 years and still GOing strong. Nucleic acids research 2019, 47(D1):D330-D338.
Yu G, Wang LG, Han Y, He QY: clusterProfiler: an R package for comparing biological themes among gene clusters. Omics: a journal of integrative biology 2012, 16(5):284–287.
Ginestet C: ggplot2: Elegant Graphics for Data Analysis. Journal of the Royal Statistical Society 2011, 174(1):245–246.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR: STAR: ultrafast universal RNA-seq aligner. Bioinformatics (Oxford, England) 2013, 29(1):15–21.
Hinton GE: Visualizing High-Dimensional Data Using t-SNE. Journal of Machine Learning Research 2008, 9(2):2579–2605.

No competing interests reported.

Supplementary20220221.docx

PDF

Reviewers agreed at journal
16 Apr, 2022
Reviewers invited by journal
06 Apr, 2022
Editor assigned by journal
17 Mar, 2022
Editor invited by journal
17 Mar, 2022
Submission checks completed at journal
17 Mar, 2022
First submitted to journal
16 Mar, 2022

You are reading this latest preprint version

A computational framework to trace tumor tissue-of-origin of 19 cancer types based on RNA sequencing

Status:

Version 1

Abstract

Figures

Introduction

Results

Dataset preparation

Random forest algorithm performed well in frozen datasets during cross validation

Discussion

Materials And Methods

Data preparation

Outlier removal

Feature selection

Functional annotation

Cross-validation

Classification

Probability filtration

RNA sequencing data alignment and quantification

Heatmap visualization

t-SNE visualization

Conclusions

Declarations

References

Competing Interests

Supplementary Files

Status:

Version 1