Cancer gene expression (mRNA, miRNA and lncRNA) and DNA methylation profiles generated by the Cancer Genome Atlas (TCGA) project were downloaded via the cBioPortal for Cancer Genomics (Cerami et al., 2012). For TCGA gene expression and DNA methylation data, only level one data was employed in our analysis. RNA-SeqV2 was used, which takes transcript length into account and is suggested to provide more accurate results. This work only contains the data of solid tumors, and a data quality control were conducted. For each cancer type, the dataset should have a sufficient number of samples (>100). The DNA methylation profiles were measured by the Infinium HumanMethylation450 platform, and we removed those CpG sites with more than 30% missing sample values. The remaining missing values were calculated using the K-nearest neighbor method. In this work, we adopted the same cohort for each dataset, and a total of 6,738 tumor samples for mRNA, miRNA, lncRNA and DNA methylation-based profile were collected spanning 20 cancer types.

In this work, we employed eight machine learning classifiers for tumor tissue origin detection (Rosenfeld et al., 2008; Moran et al., 2016; Hao et al., 2017; Soh et al., 2017; Tang et al., 2017). These methods differ in their underlying methodology, and detailed descriptions of these models appear below. All these models were implemented in Python packages (v3.9).

Random forest (RF) is an ensemble learning algorithm for classification that works based on a multitude of decision trees. Each tree in the forest is built from a sample set drawn from the training set with replacement. Each feature used to split an internal node in the decision trees are picked from a random subset of the entire feature set. We used the soft voting strategy, i.e., the probabilities assigned to each class are calculated by averaging the output of each decision tree. The number of trees are set to 200.

Support vector machine (SVM) classifier is used to find a hyper-plane to separate two classes through maximizing the distance between the hyper-plane and the support vectors, which are defined as the samples closest to the hyper-plane. For those which are not separable, SVM is able to classify them through mapping the points into a higher dimensional space. In this work, the linear kernel was used. For multi-class classification, we implemented the ‘One-vs-Rest’ (OvR) approach.

K-nearest neighbor (KNN) is a non-metric method for classification. KNN simply saves all the samples in training set. Each time a test sample is given, KNN calculates the distances between the sample and all the data points in the training set. The test sample is classified as the class most common among its k nearest neighbors. Here, we used Euclidean distance, and set K to 5.

Decision trees (DT) are tree-like models used for classification. Each internal node within a decision tree represents a classification rule and each leaf node represents a class label. We built Classification and Regression Trees (CART), which choose features through minimizing the Gini index at each node.

Least Absolute Shrinkage and Selector Operator (LASSO) is a linear classification model that uses L1-regularization strategy in parameter estimation. The probabilities of each class are calculated via logistic function. To avoid over fitting, the one norm of the coefficient was added in the loss function, and coefficients were calculated through minimizing the loss function.

Neural network (NN) is based on a collection of connected nodes called neurons. Each neuron receives the input signals of other neurons through weighted connection, and produces output through activation function. If the weighted sum of input signals exceeds a cutoff, the neuron will be activated and outputs a non-zero value. Here we use the rectified linear unit function (ReLU) as the activation function. We used multi-layer feed-forward neural network, where all the neurons are connected with the next layer, and neurons within a layer are not connected with each other. The network was trained using error BackPropagation (BP) algorithm.

Naïve Bayesian classifier (NBC) is a classifier based on Bayes’ theorem and the assumption of independence among all features. Assuming that all feature is independent, the joint distribution equals the multiplication of marginal distribution. The probabilities assigned to each class is calculated through Bayes’ theorem, and the predicted class is the class with the largest probability value.

Linear discriminant analysis (LDA) is a linear model used in classification. Given a dataset with two classes, LDA projects all the sample points to a line, trying to maximize the distance of the centers of two classes and minimize the dispersion of points within the same class. The covariance matrix is used to measure the dispersion within a class. Here, we took ‘One-vs-Rest’ (OvR) strategy to construct multi-classifier.

We compared the predictive performance of these models by tracing their overall accuracy. Overall accuracy measures how often a machine learning model correctly predicts the outcome. We calculated the overall accuracy by dividing the number of correct predictions by the total number of predictions. To further evaluate our models, 5-fold cross-validation was performed. Briefly, we randomly divided the data into five sets with approximately equal size, and used four of the five sets as the training set and the remaining set as the testing set to identify the positives and negatives. We considered precision and recall for specific cancer type i: Precision i = N u m b e r   o f   s a m p l e s   c o r r e c t l y   c l a s s i f i e d   a s   c a n c e r   t y p e   i N u m b e r   o f   s a m p l e s   c l a s s i f i e d   a s   c a n c e r   t y p e   i Recall i = N u m b e r   o f   s a m p l e s   c o r r e c t l y   c l a s s i f i e d   a s   c a n c e r   t y p e   i N u m b e r   o f   s a m p l e s   o f   c a n c e r   t y p e   i

Due to the high dimensionality of DNA methylation profiles, the dimensionality reduction step is necessary before the classifier construction. Principle component analysis (PCA) is a statistical procedure to reduce the dimensionality of a dataset with a large number of interrelated variables by creating a new set of variables called principal components. The greatest variance by some projection of the data comes to lie on the first coordinate (w1), the second greatest variance on the second coordinate (w2), and so on. The principal components were selected based on cumulative percentage of total variations. We selected the number of principle of components taken together explaining more than 95% of the variance.

At first, we identified tissue-specific DNA methylated sites to reduce the considerable redundancy of the original data. We calculated differential methylation values (β value) of CpGs for the corresponding cancer type compared with other cancers using Student’s t-test with a threshold of F.D.R. < 0.01. Next, a recent proposed feature selection method was employed, named Maximum-F-statistic-Maximum-Distance (MFMD), to further detect the tissue-specific CpG sites. Briefly, we calculated the analysis of variance (ANOVA) to compare the DNA methylation levels among cancer types. In ANOVA, F-statistic is the ratio of the variance among the means to the variance within the samples. F-statistic is used to measure the difference among cancers. Euclidean distance (ED) was used to measure the data redundancy. The criterion of MFMD is redefined as follow: M F M D = ⁡ max   w s × F s t a t i s t i c + w d × E D the variable w _s (0 < w _s ≤1) and w _d (0 < w _d ≤1) are the weights of F-statistic and distance, respectively. We ranked the CpG sites according to the MFMD values. The final feature set will have lowest ED values and highest F-statistic values. Then, top-ranked CpG sites were selected as features to construct classifiers and evaluate the classification accuracy. The top-ranked CpG sites with highest accuracy were selected as the final features.

We next attempted to determine whether tumor tissue origin can be predictive using small number CpG markers. Selecting true tumor tissue specific features is important to construct classifiers that performs well at predicting tumor origin sites. To this end, we proposed a method called Maximum-F-statistic-Maximum-Distance (MFMD) to measure the tumor tissue specificity and redundancy of CpG sites. The feature candidates were ranked based on MFMD score, and the top-ranked features were used to construct classifiers to evaluate the classification accuracy. Figure 2 summarized the performance of the classifiers as a function of the number of CpG markers. The result showed a sharp increase in the overall accuracy of the classifiers at the initial stage when the number of CpG sites is small. There are diminishing increases of overall accuracy with the involvement of additional CpG sites. When the number of CpG sites used reaches 1,000 (Supplementary Table S1), it is enough to achieve an overall accuracy >90% and the overall accuracy of the classifier starts to level off. This result indicated that a competitive performance of the classifiers can be achieved using a small number of CpG markers. We further measured the locations across the chromosome of these CpG sites, and found that most of the CpG sites are located at introns (45.7%) and promoters (21.4%).