**3. Results**

We collected mRNA sequencing, exome sequencing, and whole genome sequencing data from previously published or publicly available resources as depicted in Table 3. We applied a combination of methods to identify biological mechanisms of interest in the samples that had patient-derived transcriptomic data. Briefly, we employed a transfer learning-inspired approach to group transcripts into latent variables (LVs), and then selected the LVs that best separate out tumors by tumor type using an ensemble of random forest models. We evaluated the selected LVs for patterns of immune cell gene expression, protein activity, and gene variants.

**Table 3.** Summary of individuals and samples for the NF1 nerve sheath tumors used in this study. All samples have gene expression data and a subset have genomic data derived from whole-exome sequencing or whole genome sequencing. Some neurofibromas did not have more specific pathologic subtyping information available, and therefore were classified as "undefined neurofibromas" or NFs.


#### *3.1. Pan-NF Transcriptomic Analysis Identified Most Variable Latent Variables in NF1*

To measure the diversity of the nerve sheath tumors at the transcriptomic level, we re-processed RNA-seq data from three published datasets [37–39] and one unpublished dataset. Despite having four types of nerve sheath tumors (cNFs, pNFs, NFs, and MPNSTs) across four datasets, we observed confounding batch e ffects (Figure 1A) as some tumor types (e.g., cNF and NF) were derived from separate studies. These batch e ffects, together with a lack of transcriptomic data from normal tissue from NF1 patients, precluded us from carrying out any meaningful di fferential gene expression analysis and motivated us to seek alternate approaches to pathway identification.

**Figure 1.** Transfer learning reduced dimensionality and added additional context to gene expression datasets. ( **A**) Principal components analysis (PCA) of gene expression data indicated that counts-level data may have been batch confounded. (**B**) The relative distributions of latent variable expression across the four tumor types using a density plot indicated that the majority of latent variables (LVs) had an expression value near 0 and that the four tumor types had similar latent variable expression distributions. ( **C**) PCA of LVs indicated that batch e ffects, although reduced, may still have existed in the LV data ( **D**) A look at the 5% most variable LVs across the cohort of gene expression data indicated that the latent variables represented a wide swath of biological processes, as well as some LVs that had no clear association to a defined biological pathway.

To reduce the dimensionality of the RNA-seq data and increase the biological interpretability, we applied MultiPLIER, a transfer learning approach trained on the recount2 dataset, an independent dataset comprising thousands of gene expression experiments. This analysis quantified the expression of 962 latent variables (LVs) in all 77 samples using their gene expression data (Figure 1B, Figure S1, Table S1). We anticipated that as long as technical confounders in our dataset were independent of the technical confounders in the recount2 dataset, we would be able to find some latent variables that were independent of batch effects. We observed that expression of the latent variables was uniformly distributed across the four tumor types (Figure 1B). Furthermore, principal components analysis suggested that this method reduced the size of the previously observed batch effects (Figure 1A,C) as the within-cluster distances were significantly reduced (Figure S3, *p*-value = 2.1 × <sup>10</sup>−182, Wilcoxon test).

We then sought to identify which LVs characterize the differences between individual tumors. To find these, we examined the LVs with the 5% largest standard deviation across all LVs. Many of the resulting latent variables (Figure 1D) were found to be associated with known biological mechanisms including metabolic, immune, and transcription-related signatures. However, this analysis did not explore the association between tumor types and latent variables. Therefore, we performed additional analyses to (1) identify the gene expression features that best define each tumor type and (2) elucidate the biological mechanisms underlying the most important latent variables.

#### *3.2. Ensemble of Random Forests Identified Latent Variables That Robustly Describe Individual Tumor Types*

To identify these expression features, we used random forest models, an approach that is known for its ability to characterize complex decision spaces better than basic clustering approaches [28]. Random forests achieve this by building decision trees that account for both the values of the features and conditional logic at each "branch" of the tree. Although random forest models are generally used to classify samples, in this study we leveraged the capacity of random forest models to identify features in the data that provide meaningful information. Therefore, in this study, our focus was to identify important features rather than building the most accurate classifier.

We generated a random forest using NF-specific latent variable scores (Table S1) as described in the Materials and Methods section and depicted in Figure S2. Given the limited number of samples in our dataset, particularly for some tumor types, we anticipated high variability in the performance of our supervised learning model. To obtain a distribution of model performance, we generated 500 training sets (using stratified sampling from the full dataset without replacement), each training their own random forest. We used the distribution of model accuracy measurements and feature importance scores from these forests to estimate overall model performance and feature importance (Figure 2A). We then compared the class-specific accuracy measurements (median accuracy scores: cNF = 0.97, MPNST = 0.67, NF = 0.67, pNF = 0.57) (Figure 2A) to those from forests built with only the top 40 LVs from each class (a total of 98 LVs) from the first forest ensemble. Using the top 40 set of LVs, we were able to improve model performance for each tumor type on an independent test set (median accuracy scores: cNF = 1.00, MPNST = 0.86, NF = 0.86, pNF = 0.75) (Figure 2B). For each tumor type, reducing the feature set to these 98 LVs significantly improved the median accuracy scores for the models (Mood's median test *p*-values: cNF < 2.2 × <sup>10</sup>−16, MPNST = 3.789 × <sup>10</sup>−11, NF = 7.297 × <sup>10</sup>−09, pNF < 2.2 × <sup>10</sup>−16). This suggested that the 98 selected LVs with high importance scores for each tumor type were sufficient for characterizing the specific tumor types (Figure 2C).

**Figure 2.** An ensemble of random forests selected the most important latent variables for classifying different tumor types in NF1. (**A**) Density plot showing the distribution of F1 scores of 500 iterations of independent random forest models using all latent variables. (**B**) Density plot showing the distribution of F1 scores of 500 iterations of independent random forest models trained using only the top 40 features with high importance scores for each class obtained from models included in (**A**). (**C**) Ridgeplots of top 20 latent variables selected by the random forest for each tumor type and their importance scores for each class that were selected for later analyses.

#### *3.3. Selected Latent Variables Represented Distinct Biology of Nerve Sheath Tumor Types*

To further probe the biology underpinning the distinction between NF1 tumor types, we focused on the 98 latent variables selected by the ensemble of random forests, depicted in Figure 3 and listed in Table S2. We evaluated the tumor-wise LV expression as well as the contributing genes (loadings) for each LV. For example, some latent variables that were selected by the random forest as relevant to predict all four tumor types (Figure 3A) showed differences in expression between various tumor types (dot plots in Figure 3Biii,Ciii), whereas others were less distinct (Figure 3Bi,Ci). By investigating the loadings of each LV, we tried to map them to known biological pathways. For each LV, the loadings of the constituent genes denoted the contribution of that gene to the particular LV. These gene lists associated with the selected LVs can be used to find testable candidates for each tumor type

for downstream analysis. For example, one of the four latent variables predictive of all four tumor types was enriched in genes associated with neuronal signaling (Figure 3Bii), suggesting that this measurement is required for the model to distinguish various tumors as the presence of neuronal tissue likely varies across tumor types. Other functionally enriched latent variables, such as those depicted in Figure 3Bi,iii, implicate other biological pathways in tumor growth. However, many of the selected latent variables still remain uncharacterized (few examples shown in Figure 3Ci–iii).

**Figure 3.** Selected latent variables (LVs) represented gene combinations unique to each tumor type. (**A**) Venn diagram showing the distribution of the top 40 LVs from each tumor type. (**B**,**C**) Total values of the LVs as measured by multiPLIER across samples are represented in the dot-plots, where color of the dots represents the tumor type ("Class" label colors described in the lower left). Loading values for the top 10 genes for each LV are represented in bar-plots below. The higher the loading, the greater impact that the gene expression had on the total multiPLIER value. (**B**,**i**–**iii**) Genes constituting the latent variables associated with known cell signaling pathways. (**C**,**i**–**iii**) Genes constituting the uncharacterized latent variables.

#### *3.4. LV Scores May Be Attributed to Specific Gene Variants for Specific Tumor Types*

To characterize the 73 LVs with no associated pathway information, we focused on the 40 samples from our dataset (Table 3) that had matched gene variant data (WGS or exome-Seq) to assess if there were any genes that, when mutated, caused a significant change in LV expression. The results of this calculation are found in Table S2.

We identified 22 latent variables that were significantly (Benjamini–Hochberg adjusted *p* < 0.01) associated with single gene variants (Table S3). These latent variables, along with the genes whose variants are associated with altered expression, are depicted in Figure 4A. This approach failed to identify any genes that were mutated across multiple tumor types. In the list of genes in which variants were associated with changes in latent variable expression, we identified nine genes with variants in cNFs and two genes with variants in one neurofibroma sample. An example of how these variants are associated with the expression of a latent variable is depicted in Figure 4B. Mutations in the nine cNF-variant genes are associated with lower expression of LV 851. Because most of these variants occur in cNF samples but not in the other tumor types, it is not surprising that this latent variable is down-regulated across all cNFs (Figure 4C). Figure 4D shows that the MAP1B gene, which has been reported to play a role in glioblastoma [51] but has not been studied in NF1-linked tumors, is a major contributor towards this LV.

**Figure 4.** Some genes significantly distinguished expression of latent variables. (**A**) Latent variables (*y*-axis) whose values are significantly altered by mutations in specific genes. (**B**) MultiPLIER value of LV 851 across tumor samples. (**C**) MultiPLIER value of LV 851 across all samples. (**D**) Loading values of the top 20 genes that comprise LV 851.

#### *3.5. Selected Latent Variables Represented Specific Immune Cell Types in the Tumor Microenvironment*

Given the limited representation of different NF1 tumor types in our genomic variant dataset, we used alternate gene expression metrics to assess the biological underpinnings of the 73 uncharacterized latent variables. We performed tumor immune cell deconvolution analysis [29–31] to identify potential immune infiltration signatures present in the individual tumors using CIBERSORT and MCP-counter (Table S4). Specifically, CIBERSORT deconvolution indicated the presence of activated mast cells and M2 macrophages in all tumor types (Figure 5A). The analysis further suggested that all of the tumors have a population of resting CD4+ memory T cells. The results from MCP-counter, depicted in Figure 5B, show a complementary view of the tumor cell types due to the slightly different categorization of cell types. Specifically, we found the presence of cancer-associated fibroblasts across all tumors that were not captured in CIBERSORT.

We then used the immune scores from CIBERSORT and MCP-counter to probe some of the latent variables selected by the random forest model to better understand their role in NF1 tumor biology. We first searched for tumor deconvolution scores that were correlated with latent variable expression across all tumors. Figure 5C,E show specific LVs where the immune scores of the samples highly correlated with a predicted immune cell type. LV 546, for example (Figure 5C), was found to have a high correlation of activated mast cells in cNF samples. Figure 5E shows that LV 540 has immune scores strongly correlated with T cells in all NF1 tumor types. Complete results of immune scores across all tumors can be found in Table S4. Together, these results sugges<sup>t</sup> that the selected LVs capture signatures from the tumor microenvironment in NF1 samples.

**Figure 5.** Various immune cell signatures correlated to specific LVs that differentiate tumor types in NF1. (**A**) CIBERSORT deconvolution of bulk nerve sheath tumor expression data predicted the presence of activated mast cells and M2 macrophages and resting CD4+ memory T cells in all of the tested tumor types. (**B**) MCP-counter based deconvolution of bulk nerve sheath tumor expression data predicted the presence of cancer-associated fibroblasts across all tumor types, and diversity in T cell population across tumor types. (**C**) Correlation of CIBERSORT immune score (*x*-axis) with expression of latent variable 546 highlighted the increased presence of activated mast cells and resting dendritic cells in cNFs (circles). (**D**) Top 20 gene loadings of LV 546. (**E**) Correlation of MCP-counter score of Tell infiltration (*x*-axis) with LV 540. (**F**) Top 20 gene loadings of LV 540.

#### *3.6. Selected Latent Variables Captured Protein Regulatory Networks in NF1 Tumors*

Of the 98 latent variables selected by the random forest, 55 could not be characterized via enrichment of the gene loadings (Figure 3), mutational patterns (Figure 4), or immune subtypes (Figure 5). To characterize them, we applied the metaVIPER algorithm [32] to identify specific protein activity measurements that correlated with latent variable expression. This algorithm leverages

previously published regulatory network information [52] to infer protein activity in each sample to assign numerical scores of activity across 6168 proteins for each of the 77 samples (Table S5). We then measured the correlation of these scores and latent variable scores (Table S6, Figure 6A) to identify functional aspects of the latent variables.

**Figure 6.** Integration of protein activity information with LVs can identify candidate drug targets for different NF1 tumor types. (**A**) A heatmap of correlation scores of known proteins with regulatory networks (or regulons) that are represented in the characterized and uncharacterized LVs selected above. The green bar across the top depicts how many protein activity scores had a Spearman correlation greater than 0.65. (**B**) Clustering of the LV-correlated VIPER proteins highlighted five clusters of latent variables with similar VIPER protein predictions, suggesting that these five clusters may have functional overlap. (**C**) Mean LV expression within the clusters highlighted differential expression within the clusters across tumor types. Tumor type is indicated by colors on the right. (**D**) Drug set enrichment analysis of the average VIPER protein correlation of cluster 2 identified some drugs and preclinical molecules that are enriched with targets in this cluster.

As seen in previous work [43], we observed that multiple latent variables exhibit similar correlation patterns with active proteins (Figure 6A). We clustered the correlation scores to see if we could group the latent variables with similar protein activity predictions (Figure 6B). In clustering the latent variable–protein activity scores, we observed five distinct clusters of latent variables with similar predicted protein activities (Figure 6B). Aggregation of each cluster of latent variables (mean expression within each cluster) demonstrated that these functional clusters were di fferentially expressed in the di fferent NF1 tumor types (Figure 6C).

We then assessed the druggability of each of these five clusters by taking the average correlation of each protein within the cluster and performing gene set enrichment analysis against a database of small molecules with known biological targets [49]. This enabled the identification of drugs and drug-like compounds that are significantly enriched for targets in each cluster (Table S7). For example, cluster 2, which is expressed in pNF, NF, and MPNST more than it is expressed in cNF, has correlated VIPER proteins that are enriched for both clinically approved drugs (dovitinib) and drug-like small molecules (Figure 6D). Furthermore, we found that cluster 3 was enriched for compounds that a ffect cell cycle progression (e.g., dinaciclib, abemaciclib), whereas clusters 1 and 5 were enriched for compounds such as CUDC-101 (7-[4-(3-ethynylanilino)-7-methoxyquinazolin-6-yl]oxy-N-hydroxyheptanamide) and analogs that inhibit histone deacetylases (Figure S4).
