{'fraction_non_alphanumeric': 0.02790103297717511, 'fraction_numerical': 0.02058854850733023, 'mean_word_length': 4.927203871239217, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A new parsimonious method for classifying Cancer Tissue-of-Origin Based on DNA Methylation 450K data

DNA methylation is a well-studied genetic modification that regulates gene transcription of Eukaryotes. Its alternations have been recognized as a significant component of cancer development. In this study, we use the DNA methylation 450k data from The Cancer Genome Atlas to evaluate the efficacy of DNA methylation data on cancer classification for 30 cancer types. We propose a new method for gene selection in high dimensional data(over 450 thousand). Variance filtering is first introduced for dimension reduction and Recursive feature elimination (RFE) is then used for feature selection. We address the problem of selecting a small subsets of genes from large number of methylated sites, and our parsimonious model is demonstrated to be efficient, achieving an accuracy over 91%, outperforming other studies which use DNA micro-arrays and RNA-seq Data . The performance of 20 models, which are based on 4 estimators (Random Forest, Decision Tree, Extra Tree and Support Vector Machine) and 5 classifiers (k-Nearest Neighbours, Support Vector Machine, XGboost, Light GBM and Multi-Layer Perceptron), is compared and robustness of the RFE algorithm is examined. Results suggest that the combined model of extra tree plus catboost classifier offers the best performance in cancer identification, with an overall validation accuracy of 91% , 92.3%, 93.3% and 93.5% for 20, 30, 40 and 50 features respectively. The biological functions in cancer development of 50 selected genes is also explored through enrichment analysis and the results show that 12 out of 16 of our top features have already been identified to be specific with cancer and we also propose some more genes to be tested for future studies. Therefore, our method may be utilzed as an auxiliary diagnostic method to determine the actual clinicopathological status of a specific cancer.


DNA methylation and cancer identification
DNA methylation is an epigenetic modification of the genome that is crucial for the normal regulation of gene transcription [1]. It is a biological process involving the transfer of a methyl group onto the C5 position of the cytosine to form 5methylcytosine [2]. In normal cells, this modification results in different interaction properties assuring the proper regulation of gene expression and of gene silencing [3]. whereas hypermethylation of the CpG islands is associated with silence of cancer suppressor gene [4].) With the DNA methylation profile, it is expected that classification of cancer can be achieved based on different levels of DNA methylation, which may be helpful in tumor diagnosis and drug development [1,5].
Generally, it is known that cancer is a group of diseases that are driven by progressive genetic abnormalities that include mutations in tumor-suppressor genes and oncogenes, as well as chromosomal abnormalities [6]. While the role of genetic mutations has been highlighted by many researchers [7][8][9][10][11], epigenetic alternations such as DNA methylation have also been found highly related to cancer development [12].
Traditionally, cancer identification is carried out through medical tests such as lab tests, imaging tests and biopsy [13]. Previous work of Masilamani et.al [14] suggests that lab tests such as blood tests for tumor markers have great potential for early detection of cancer. However, as pointed out by National Cancer Institutes (NCI), abnormal lab results are not a sure sign of cancer and most of them may only be helpful in assessment of cancer risk rather than disgnosis [15]. Likewise, imaging tests are helpful in the early detection of cancer. They can help locate and find out the stage of the tumor. However, due to some similar appearance between cancer and other types of diseases, they are not suitable for making a final diagnosis.
Additionally, since it takes millions of cells to make a tumor big enough to appear on an imaging test, imaging test may only be helpful in detecting large groups of cancer cells [16]. In terms of biopsies, although they provide the most accurate diagnosis of cancer, a sample is needed to run such tests. Usually, the biopsy sample is obtained with a needle, endoscopy or even a surgery and these may cause great pain to patients [13].
Considering the limitations of current techniques, we hope to provide additional tools for cancer indentification and classification from a molecular level using sequencing DNA methylation data. The genetic analysis of methylated DNA data may also provide new insights into drug development.

Sequencing data improves cancer classification
Accurately identifying specific cancer holds great promise for further treatment selection and prediction of prognosis in cancer [12]. However, due to the limitions of surgical tolerance in patients, it might be difficult for the surgeons to use complex anatomy to obtain cancer diagnosis in high accuracy [17]. Sometimes, even when the tumors are successfully obtained, similar histopathological appearance of different tumors may still lead to misdiagnosis and consequently treatment may not achieve effective therapeutic effect [18]. The uncertainty of this kind of prognosis propels the development of other diagnostic aprroaches with higher certainty. With the improvement of the Next-Generation Sequencing (NGS) technology, more and more sequencing data e.g., DNA-seq, RNA-seq are available to researchers which require advanced methods to mine these data and provide a more general and reliable information for clinical diagnosis [19] through specific biological insights. A great number of research have been carried out for classifying cancer types and subtypes using sequencing data. For instance, Li(2017) and his team, who used the RNA-seq Data to extract 20 features in the classification of 33 different types of tumors with the mean of ML, with the final accuracy rate exceeded 90% [20]. Similar thoughts could be seen from Zhang (2018)'s project, where they proposed a prognosis model based on Bayesian network classification of breast tumor with 888 features extracted from DNA methylation Data [21]. Raweh (2018) and his team used TCGA data sets and RnRead software to extract the 512 useful feature from the raw data, then utilized the method of fast Fourier transform algorithm for seven types of cancer respectively extracted less than 10 corresponding characteristics and finally used different classifiers (such as random forest) and support vector machine to get an overall accuracy which was above 97% [22]. However, despite the fact that many techniques have been applied to analyze gene expression data, it has been demonstrated by previous studies that some common machine learning methods such as Support Vector Machine may not perform well due to the high dimensionality and small sample size of Gene Expression Profiles (GEP) dataset [23].
Thus, the design of dimensional reduction method to drastically reduce the number of features in sequencing data would be the key to acquire a model in high accuracy.

Our work
In this study, we use the DNA methylation 450k data from The Cancer Genome Atlas to evaluate the efficacy of DNA methylation data on cancer classification for 30 cancer types. Since the dimension of the data is extremely high(over 450 thousand ), we first use the variance selection method to reduce dimension and then utilize the recursive feature elimination algorithm to extract the features. Enrichment analysis is then carried out to explore the biological meaning of the selected genes.

Motivation
The motivation of this paper is two-fold: (i) First, inspired by Darst who acknowledged the positive impact of Random Forest-Recursive Feature Elmination algorithm in extracting features from smaller data sets, our team apply the algorithm into high dimensional data [24].
(ii) Second, the 450k methylation value provides a genome-wide quantitative representation of DNA methylation value in a convenient format, which makes the data preprocessing process for analysts easier and more feasible and thus arouses our interest to explore its correlation with cancer [25].

Novelty
The novelty of this paper is three-fold: (i) First, many previous works [20,21] dived into the area where the sequencing data was used to diagnosing cancer type, but to the best of our knowledge, none of them focused on the DNA methylatyion 450k data to carry out the cancer classification.
(ii) Second, instead of building a model with hundreds of identifiers like  and Raweh (2017)did, our team focus on the feasibility to construct a parsimonious model with relatively fewer number of identifers (within 50 features).
(iii) Third, we prepose a novel method which includes variance selection, recursive feature elimination and boosting classifiers in extracting features from high dimensional data.

Key contributions
The key contributions of this paper are two-fold: First, we have demonstrated the efficiency of our parsimonious model, with an overall accuracy of 91% , 92.3%, 93.3% and 93.5% for 20, 30, 40 and 50 features respectively, achieving the same level of accuracies of over 90% as other studies which use DNA micro-arrays and RNA-seq Data.
Second, 12 out of 16 top features in our studies have been already shown to play a part in cancer development, and we thus recommend some more genes for further testing.
This paper is divided into four sections as follows: 1) Section I is the introduction. It briefly describes the background of the cancer identification and the significance of DNA methylation for cancer identification .
2) Section II contains methodology. This part mainly introduces the source of the dataset and how our work is conducted in the two phases: dimensionality reduction and feature selection.
3) Section III Results and Discussions. This chapter presents the result of our alogirthm and compared the performance of models with different number of features and classifiers. At the same time, it also includes enrichment analysis and gene annotations to explore the biology meaning of our works 4) Section IV contains the conclusion of this article and the work we expect to do in the future.

Dataset
The Cancer Genome Atlas (TCGA) is a landmark project led by the National Cancer Institute (NCI) and National Human Genome Research Institute (NHGRI) that aims to catalogue major cancer-causing genetic mutations [26]. The TCGA dataset, molecularly characterized over 20,000 primary cancer and matched normal samples spanning over 30 cancer types, is publicly available [27,28].
This study focuses on the DNA methylation data from the TCGA dataset provided by UCSC Xena. DNA methylation profile was obtained using the Illumina Infinium HumanMethylation450 platform which includes more than 450 thousand CpG site probes(HumanMethylation450), which provides quantitative methylation measurement at the single-CpG-site level. Featuring more than 450 thousand methylation sites, the Infinium HumanMethylation450 BeadChip offers a combination of comprehensive, expert-selected coverage and high throughput at a low price, thus enabling cost-effective DNA methylation analysis for a variety of applications [29]. Due to its comprehensive genome-wide coverage as well as incomparable cost-effectiveness, the HumanMethylation450 datasets have become a popular choice for the study of epigenetic changes in many diseases processes [29][30][31]. Saumya used DNA Methylation 450K data to reveal the peripheral blood differential methylation in male infertility and Bonnie utilized the same data to make the MGMT mrthylation assessment in glioblastoma [32,33].
In this study, the methylation 450kdata we use derived from UCSC Xena has already gone through quality control such as filtration of low quality probes as well as normalization process and DNA methylation level is measured in beta value defined as the ratio of the methylated probe intensity and the overall intensity (sum of methylated and unmethylated probe intensities) [34].
which is a continuous variable between 0 and 1. Higher beta value represents higher level of DNA methylation whereas lower beta value represents lower level of DNA methylation.

Dimension Reduction
The study focuses on 30 DNA methylation datasets extracted from TCGA (listed in Table 1). The original data consists of over 450 thousand identifiers representing the ID number of the probes and 9743 samples in total, with a number of features containing missing value due to low quality of certain Methylation450 BeadChip probes.
Inspired by the idea of neural network, the feature selection from the original 480,000 CpG sites is implemented in two-phases ( Figure 1). we propose a 'hiddenlayer' (Medium Layer)size features which is the output of the variance selection (Dimension Reduction) as well as the input of the latter RFE [35] process (Feature Selection).
In the first step, by variance selection method (Dimension Reduction), we reduce the number 450,000 to 2000 features.
In the second step, the first phase's output becomes the input of the RFE process (Feature Selection) and we extract 20 features from the 2000 CpG sites, which will be later discussed in 2.3. The reason we construct a 'hidden-layer' here is that though RFE can provide precise solution, the time overhead is a defect when the input features are enormous [35]. In addition, if we direcly apply RFE algorithm to 480,000 featrues, the algorithm's stability is questioned since redundant and noisy genes are not removed at first as Tang proved [36].
Before the dimensionality reduction and feature selection, features with missing values in the samples were removed and then 30 datasets were merged into one overall dataset.
After that, Dimension reduction was conducted based on variance between and within different cancer groups . Afshar(2020) demonstrated the effectiveness of removing features with low variance so as to select important features in high dimensional data [37]and Model(2001) showed that genes with large variance are more important features since they explain most of the total variances [38].
Therefore, we combined their thoughts and applied the variance selection in this study.
The following are the formulas of the within groups variance !,# $ and between

Features selection
In this part, we mainly utilize the feature selection algorithm recursive feature elimination method (RFE) method [35] to obtain the most important features from a total of 2000 identifiers after the first step. As the work of Guyon et.al [35] illustrates, this algorithm achieve high accuracy in genetic diagnosis and gene selection. More specifically, RFE is a method for feature selection which selects features by recursively obtaining increasingly smaller sets of features until gets the required number of features. An estimator, namely a classifier gives the impotrance of the features, is an essentail part of this algorithm, which computes the significance of the features. In each loop, it ranks the significance of the features and eliminates the worst feature. Here, we choose accuracy as the indicator for the RFE and mainly tune the number of chosen features to compare the results.
Since the RFE method is estimator bounded, the choice of it can be vital. According to Deng(2012), tree-based classifiers are frequently chosen as the estimator since they can provide variable significance scores and perform relatively strong [39].
Similar points can be seen in Touw(2013) 's paper, which demonstrates the significance of random forest classifier to the bioinformatics area [40]. Therefore, random forest are chosen in this paper as the estimator and we also proposed its similar algorithm decision tree [41] and extra tree [42] to be tested.
On the other hand, Moon(2016) demonstrated that the L1-norm SVM performed well as the feature selection estimators when dealing with the biomedical data [43].
Reasons are that the figure for the biostatisticians tend to be centered in an interval which can be hard to distinguish while L1-based classifiers are able to make them sparse and easier for detect. As a result, since the data we have are the decimals between 0.0 and 1.0(mainly 0.1 to 0.9 according to the information on https://xenabrowser.net) , we follow Moon and his team to use the L1-norm SVM as the estimator Therefore, we aim to use these four estimators to get the scores for each feature and then select the top several to feed into the classifiers. To mention above, fivedifferent classifiers are used in this paper. Following Sahu(2017) 's research [44] on the DNA microarray data, Support Vector machines (SVM) [45], k-Nearest Neighbors (KNN) [46] and Multilayer Perceptron (MLP) 2.3 [47] are chosen. Besides, another two ensemble classifiers LightGBM [48] and Xgboost [49] are also used.
The feature set filtered by variance was first divided into 10 parts for 10 folds cross validation( Figure 2). In each iteration, keep one fold of the data as the validation set and use the other 9 parts as the training set. Then we use the RFE method to fit on the training set and achieve the selected features. Thus, a new subset of the data, written by using the selected CpG sites as features, enters into the classifier to be evaluated. Here, the key point is that we need to check the performance of the selected features as well as choosing the best group (features plus the classifier) as the final model.
Since there are four estimators as well as five various classifiers in the experiment, a total number of 20 groups are made and each group contains 10 iterations since a ten folds cross-validation is contained in each group. As shown in the graph( Figure   3), the training set and test set were tested respectively with the models (the motion for testing the training part is to check whether the model is overfitting), and the accuracy was recorded at intervals of 10 features from 20 to 50. In this way, two lists each containing train score and test score respectively are obtained. Therefore, for each of the 20 models, we obtained accuracy for each folds on the validation set and averaged them to get one mean set accuracy for the different number of features as well as calculated the minimum and maximu value of the 10 groups.

Results and Discussions
In this section, we apply our feature selection method to extract important features from the dataset and then analyze the results. The biological meaning of the figure will also be included.

.1 Extra Tree as the Estimator
In the methodology part, we have mentioned that an overall 20 models are used to look for the best estimator for the RFE algorithm, and the results are as followed. We did the 10 fold cross validation to each of these combined models and same to other paper [50], the accuracy here is based on its performance on the test set. The number of features we first selected is 20 as our novelty here is to look for models with least features but also keep a relatively decent accuracy.

Figure 4 test accuracy on each model
From Figure 4, extra tree algorithm is selected as the estimator of the RFE algorithm since the model with the extra tree have higher 10 cross validation accuracy compared to others. Also, since boosting classifier as lightgbm and xgboost performance well with our model, we additionally test the score of the Catboost classifier [51].

Model performance for 20 points
We first make our experiment on 20 points since it is the proposed lower bound for our feature numbers. According to 1.1, we obtain the 20 features as below.

Max Min Avg
In the research, the sklearn metrics to access the performance of the model includes true positive, true negative, false positive and false negative. Several criteria such as recall score, precision score, f1 score and the overall accuracy will be used as the main indicator for the paper. The task is a multi-class classification, so we choose the parameter 'weighted' to be the key parameter for the indicators since the samples of each class in the dataset is not totally the same. In the following formulas, * represents the number of sample that the truth is class i and the test predicts class i ; * represents the number of sample that the truth is not class i and the test predicts not class I ; * represents the number of sample that the truth is class i but test predicts not class i ; * represents the number of sample that the truth is not class i but test predicts class i [52]. . The F1 score is the harmonic mean of the precision and recall. [52][53][54] And then we get the model performance of the three classifiers with only 20 features based on these indicators.   Figure 5).

Figure 5 Confusion metrics for 20 points
The element in the diagonal of the matrix represent the sample which can be predicted correctly from the model, and all the other element in the confusion metrics are the wrong predicted samples [55]. On the other hand, similar to Li(2017)'s work, though they achieve an overall high accuracy on the test set, they still fail to classify three tumor types, such as READ(rectum adenocarcinoma) and COAD(colon adenocarcinoma).
Our team here again showed the difficulty of classifying between READ and COAD class and we may test whether increasing the features can solve the problem in 3.3 From figure 6 , the performance of the model increases continuously until reaching 40 points, the slope here starts to decrease. Since our aim is to look for the solution with best efficiency (the lowest number of features with high accuracy), we prefer the solution with 40 features. We then analyze the confusion metrics of the models as below (Figure 7, Figure 8 and Figure 9). The UCS class has an obviously higher accuracy compared to 20 points and so does the ESCA class, so we point out the possibility that these two classes can achieve a decent overall cross-validation accuracy with enough number of features. However, since our paper mainly focuses on classifying these cancer types based on tiny number of features, we will not include further discussion here. Meanwhile, 21 of 30 class with 40 features can achieve accuracy over 0.9 compared to 17 for the 20 features and. In a nutshell, our solution using only 40 identifiers is able to predict the 30 cancer types with 93.3% overall accuracy using the data of DNA methylation 450 values.

Test robustness of the model
This section works as an inspection for the model, we need to examine the convergence of the RFE algorithm.

RFE Inspection
First, our team examine the intersection of the features. Since RFE algorithm is a backward stepwise selection method which means that the solutions of fewer feature should better be included in the larger feature set. Then we calculate the exact intersection set as followed.  Table 3

elements of intersection set
Therefore, the convergence of our model is totally good since the inclusion rate is high. To be mentioned, since the 2000 filtered identifiers need to be reduced to 20, the first twenty selected features should be most robust, then the second twenty features . So it is reasonable that the overlapped features of 20 points with the 50 points are greater than the 20 points with the 40 points.
Second, we dive into the performance of the intersection points. points) do follow the trend. One point to be noticed is that the indicators decrease for features fewer than 20, which demonstrate that 20 features can be a good starting point for the experiment.

Enrichment Analysis
After extracting the identifiers(CpGs) from the data, we map back to find the related genes(Prefix 2) and enrichment analysis has been carried out with DisGeNET [56] and GO Biological Processes [57]by using the metascape [58].

Figure 12 Network of enriched terms colored by cluster ID
Next, we identified all statistically enriched terms accumulative hypergeometric pvalues and enrichment factors were calculated and used for filtering. Remaining significant terms were then hierarchically clustered into a tree based on Kappastatistical similarities among their gene memberships. Then 0.3 kappa score [65] was applied as the threshold to cast the tree into term clusters. We then selected a subset of representative terms from this cluster and convert them into a network layout( figure 6). More specifically, each term is represented by a circle node, where its size is proportional to the number of input genes fall into that term, and its color represent its cluster identity. Terms with a similarity score > 0.3 are linked by an edge [66][67]. As can be seen from figure 6, it is clear that the first polygan gather most terms and thus need a deeper analysis. The table below (  Table 4 Detailed information of enriched terms in Figure 12

Gene Annotation
In this section, we will give gene annotations to the overlapped features(intersection of 20, 30, 40, 50 solutions  family of proteins. The encoded protein regulates the activity of various transcription factors, including the androgen receptor, Smad3/4, and p53 [69]. Rogers(2013) demonstrated that the ZMIZ1 might confer selective advantage to tumor cells, which indicated the casual link between ZMIZ1 and tumorigenesis. [70] GBAP1 is a pseudogene which is defined as nonfunctional gene bearing close resemblance to the known genes at different loci [71]. Research has been conducted to uncover the epigenetic regulation, biological function, and clinical application of GBAP1 in Gastric Cancer risk [72].
HOXB1 plays a important role as an oncogene in diverse tumors as result suggests that HOXB1 gene is a tumor suppressor which is regulated by miR-3175 in Glioma. [73] FGF18 contributes to breast cancer cells following activation of the Akt/GSK3β/βcatenin pathway, which makes FGF18 a possible candidate target for breast cancers. [74] TBC1D16 have been found associated with tumor progression and metastasis in multiple cancer types. Hypomethylation of TBC1D16 was observed in breast cancer and metastatic tumors [75].

SH3PXD2B encodes an adapter protein that is characterized by a PX domain and four
Src homology 3 domains [76]. The encoded protein is required for podosome and invadopodia formation and is involved in cell adhesion and migration of numerous cell types. It has been demonstrated that inhibitors of podosome and invadopodia formation might have utility in the treatment of vascular diseases and cancer [77].
HOXA9 which is down-regulated in lung cancer and acts as a tumor progression suppressor is closely related to the aggressive growth of lung cancer cells [78].
HOXA9 contains a homeobox domain that binds to DNA and regulates downstream gene expression [79]. In addition, it is demonstrated that exogenous up-regulation of HOXA9 inhibits tumor cell invasion and migration, inhibits the expression of zinc finger 2 (SNA12/SLUG) by inhibiting the activity of the nuclear factor (NF) -kb [80].
PRKCE is generally related to the expression of protein kinase C epsilon (PKCɛ).
Encoded by PRKCE, PKCɛ is an enzyme often related to cell transformation and tumorigenesis. Under the action of PKCɛ, the Ras/Raf pathway is activated, resulting in the transcription of genes involved in cell proliferation and growth [81]. there is significant DNA hypermethylation in the EPAS1 regulatory region, which is associated with a decrease in EPAS1 mRNA level in primary cancerous tissues [82].
NKX6Bis a murine-homeobox-containing gene localized distally on Chromosome (Chr) 710q26, a region where frequent loss of heterozygosity has been observed in many brain tumors.It is suggested that NKX6B is a possible candidate tumor suppressor gene for brain tumors, particularly for oligodendrogliomas. [83] PLAC8 is a multi-functional protein that is highly expressed in the intestine, lung, On the other hand, our team also pointed out the problems of distinguishing between COAD(colon adenocarcinoma) and READ(rectum adenocarcinoma) class using the methylation value, which agrees with Li (2017). Since Li uses the RNA-seq data which is somewhat different from us, it is necessary for more research to be done into the area in the future.