Identification of Targetable Pathways in Oral Cancer Patients via Random Forest and Chemical Informatics

Treatment of head and neck cancer has been slow to change with epidermal growth factor receptor (EGFR) inhibitors, PD1 inhibitors, and taxane-/plant-alkaloid-derived chemotherapies being the only therapies approved by the U.S. Food and Drug Administration (FDA) in the last 10 years for the treatment of head and neck cancers. Head and neck cancer is a relatively rare cancer compared to breast or lung cancers. However, it is possible that existing therapies for more common solid tumors or for the treatment of other diseases could also prove effective against oral cancers. Many therapies have molecular targets that could be appropriate in oral cancer as well as the cancer in which the drug gained initial FDA approval. Also, there may be targets in oral cancer for which existing FDA-approved drugs could be applied. This study describes informatics methods that use machine learning to identify influential gene targets in patients receiving platinum-based chemotherapy, non-platinum-based chemotherapy, and genes influential in both groups of patients. This analysis yielded 6 small molecules that had a high Tanimoto similarity (>50%) to ligands binding genes shown to be highly influential in determining treatment response in oral cancer patients. In addition to influencing treatment response, these genes were also found to act as gene hubs connected to more than 100 other genes in pathways enriched with genes determined to be influential in treatment response by a random forest classifier with 20 000 trees trying 320 variables at each tree node. This analysis validates the use of multiple informatics methods to identify small molecules that have a greater likelihood of efficacy in a given cancer of interest.


Introduction
Head and neck cancers can be defined as cancers in the upper airway and/or digestive tract found in oral cavity and laryngeal, pharyngeal, oropharyngeal, and hypo-pharyngeal tissues. Such cancers make up 3% of cancers diagnosed each year. 1 Head and neck cancer incidence has declined from 25 cases per 100 000 in the 1990s to 15 cases per 100 000 in the present day. 2 While the decrease in head and neck cancer incidence may be due to a drop in tobacco use, 3,4 the mortality associated with head and neck cancers has not changed significantly in the last 20 years. Treatment of head and neck cancer has also been slow to change with epidermal growth factor receptor (EGFR) inhibitors, PD1 inhibitors, and taxane-/plant-alkaloid-derived chemotherapies being the only therapies approved by the U.S. Food and Drug Administration (FDA) in the last 10 years for the treatment of head and neck cancers. Head and neck cancer is a relatively rare cancer compared to breast or lung cancers. However, it is possible that existing therapies for more common solid tumors could also prove effective in oral cancers. Many therapies have molecular targets that could be appropriate in oral cancer as well as the cancer in which the drug gained initial FDA approval. Also, there may be targets in oral cancer for which existing FDA-approved drugs could be applied. This study describes informatics methods that use machine learning to identify influential gene targets in patients receiving platinum-based chemotherapy, non-platinum-based chemotherapy, and genes influential in both groups of patients.
Drugs approved by the FDA for oral cancer are methotrexate, cetuximab, pembrolizumab, nivolumab, and docetaxel. These therapies are combined to be used in conjunction with platinum-based chemotherapies such as cisplatin or carboplatin unless those therapies are contraindicated due to comorbidities such as renal disease. The small number of new oral cancer drugs could be attributed in part to the low overall burden of oral cancer in comparison with other cancers. The current timeline for FDA approval of a novel small molecule or biologic is 10 years or more. Repurposing existing FDA-approved drugs is a popular method used to shorten this process to 3 to 4 years. Cetuximab is an EGFR inhibitor that has been shown to decrease the rate of progression of oral cancer used in conjunction with Cisplatin. Pembrolizumab and nivolumab are both PD1 inhibitors that used T-cells to attack cancer while it progresses. These therapies use the body's immune system as another treatment modality to reduce the burden of oral cancer. Current literature provides support for the role of ligand channel gating, hedgehog signaling, 5-7 NOTCH, B-WICH, 8 inflammasome, 9 WNT, 10 and Calci neurin pathways in cancer. [11][12][13] The role of these pathways and targeting specific genes within them has been pursued in other cancers, but have with few exceptions, not yet been examined in oral cancer. Possible gains from targeting these pathways would be initiating immune response, targeting cancer metabolism, targeting signaling for metastasis, and targeting inflammation pathways that may drive progression. If there is a synergistic effect to attacking 2 Cancer Informatics multiple hallmarks of cancer simultaneously, then the net gain to the patient would be a net gain in overall post diagnosis survival time for the oral cancer patient.
Identifying a means by which drugs may be prioritized for further screening and validation for a specific cancer type would be desirable. Databases linking genes, the proteins they express, ligands corresponding to those proteins, and structural data that can be analyzed all exist in varying forms or completeness across different publicly available databases. This study describes how integration of analyses of these databases can be used to select gene targets in a specific cancer and how therapies can be prioritized for screening based on existing structural information for the ligands associated with genes and the proteins they express.
There are several hurdles to the analysis of high-dimensional genomic data using traditional regression analyses. Random forest analysis is a machine learning approach that is less hindered by datasets with large predictor to observation ratios. In this study, we apply random forests analysis to gene expression data to identify those genes and pathways that are most predictive of post diagnosis survival across treatment strata. This is the first application of this approach to head and neck cancer patient data in The Cancer Genome Atlas (TCGA). National Comprehensive Cancer Network (NCCN) 14 guidelines recommend that node-positive patients with tumors of clinical stage 3 and greater receive chemotherapy. Platinum-based chemotherapy with radiation and surgery is the current standard of care recommended to these patients. Patients that do not receive chemotherapy recommendations by NCCN are node negative with clinical stage 2 and lower. Following standards set by NCCN guidance on treatment, this study chose to identify patients receiving and not receiving platinum-based chemotherapy as separate groups. Analysis of influential genes in each group will improve knowledge of possible mechanisms driving treatment response for early stage and more advanced tumors.
Random forest is a machine learning approach to identifying the most important predictors in high dimensional datasets.15,16 This approach is uniquely suited for classification of observations in datasets where P (predictors) are greater than N (number of observations). Random forest randomly selects predictors from a large group of predictors and then applies those predictors to a decision tree predicting overall survival. Random forest does not pay a statistical penalty when the number of observations is small. Instead the strength and limitation of this method is its reliance on computational intensity. That is, as the number of decision trees in a random forest increase, so does classification accuracy. Accuracy is also dependent on the number of predictors tried at decision tree nodes. As node size and forest size increase, so does forest classification accuracy. However, there is a rate of diminishing returns in the accuracy gained from each tree added to a forest. Therefore, computational time and cost must be factored into all random forest analysis plans to measure project feasibility. Random forest has been successfully applied to predicting cancer diagnosis and treatment response for a variety of cancers. [17][18][19][20][21] For this study, we have selected to apply random forest analysis to the gene expression values of oral cavity cancer patients to identify the upregulated pathways most predictive of improved treatment response across gender and environmental exposure subgroups like alcohol and tobacco. RNAseq data are inherently high dimensional, applying typical regression models to such data can be costly as large sample sizes are required to identify even moderate effect. Identifying gene interactions can be even more costly in terms of the required statistical power. Stratified pathway analysis via random forest methods has been shown to be successful in identifying single influential genes (within the context of larger pathways) that are predictive of overall survival with limited sample size. 22 This approach has not yet been applied to identification of influential genes and gene interactions within oral cancer patients stratified specifically by treatment. In this way, the importance of pathways and genes of interest can be compared across strata to assess which subgroups may be most sensitized to changes in gene expression within a given pathway.

Methods
This study focuses on the identification of the role of gene expression in oral cavity cancer patients and applying machine learning approaches like random forest to determine genes that are important in influencing treatment response. Reference ligands known to bind to proteins expressed by genes deemed influential by random forest can be sent through a virtual screening pipeline to identify small molecules with greater likelihood of acting as protein agonists/antagonists. Ligands that have a strong shape similarity to known binding ligands have greater potential for success in high-throughput screening endeavors. As shape similarity alone is insufficient in identifying new drug leads, all leads will be validated with existing literature, and those leads without previous biological validation will be presented as such.
By using a stratified random forest analysis, we will be able to rank genes within the strata of chemotherapy treatment status. This approach will allow for the identification of those top ranked genes that are unique to each stratum. This will be done by identifying common and unique genes between sets of genes influencing the treatment response in patients receiving platinum-based chemotherapy and those that do not. The result will be the identification of oral cavity cancer pathways influencing treatment response which will inform researchers on mechanisms driving treatment response in specific groups such as late-stage, node-positive patients who are more likely to receive chemotherapy treatment. This analysis will illustrate and support existing studies showing the strength of machine learning methods as an alternative method in identifying gene expression values influencing treatment response. This study is focused not only on the predictive power of an aggregated panel of gene expression values, but rather to integrate random forest with chemical informatics and thus describe methods to shorten the pursuit of novel therapies treating cancers with relatively lower incidence.

Retrieval of Public Data
This study used clinical and genetic data obtained from TCGA. Genetic data included raw counts per million (CPM) of RNA sequence expression values for 523 patients posted to TCGA. Of these 523 patients, 313 were diagnosed with Oral Squamous Cell Carcinoma (OSCC). Oral squamous cell carcinoma patients included tongue, buccal mucosa, alveolar ridge, general oral cavity, and soft palate tissues. Of these 313 patients, 267 were included based on complete survival time. Of these 267, 109 received either Carboplatin, Oxaliplatin, or Cisplatin, while 158 patients received a treatment other than the platinum-based chemotherapy treatment regimen. All tissue samples were collected prior to start of treatment. Clinical data on tumor stage, necrosis, size, and nuclei were also retrieved from TCGA. We obtained demographic data on ethnicity, race, and gender from TCGA clinical files. In addition, this dataset had information on environmental exposures like tobacco history (ever/never smoke) and number of alcohol containing drinks in a day (greater than 2 drinks consumed per day, 2 drinks or less consumed per day). Overall survival time in months was extracted as a measure of treatment response.

Description of random forest approach
Stratified pathway analysis considers important covariates in data analysis. In this analysis of head and neck cancer gene expression data, this study used information regarding patient age, sex, alcohol, and smoking exposures. In the first round of analysis, the sample of 267 OSCC patients was organized into 109 and 158 subsamples based on whether patients received platinum or nonplatinum-based chemotherapy, respectively. For each group, this study built a random forest 15 to predict survival time based on the gene expression levels within each subgroup.
To better communicate the function of random forests, understanding of a decision tree construction is needed. A decision tree is constructed by the following steps: Step A: A bootstrap sample was taken from the original sample.
Step B: A decision tree is grown for each bootstrap sample.
Step C: At each tree node, a predetermined number of predictors were applied randomly to create branches within the tree.
Step D: A branch is formed using the predictor from step C.
Step E: Steps C and D are repeated until the end of every tree branch contains samples above or below the same survival threshold or contain only one sample.

Random forest result measurements
Random forests build many decision trees to comprise a forest. Each tree is put together by using a random bootstrap sample of the original data and applying a random number of predictors at each node of a decision tree. The SRCRandomForest R package 23 employed in this study sets aside half of the data to be used for validation purposes to measure the accuracy of the random forest model. The SRCRandomForest package was chosen for this study because it applies a multivariate Cox regression model to produce each decision tree within each Forest. In this way, censored and non-censored data can be considered when carrying out this analysis. In this way, the model can measure the degree of influence of each gene on patient survival. The P values yielded through this analysis are defined as the "proportion of cross-validation errors smaller than the cross-validation errors obtained from 500 iterations of random forest runs of randomly permuted labels of patients." This list of genes can be used to identify pathways that are enriched with the influential genes identified through random forest at odds that would be greater than can be attributed to chance alone with a P value of .05. This analysis will present pathways common to, and unique to, each chemotherapy treatment strata. Such analyses may identify plausible biological mechanisms that enhance understanding of observed differences in survival. To reiterate, the focus of this study is not to pursue a diagnostic tool but to identify those gene expression values exerting a strong influence on treatment response. This study adopted random forest as the machine learning method of choice due to its superior interpretability and scalability. Random forest and many applications incorporating random forest have been developed more than 30 years. Random forest was of interest to this study as an application developed specifically for the use of survival data had been developed as an open-source R package. Although random forest has been shown to be less efficient than extreme gradient boosting applications, the long history of this approach has allowed it to evolve to produce an application well suited toward achieving the analytic goals of this study.

Random forest tuning parameters
Our random forest model used 20 000 tree forests for a forest size, with 320 variables tried at each node in each decision tree; 20 000 trees were the point at which we could identify no significant increase in our ability to predict patient survival; 320 tries at each node in each decision tree were double the recommended number of tries given by the author of our R software package RandomForestSRC. The author H.
Ishwaran et al 23 refers to a generally accepted practice of using the "square root of the total number of predictors as a starting point for the number of variables tried at each node." It was for this reason that we applied 20 000 trees and 320 variables tried at each node in each decision tree for every stratum in our analysis. This approach was applied to our entire final sample of 267 OSCC patients. We then divided this sample by whether a patient received platinum-based chemotherapy or not. All random forest analyses were carried out on 64 computational nodes in the University of California High Performance Computational Cluster. Random forest analysis was made more efficient by running analysis in conjunction with the R package parallel. The output of each group's analysis produces a list ranking each gene. This analysis identified common and unique pathways between the entire dataset and each chemotherapy treatment group. We then identified the unique and common genes between chemotherapy groups. This analysis will allow us to observe the difference in gene importance and corresponding pathways in relation to overall survival.

Description of virtual screening approaches
Using chemical informatics techniques, ZINC drug database 24 of 1379 FDA-approved drugs (FDA) and ZINC Traditional Chinese Medicine (TCM) database of 39 894 small molecules can be used to apply three-dimensional chemical informatics approaches to the identification of small molecules that are the best candidates for inhibition of proteins expressed by those genes influencing treatment response. Reference ligands for each protein are obtained from the Royal Chemistry Society Protein Database 25 and then virtually screened against FDA and TCM small molecule libraries. Molecular shape overlay is an approach for measurement of the similarity of one molecule in comparison with another. A Tanimoto coefficient is used to measure the degree of similarity between 2 molecules.
The Tanimoto coefficient between 2 points, a and b, with k dimensions is calculated as follows The Tanimoto similarity only applies to binary variables, for binary variables the Tanimoto coefficient ranges from 0 to 1 (where 1 is the highest possible similarity). 26 A goal of this study is to perform searches of 2 small molecule databases FDA and TCM using a maximum common substructure measurement of Tanimoto similarity from the R Rcpi package 27 that has been shown to perform robustly across a variety of molecule types. Table 1 provides information on those genes that are significantly enriched within our gene set beyond what would be expected by chance alone. The significance of enrichment is calculated as the odds of randomly selecting the number of genes in the submitted set of genes by randomly selecting from 20 530 genes more than 100 times. The false discovery rate (FDR) reported by the Pathway Reactome application used for this study represents an adjustment for multiple comparisons across all pathways. Benjamini-Hochberg FDR is calculated as the P value ranking (smallest being 1 and all following having greater or equal rank dependent on size of P value) divided by the number of tests performed and multiplied by the significance criterion. In this study, .05 is the criterion used to measure significance. The FDR can be interpreted as the proportion of tests within a set of tests that falsely rejects the null hypothesis. If an FDR is .5, then 50% of those pathways identified will falsely reject the null hypothesis. It is important to note that the FDR calculation used by Pathway Reactome defaults to perform a large number of analyses/tests. In addition, there are many pathways examining similar genes and gene types. Unfortunately, this thorough examination strategy also inflates the number of analyses and causes the FDR to become overly conservative.

Pathway analysis
To identify those genes that are most likely to be connected to influential pathways, 2 filters were applied to gene selection. First, the gene had to be in the top 5% of influential genes identified via random forest. Second, if the gene was not within the top 5% of genes, it could still be included within the analysis if it was within the top 40% of genes and was known to be connected via past experimental studies supporting gene network connections found in the Cytoscape/Pathway Reactome plugin database. Gene topology can be accessed by uploading topology from a given gene in Cytoscape and merging it with that gene's corresponding random forest importance values. Gene topology refers to the connectedness of a gene (supported by experimental results referenced in cytoscape) to other genes within a Cytoscape network. Thus, the selection of genes included those ranked as the most important by random forest (top 5%) or of moderate importance and high topology (shown through past experiment or literature search to be connected to 50 genes or more). For each set of analyses, 1000 influential genes were selected and 100 high topology genes were selected. In this way, integration of Cytoscape network analysis with random forest results allowed for identification of pathways significantly enriched with the genes identified by the random forest model. For a complete overview of the analyses in this study please reference the online Supplemental Materials.

Top pathways
This section will describe pathways uniquely influential in patient response to platinum-based chemotherapy and influential in Schomberg 5 response to non-platinum-based therapy. Influential pathways shared by both platinum-based chemotherapy users and nonusers will also be presented. An influential pathway will be defined as a pathway that is significantly enriched with genes that were in the list of top 1000 (5%) of most influential genes yielded by random forest analysis for platinum-based chemotherapy users, nonplatinum-based chemotherapy users, and for those pathways enriched with genes shared in common in lists of top 1000 genes for platinum-based chemotherapy users and non-users. The top 5% of genes were conservatively selected to produce greater certainty of the link between highly ranked genes, the pathways in which they were enriched, and the link between treatment response and significant pathways identified through gene enrichment analysis. Top pathways could also be enriched with those genes in the top 40% of important genes identified by random forest if they also had a high amount of connectedness (a gene was connected to 25 gene nodes in a gene network) reported by the Cytoscape/Reactome application.
Pathways significantly enriched with genes identified as the most influential (top 5%) in predicting treatment response for platinum-based chemotherapy users were those related to calcium channel gating, hedgehog signaling, histone acetylation, elastic fiber production, tRNA acetylation, hexokinase deficiency, inhibition of adenylate cyclase, and CLEC7A inflammasome pathway. It has been reported that calcium channel gating has been associated with multiple cancers. [28][29][30] There have also been recent studies evaluating the benefit of targeting histone deacetylation pathways in oral cancer. [31][32][33] The hedgehog signaling pathway has also been shown to signal progression in other cancers "Hh signaling has been shown to regulate the selfrenewal of CSCs in breast, glioma and multiple myeloma, and more convincingly in the maintenance of chronic myelogenous leukemia (CML) stem cells." [34][35][36][37][38][39] All significant pathways for platinum-based chemotherapy users are in Table 1.
Significant pathways for patients not using platinum-based chemotherapy were those related to B-WICH complex, TP53 pathway, fibroblast growth factor receptor (FGFR) pathway, potassium channel gating, and RNA polymerase chain elongation pathways and their epigenetic regulation. TP53 and FGFR pathways represent the expression of canonical oncogenes which have been shown to be cancer drivers and associated with the production of all cancers. [40][41][42][43] The B-WICH complex has been found to be linked to maturation of invadopodium in breast cancer and has been suggested as both a biomarker and target for cancer invasiveness. 8,44 Potassium channel gating has also shown to be a potential target for head and neck cancers due to this pathways association with immune response and treatment response. [45][46][47][48][49] RNA polymerase chain elongation and its role in transcription is a logical contributor to cancer progression and differentiation; however, the lack of specificity makes this a difficult pathway to target specifically in cancer cells. All significant pathways for non-platinum-based chemotherapy users are in Table 2.
For those pathways enriched with genes shared by users and non-users of platinum-based chemotherapy, it was found that pathways related to g-protein beta folding, nuclear factor of activated T-cells (NFAT) activation, and repression of Wnt pathway genes were enriched with genes influential in treatment response in both groups of patients. Repression of Wnt pathway genes may be done through targeting the sonic hedgehog pathway as previously outlined or through more direct means which have been researched in multiple other cancers. [50][51][52] Nuclear factor of activated T-cells proteins have been found to be associated with cancer progression in blood and  [53][54][55] These pathways influencing treatment response in both users and non-users of platinum-based chemotherapy can be seen in Table 3.

Important genes and biological implications
A visualization of pathways overlapping between users and non-users of platinum-based chemotherapy highlights the importance of several genes in a way that random forest analysis alone could not. By visualizing the 4 common pathways, it becomes possible to identify not only highly influential genes but also those genes that have the highest degree of connectivity to influential genes. Using annotation built into Cytoscape, 56,57 we can also identify existing small molecules used in cancer therapy that are not yet commonly used in oral cancer, and we can also observe those genes previously found to be associated with oral cancer. Genes found to be influential in oral cancer for patients receiving platinum-based chemotherapy with existing literature supporting the targeting of these genes in cancer were INSR, BRAF, and PSMB7 which are targeted by Ceritinib (regorafenib and dabrafenib) and bortezomib, respectively.
These drugs are not currently FDA approved for treatment in oral cancer. Genes found to be influential in oral cancer for patients not receiving platinum-based chemotherapy with existing FDA-approved chemotherapy drugs targeting the products of said genes are FGFR3, EGFR, PRKAA2, CSNK2A1, INSR, MET, CAMK2A, PSMB5, and PSMB1. There are multiple chemotherapy drugs targeting these pathways with 14 different drugs targeting EGFR alone. It should be noted that EGFR is a gene pathway being targeted in current oral cancer treatment. Sarafenib Tosylate, Pazopanib Hydrochloride, and Vadetanib all target the FGFR pathway specifically. Sunitinib Malate is unique in that it has been found to act on 4 different genes that were found by random forest to be influential in treatment response: FGFR3, CDNK2A1, PRKAA2, and CAMK2A. Again, we see that Certinib acts on a gene that is influential in both patient treatment groups, and that gene is INSR. Bortezomib, Carfilzobib, and Ixazomib all act on PSMB5, which is an influential gene in both Platinum and non-platinum-based therapy. PSMB5 and PSMB1 are both found to be within the top 5% of influential genes in random forest analysis and are genes that are significantly enriched within the sonic hedgehog pathway. Chemotherapy drugs and their relationship to genes in common influential pathways between users and non-users of platinum-based chemotherapy are visualized in Figure 1 for platinum-based chemotherapy users and Figure 2 for non-platinum-based chemotherapy users.
In addition to analysis of the intersection of existing cancer drugs and genes deemed influential by random forest, this study also looked at the intersection between gene topology within a pathway and random forest influence. This analysis identified CTNNB1, PLCG2, SHC1, UBA52, UBB, UBC, and HDAC3 as genes that meet filters of belonging to one of the 4 common enriched pathways, being a gene that is one of the top 5% of influential genes listed by random forest analysis and being connected to more than 100 genes within the 4 interconnected  Hexagon shapes are genes. Dark red genes are of greater influence based on random forest analysis results (within top 5% of influential genes), white genes do not fall within the criteria of being in the top 5% of influential genes. All 4 common pathways enriched in patients receiving platinum and non-platinum therapy (described in Table 3) were merged into one pathway using Cytoscape. genes were clustered according to the Reactome Pathway Plugin available via Cytoscape.

8
Cancer Informatics pathways (Figures 1 and 2). CTNNB1 mutations have been found to be predictive of lung and other thoracic cancers, [58][59][60][61][62] and PLCG2 and calmodulin knockdown have been shown to induce paclitaxel sensitivity in cervical cancer tumors. This may prove of use to oral cancer patients, which may be assigned to paclitaxel or other taxol regimen. 63,64 SHC1 has been shown to be a regulator of EGFR function and thus a potential target for multiple cancer types where EGFR is a key driver. [65][66][67] Ubiquitin genes UBA52, UBB, and UBC have been shown to be associated with several cancers and research is currently being pursued in targeting ubiquitin ligases to improve treatment response. [68][69][70][71][72] Histone-deacetylase genes specifically HDAC3 is shown to be a hub to several genes that are influential in platinum-based chemotherapy response genes that have also been associated with metastatic invasion in breast and pancreatic cancer ( Figure  1). Inhibition of HDAC3 was shown to impact signaling to cancer stem cells. This gene has been shown to be a regulator of apoptosis control exerted by TP53. [73][74][75][76] These genes were not shown to currently have any antibody or small molecule therapies targeting their action. High-level topology (more than 100 connected genes) and high random forest importance ranking should provide impetus for further research into the targeting of gene action in oral cancer.

Chemical informatics analysis of drug targets and leads
For those genes meeting topology and random forest filters, the known ligands of proteins expressed by each gene were identified though the RCSB protein data bank. Structural files of ligands were downloaded as .sdf files and uploaded into the chemical informatics R package Rcpi. Once loaded, each ligand had to undergo virtual screening against all FDA-approved drugs to identify existing FDA-approved drugs that may prove efficacious as therapeutic agents. Only those molecules with a Tanimoto similarity score >50% were included in results. A TCM small molecule database was also used as biologicderived small molecules are known to provide better shape overlay when screened against other biologic small molecules.
In addition, the molecules in the TCM database have been Hexagon shapes are genes. Dark red genes are of greater influence based on random forest analysis results (within top 5% of influential genes), white genes do not fall within the criteria of being in the top 5% of influential genes. All 4 common pathways enriched in patients receiving platinum and non-platinum therapy (described in Table 3) were merged into one pathway using Cytoscape. genes were clustered according to the Reactome Pathway Plugin available via Cytoscape. loss of ring structure is indicative of differences in influence of genes between patients receiving platinum and non-platinum therapies.

9
shown to be generally safe in people by merit of its long historical use in human populations. For CTNNB1, several ligands were identified via RCSB PDB (2s)-3-{[{[(2s)-2,3-dihydroxypropyl]aoxy}(hydroxy)phosphoryl]oxy}-2-[(6e)-hexadec-6-enoyloxy]propyl (8e)-octadec-8-enoate was the single ligand associated with CTNNB1 that was used for virtual screening against the FDA and TCM libraries. Unfortunately, neither library yielded a small molecule candidate with greater than a 50% Tanimoto score. The ligand of PLCG2 did yield several interesting drug leads in both the FDA-approved library and the TCM library. Fludarabine has been tested in oral cancer cell lines and found to be effective in inducing cell apoptosis. 77 Ganciclovir an HIV drug has also been tested in oral cancer and found to have a clinical effect on cell differentiation. 78 Entecovir and didanosine are drugs used in the treatment of HepB infection and HIV and have not yet been tested on oral cancer. Small molecules in the TCM database meeting shape overlay filters were [4-[2,6-dimethylmorpholin-4-yl)sulfonylphenyl]-[4-(2-phenoxyethyl)piperazin-1-yl]-methanone which has not yet been used on oral cancer cell lines. There was overlap between drug leads for ligands of UBB and PLCG2. This is due to the structural similarity between cytosine and guanine ligands used as reference molecules for similarity matching. Drug leads for UBB included Cytarabine (Cancer), Fludarabine (Cancer), Azacitidine (myelodysplastic syndrome), Gemcitabine (Cancer), and Lamivudine (HIV). Gemcitabine is unique in that it is the only drug of those listed, which has been approved by the FDA for use in oral cancer patients. For HDAC3, there were no matches exceeding a Tanimoto score threshold of 50% of the reference molecule when using a library of FDA-approved drugs. The TCM database did yield a match with an extract from Mallotus philippensis, a member of the Euphorbiacae plant family. An extract of this plant known as Rottlerin has been found to inhibit growth of colon cancer and breast cancer cells. 79,80 Quercetin and Diosmetin were other phenols found in citrus that were also identified as matches meeting Tanimoto thresholds. There are no reports of the effect of quercetin or diosmetin in oral cancer. It should be noted that Quercetin and Rottlerin have been noted in literature as promiscuous ligands that are often found in natural product in silico screenings. 81 A recent study reported quercetin as the number one natural product in terms of number of occurrences within the database. 82 The aforementioned studies do show that Rottlerin is related to the metastatic potential and viability of colorectal cancer cells. 79 Caution and validation of results with existing literature or carefully designed follow-up experiments should always be pursued to justify the results of promising in silico analyses. These analyses enhance knowledge of genes influencing treatment response in oral cancer. This study describes how pathway, network, and chemical informatics analysis can be paired with literature review to identify drug leads for oral cancer treatment. Reference ligands associated with influential genes in the Royal Chemistry Society Protein Database are listed along with drug leads and their corresponding Tanimoto similarity scores in Table 4.

Discussion
A machine learning approach known as random forest was used to identify genes influencing oral cancer treatment response specific to the platinum-based chemotherapy treatment type and the non-platinum-based chemotherapy treatment type. This article emphasizes the benefits of integrating the results of this line of analyses with pathway, network, and chemical informatics analysis to identify promising gene targets and drug leads. Biological plausibility of these findings was highlighted with a review of existing literature supporting the findings for pathways, genes, and small molecules that our reported approach identified as influential in oral cancer. The results of this work identify pathways influencing treatment response in platinum-based chemotherapy users, non-users, and those common to both users and non-users. Network analysis via Cytoscape allowed for the identification of those influential genes within each treatment modality group within the context of interconnected gene networks. The utility of random forest was underscored in that in addition to pathways it also provides a rank to each gene in its influence on treatment response. This approach is a low-cost method of prioritizing gene targets and drug leads. These methods are validated in that the genes identified have been shown to be associated with cancer progression in oral cancers and other cancers. Several drug leads identified were also shown to be effective in inhibiting oral cancer cells and were reported to be in different phases of the drug approval pipeline. A possible criticism of the method outlined in this study is that there is uncertainty in the degree of trust that should be extended to random forest measures of gene influence and the inference of importance to the pathways in which "influential" genes reside. To further such criticisms, a point could be made that the Tanimoto threshold of >50% similarity could be perceived as low and the 50% difference in the molecules compared may prevent activity and may also be shown to have toxicity for a given disease state. Given such uncertainties, it may seem that the evidence supporting these methods is tenuous. This study recognizes these criticisms; however, the counterpoints must be made that gene influence is not observed in a single sample of the data, but rather in more than 20 000 permuted samples of the data in which the top-ranked genes were found to be more influential than thousands of other genes. The computational intensity provided in this study (20 000 trees and 320 tries at each node of each decision tree) provides justification of the trust provided for each gene influence value. In respect to the results yielded by chemical informatics analysis, it is important to note that the outlined chemical informatics method was able to identify gemcitabine as a drug that has been approved for use in oral cancer by the FDA. This method also identified Fludarabine and Ganciclovir which have both been reported as providing significant reduction in oral cancer cell line progression and viability.
Drug leads were identified in both FDA and TCM libraries, and the benefit in expanding the number of libraries is that it increases the probability of finding a match meeting the Tanimoto threshold of >50%. The negative aspect of adding libraries is that if computational resources are not planned for accordingly, then the amount of time required to screen against each reference molecule will scale upward with library size. The tools used in this study were all open source and freely available, and a limitation to the adoption of this pipeline is that tools and their dependencies are distributed across different R repositories that may or may not be kept up to date. Combining these tools into a single package that allows for the identification of both gene targets and drug leads may enhance the pace of drug discovery pipelines. We have shown in this study that random forest is well suited to datasets with small observations and high number of features. Gene targets that have been shown (through literature review) to be associated with treatment response and cancer progression were identified through this study's use of random forest analysis. Stratifying this analysis by the type of chemotherapy received allows for interpretation of influential genes and pathways within the context of treatment. Indeed, the lack of overlap in the importance of genes from one treatment modality to another highlights the idea that the gene expression patterns influencing platinumbased treatment response differ from those gene expression patterns influencing non-platinum treatment response.
It is likely that there is bias inherent to the stratification of patients by chemotherapy treatment type. Chemotherapy treatment is associated with clinical variables like clinical stage, tumor size, and tumor grade, as well as gender and socioeconomic quintile. [83][84][85] This study attempted to address these confounders by including them within the bag of randomly selected features available for construction of the random forest model. By integrating chemical informatics analyses, random forest results can be translated into lists of drug leads for each target gene. This method identified drug leads that have already entered or passed phase 3 trials. Our review of identified drug leads and comparison with existing annotations show that the chemical informatics methods described can identify small molecules with therapeutic potential. This study provides the impetus for further exploration of the role of the identified small molecules in oral cancer treatment response and the targeting of those genes identified as most influential by our series of analyses. This study also serves as a model for researchers identifying gene targets in rare cancers where the number of cases is limited.

Author Contributions
JS was the sole contributor to manuscript composition, experimental design, analysis, and production of all tables and figures in this manuscript.

Supplemental Material
Supplemental material for this article is available online.