Hydra: A mixture modeling framework for subtyping pediatric cancer cohorts using multimodal gene expression signatures.

Precision oncology has primarily relied on coding mutations as biomarkers of response to therapies. While transcriptome analysis can provide valuable information, incorporation into workflows has been difficult. For example, the relative rather than absolute gene expression level needs to be considered, requiring differential expression analysis across samples. However, expression programs related to the cell-of-origin and tumor microenvironment effects confound the search for cancer-specific expression changes. To address these challenges, we developed an unsupervised clustering approach for discovering differential pathway expression within cancer cohorts using gene expression measurements. The hydra approach uses a Dirichlet process mixture model to automatically detect multimodally distributed genes and expression signatures without the need for matched normal tissue. We demonstrate that the hydra approach is more sensitive than widely-used gene set enrichment approaches for detecting multimodal expression signatures. Application of the hydra analysis framework to small blue round cell tumors (including rhabdomyosarcoma, synovial sarcoma, neuroblastoma, Ewing sarcoma, and osteosarcoma) identified expression signatures associated with changes in the tumor microenvironment. The hydra approach also identified an association between ATRX deletions and elevated immune marker expression in high-risk neuroblastoma. Notably, hydra analysis of all small blue round cell tumors revealed similar subtypes, characterized by changes to infiltrating immune and stromal expression signatures.

Introduction

Large cancer sequencing projects, including The Cancer Genome Atlas (TCGA) and Therapeutically Applicable Research to Generate Effective Treatments (TARGET), have facilitated the development of cancer gene expression compendia [1-6], but these compendia often lack expression data from corresponding normal tissue. Without the normal comparator, Hoadley et al. (2018) found that cell-of-origin signals drive integrative clustering of TCGA data. Strong cell-of-origin and tumor microenvironment (TME) signals may also complicate the interpretation of gene expression results for precision oncology applications, so careful modeling of the data is necessary to infer accurate conclusions.

The TME includes tumor cells, stromal fibroblasts, immune cells, and vasculature [7]. Similarities in TME composition across tumor samples have led to the identification of TME states (e.g. inflamed, immune-excluded, immune-desert). While these states are dynamic, they can still shed light on the immunogenicity of tumor cells and correlate with response to cancer immunotherapies [8]. The TME cellular composition can be inferred from tumor RNA-Seq data since host cell RNA is sequenced along with the cancer cell RNA. Tumor progression and response to therapies is associated with features of the TME. Therefore, targeting the TME therapeutically may improve treatment outcomes in some cancers.

Immunotherapies that activate the host immune system to eradicate tumors have been effective in treating several cancer types, particularly cancers with a high mutation burden [9, 10]. Pediatric cancers tend to have fewer mutations than adult cancers, and while there has been limited testing of immunotherapies in pediatric cancer patients, the currently available data suggest lower response rates than adult cancers [11, 12]. However, improved immune subtyping of pediatric cancers may identify subsets of patients that are candidates for powerful immunotherapies. In addition to infiltrating immune cells, cancer-associated fibroblasts (CAFs) assist in extracellular matrix remodeling and activation of growth factor signaling. CAFs facilitate tumor growth, metastasis, and resistance to some therapies, so identification of CAF functions within a tumor may also facilitate clinical decision making. Methods are needed to both infer and characterize gene expression subtypes that correlate with tumor microenvironment states to accelerate the development of personalized therapies for pediatric cancers.

Tumor/normal differential expression analysis in which a cohort of tumor tissues is compared to corresponding normal tissue samples is an effective approach for identifying gene expression biomarkers [13-15], but it is often not possible to conduct this analysis in a clinical setting. Sufficient biological and technical replicates are limited by tumor tissue availability, and healthy neighboring tissue often cannot be isolated. In addition, for many pediatric cancers, the cell-of-origin, and thus the appropriate reference normal tissue, is not known. Besides differential expression analysis, single-sample pathway analysis can be used to identify upregulation of biological gene sets in tumor subtypes. Among the most widely used pathway analysis approaches is gene set enrichment analysis (GSEA) [16, 17]. GSEA identifies coordinated expression of pathway genes using gene ranks and a Kolmogorov-Smirnov-like test statistic. GSEA is usually performed on differentially expressed genes to compare two cohorts or phenotypes, but single-sample GSEA is also available when there is not an obvious comparator. GSEA uses curated pathway gene sets like those in the Molecular Signatures Database (MSigDB) [18].

Cancer gene expression subtypes are traditionally identified using unsupervised clustering methods such as consensus clustering analysis [19-21]. These methods are generally underpowered because the number of genes greatly exceeds the number of samples. Dimensionality reduction approaches such as Principal Component Analysis (PCA) have been found to underestimate the dimensionality of gene expression data [22]. Lenz at al. (2016) found two cases in which PCA fails to identify a biological signal: when the size of the cluster is small and when the effect size is small. Lenz et al. (2016) suggests investigating multimodally expressed genes to improve identification of cancer subtypes. Cancer subtypes naturally lead to multimodal expression patterns because each subtype expresses a correlated set of genes at different expression levels. Expression subtypes may result from dysregulated pathway expression within cancer cells, but another source of multimodal expression comes from varying amounts of infiltrating immune and stromal cells in the TME.

Gaussian mixture models are a powerful class of unsupervised clustering algorithms that can be used to detect multimodally expressed genes [23-25]. A Gaussian mixture model is appropriate when the expression data can be modeled as a mixture of two or more Gaussian distributions [26]. One limitation of Gaussian mixture models in this context is that the number of clusters in the data is often not known beforehand, so a parameter search must be used to identify the best-performing model. However, this is a computationally expensive approach. This problem can be overcome by placing a Dirichlet process prior on the number of expression clusters. The number of clusters is then inferred while fitting the mixture model using Markov chain Monte Carlo (MCMC) sampling [26]. This approach has not been widely used in clinical cancer research because these algorithms are still computationally expensive, but recent advances in Bayesian variational inference have made this approach scalable for precision oncology applications [27].

Here, we present the hydra framework for identifying clinically relevant expression subtypes and classifying N-of-1 tumor samples using learned models. We provide an overview of the hydra framework, assess performance for detecting differential pathway expression, and apply the framework to better understand expression patterns in high-risk neuroblastoma and other small blue round cell tumors. We apply the learned models trained on publicly available cancer gene expression data to the N-of-1 setting and show that this framework can identify distinct immune and stromal expression signatures that differentiate pediatric cancer samples. Finally, we identify recurrent tumor microenvironment signatures across pediatric cancer types associated with differences in patient outcomes.

Materials and methods

Dirichlet process gaussian mixture model

Traditional parametric models, like the finite mixture model, use a fixed number of parameters (i.e. number of clusters). Over- or underfitting can occur when the parametric model does not reflect the underlying data [48]. Unlike the finite mixture model, the Dirichlet process mixture model (DPMM) represents a theoretically infinite number of clusters and can adapt the number of clusters based on prior belief and the data [26, 48, 49].

The Dirichlet process (DP) is an infinite dimensional extension of the Dirichlet distribution [50] and is commonly used as a prior distribution for infinite mixture models [51, 52]. The Dirichlet process has two parameters: the concentration parameter α and centering distribution H. The concentration parameter α, where , controls the extent to which samples from the DP resemble the centering distribution H. We model gene expression as a multivariate Gaussian distribution, so our centering distribution is a normal-Wishart distribution ().

We briefly describe the stick-breaking construction of the Dirichlet process G ∼ DP(α, H). Consider a stick of unit length. To generate an infinite number of mixing weights π1, π2, …, π for the DPMM, first break a stick of unit length at ν ∈ [0, 1] where ν is sampled from a Beta distribution, and set π1 to be the length of the first piece. We repeat this process using the remainder of the stick for each π. The DP is truncated to the number of clusters K, which was shown to accurately approximate the infinite posterior for large K [26, 48, 50, 53–55].

Next, we sample the parameters from the centering distribution H weighted by the mixing components. If we consider a probability space Θ where , then H is a measure on the partitions of Θ. For our application, we will partition the parameter space Θ into finite, measurable partitions B1, B2, …, B.

This construction generates the marginal of the Dirichlet process, which follows a Dirichlet distribution. Samples from the marginal distribution are finite, discrete, and sum to 1 [50]. Next, we outline how the DPMM groups gene expression samples x under cluster-specific parameters and where z ∈ 1, 2, …, K is the cluster index.

To improve our methods ability to scale to larger datasets, we incorporated the bnpy memoized online variational inference algorithm (moVB) [53] into our analysis framework. The moVB algorithm uses variational inference to approximate the posterior distribution and interleaves birth, merge, and delete moves to avoid local optima and remove redundant clusters [56]. We found that the moVB algorithm accurately identified the number of clustering on validation datasets (S1 Fig), whereas standard MCMC sampling procedures tended to overestimate the number of clusters.

Hydra method

We developed a Bayesian non-parametric clustering framework for identifying biological and technical variation in large cancer gene expression datasets without the need for a reference normal dataset. To our knowledge, this is the first reproducible and widely deployable implementation of a non-parametric mixture model framework designed to overcome the challenges of precision oncology gene expression analysis. The hydra pipeline is an open source software tool hosted on GitHub (www.github.com/jpfeil/hydra). A Docker container is available for deployment across environments (https://hub.docker.com/r/jpfeil/hydra).

The hydra framework contains three main command-line tools: filter, enrich, and sweep (Fig 1). The filter command is run first to isolate the multimodally expressed genes using a univariate Dirichlet Process Gaussian Mixture Model (DP-GMM). There are two methods for analyzing the resulting set of multimodally expressed genes. The enrich method, which subsets to the genes found to be significantly enriched in biological pathways, and the sweep method, which searches within user-defined gene sets for multimodal expression signatures. The underlying analysis routines can be accessed within the Docker using Jupyter notebooks to facilitate the development of user-defined workflows.

Fig 1

Overview of the hydra framework tools.

A: Suggested workflow for applying hydra framework tools to identify clinically relevant gene expression subtypes. B: The hydra filter command removes unimodally distributed genes which greatly reduces the number of genes in downstream clustering analysis. C: The hydra enrich command takes the multimodally expressed genes and returns enriched gene sets. The enriched gene set genes are used for multivariate clustering of samples. D: The hydra sweep command looks for multivariate normal clusters within user-defined gene sets. This can be used for the automatic detection of clusters in large gene set databases. Abbreviations: Tumor microenvironment (TME).

The filter command (Fig 1B) takes an expression matrix and filters the genes down to the multimodally expressed genes using the DP-GMM described above. We apply a DP-GMM to each gene, saving the model for genes with two or more expression clusters. This creates a directory of multimodally expressed gene models which can be used to predict differential expression in new samples. This analysis framework is a novel contribution to the precision medicine research community. Our approach has several beneficial properties. For example, training models on curated data sets and applying the models to new samples avoids the use of reference distributions, which overestimate the uncertainty in the analysis by not accounting for subtype expression. Furthermore, this approach identifies the set of most strongly differentially expressed genes within a disease context, which may enrich for potential biomarkers for precision medicine applications. The multimodally expressed genes are also used in downstream clustering analysis.

The enrich (Fig 1C) and sweep (Fig 1D) routines are two independent analyses to explore multimodal expression in cancer gene expression cohorts. In addition to identifying expression variation within a disease context, we also found that multimodally expressed genes that participate in a biological pathway tend to have correlated expression distributions. This insight facilitates the detection of multimodal expression signatures by enriching for genes that have multimodal expression distributions and participate in known biological processes. The hydra software comes prepackaged with popular gene sets, including the Molecular Signatures Database (MSigDB) [18], the Gene Ontology terms [57, 58], and the EnrichmentMap gene sets [59]. The gene set database is configurable, so additional gene sets can be added at runtime.

The enrich command uses a hypergeometric test [60] to discover enrichment of multimodally expressed genes within a user-defined database of gene sets. This creates a list of gene sets and a list of enriched gene set genes. The enrich method outputs a table of enriched gene sets while also clustering samples across the genes that participate in the enriched gene sets. The table of enriched gene sets may reveal surprising expression patterns and generate hypotheses for further investigation of tumor subtypes.

The implementation of the enrich method includes an important parameter known as the minimum component probability. The minimum component probability is the probability of placing a sample within the smallest expression cluster. This is an additional filter to remove multimodally expressed genes that influence a relatively small subset of tumor samples. This parameter gives the user the ability to subset the enriched genes to those that influence a greater number of patients. To aid in the exploration of minimum component thresholds, we implemented a scan sub-routine. The scan routine tunes the analysis with respect to the constraints of the available data (e.g. number of samples and number of genes), which is an important factor in pediatric cancer research since data is often difficult to obtain and so datasets are relatively small. We recommend setting this threshold such that the number of genes is less than the number of samples because otherwise the inference may become unstable [61].

The sweep routine identifies differentially expressed gene sets and can be used as an alternative to single-sample GSEA [16]. For each gene set, a multivariate DP-GMM is applied to determine if more than one expression cluster is present within the gene set. This approach is useful when curated gene sets are available for the disease of interest, but manual inspection of each gene set is not feasible. Reducing the genes to multimodally expressed genes facilitates the detection of differentially expressed gene sets. Existing gene set enrichment tools are known to under-perform when the expression is correlated [62], but our approach is designed to identify distinct correlation structures within gene expression datasets.

We have also implemented routines for cluster profiling and N-of-1 tumor analysis. These routines are accessible within the docker container using the Jupyter notebook command. Cluster profiling analysis of clusters derived from the enrich or sweep routines includes GSEA [63] to identify the pathway expression that characterizes each cluster. GSEA uses all available genes since it requires non-differentially expressed genes to assess the significance of an enrichment score. A t-statistic is calculated for each gene, comparing gene expression values of samples inside to those outside of a cluster. Cluster profiling GSEA uses the ranked gene-level t-statistics to determine gene set enrichment.

The N-of-1 tumor analysis routine classifies a new gene expression profile into one of the inferred clusters, calculates a gene-level z-score for that sample relative to the normalized expression distribution, and performs standard GSEA using a preranked list of z-score values [63]. This procedure can identify new gene expression signatures that may not be detectable using the entire expression cohort as a background reference distribution. This approach is another novel contribution to the field and may facilitate the identification of clinically relevant signatures that are being overlooked in current gene expression analyses.

Synthetic data generation and validation

We first tested the hydra framework’s ability to detect differential pathway expression using synthetic cancer data. We compared hydra sweep to two widely used gene set enrichment tools: single-sample gene set enrichment analysis (ssGSEA) and gene set variation analysis (GSVA) [64-66]. Both methods are implemented in the GSVA R package [65]. In order to accurately model correlation structures within cancer cohorts, we modeled the synthetic cancer gene expression data as a multivariate Gaussian distribution. We used the TCGA glioblastoma multiforme (GBM) cohort (N = 166) to model a background mean and covariance matrix for the synthetic data analysis. We chose TCGA GBM, a very different disease from those analyzed in the remainder of this manuscript, to avoid overfitting the hydra method to diseases of interest. This also enables us to demonstrate the flexibility of our method to analyze data from a variety of cancer genome sequencing projects.

This approach allowed us to model cancer gene expression data while also controlling for subtype-related expression variation. We downloaded the RSEM-quantified transcripts per million (TPM) normalized gene expression measurements from the UCSC Xena Browser [3]. We focus our analysis on normalized gene expression data because this data is more widely used in the cancer research community and fewer methods are available to analyze normalized counts. To reduce heteroscedasticity and the effect of outlier expression levels, we transformed the expression data to log2(TPM + 1) [67].

We defined an expression subtype as a subset of samples with a distinct expression mean and correlation structure compared to other samples within the disease cohort. To avoid biases in the synthetic data generation process, we used random sampling to select MSigDB gene sets for each subtype, the size of the subtype, and the correlation structure within the subtype. We randomly generated a covariance matrix for the cancer subtype expression data, but used the underlying covariance matrix of the TCGA glioblastoma multiforme dataset for the background samples. We tested the effect of having 10% and 25% of genes within a gene set being differentially expressed (%DEG). In addition to these parameters, we tested a range of effect sizes: 0.25 (least different), 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, and 3.0 (most different). This process was repeated twice for each gene set to create synthetic training and test data, which resulted in the generation of 640 synthetic datasets.

We then applied the hydra framework using the hydra sweep command (Fig 1C), since this method is directly comparable to the single-sample GSEA methods. The mean expression filter removed any genes with a mean expression of fewer than 1.0 log2(TPM + 1). This avoids lowly-expressed genes that may have particularly noisy expression measurements. The prior on the hydra covariance matrix was the identity scaled by 2.0 and the prior on the number of clusters was set to 2 because we expect there to be an activated cluster and a baseline expression cluster. We set the over-expressing cluster to be the cluster with the largest L1 norm.

Pediatric cancer gene expression data

We downloaded pediatric cancer RNA-Seq data for neuroblastoma, osteosarcoma, Ewing sarcoma, alveolar rhabdomyosarcoma, and embryonal rhabdomyosarcoma from the UCSC Treehouse Compendium (https://treehousegenomics.soe.ucsc.edu/public-data/). This data was produced using the same RNA-seq pipeline, so potential computational batch effects are minimized [1, 6]. Clinical data for the TARGET neuroblastoma and osteosarcoma samples were obtained from the TARGET Data Matrix (https://ocg.cancer.gov/programs/target/data-matrix). We also analyzed a set of 58 synovial sarcoma microarray profiles with matching metastasis rate data [68].

TARGET neuroblastoma analysis

We applied each hydra tool to the TARGET MYCN-non-amplified neuroblastoma cohort. We first obtained the multimodal gene models using the hydra filter tool. The hydra filter tool identified all genes with a multimodal expression pattern. We used the mean expression filter to remove genes that may have unstable measurements due to low transcript abundances. We excluded all genes with a mean expression value less than 1 log2(TPM + 1).

The hydra sweep command was applied to search for subtype expression within curated MSigDB gene sets. We included the hallmark (n = 50), BioCarta (n = 289), KEGG (n = 186), PID (n = 196), and Reactome (n = 1499) genesets [18]. We include all signatures with a minimum component probability of 10%. For example, the smallest subtype cluster considered in this analysis had 7 samples, since the total number of samples was 70. We investigated relationships among differentially expressed gene sets by clustering the gene sets by their pairwise Jaccard index. This created a similarity network that was then visualized using the Gephi software tool [69].

The hydra enrich command identified correlated expression signatures using the enriched GO term genes (FDR < 0.01). The multivariate mixture model α concentration parameter was set to 5.0; the prior on the covariance matrix was set to the identity scaled by 2.0. The prior parameter for the number of clusters was set to 5. Our synthetic data analysis found that the signal decreases below an effect size of 1.0, so we use this parameter value for all following analyses. We used the hydra scan routine to search a range of minimum component probability thresholds (see Results) and found that a threshold/probability of 20% yielded the most clusters while keeping the number of genes (p = 42) below the number of samples (n = 70).

To validate tumor microenvironment expression subtypes, we correlated the hydra enrich expression clusters with the results of tumor microenvironment profiling tools xCell [31], CIBERSORT [70], and ESTIMATE [32]. We also compared the hydra enrich approach to state-of-the-art consensus clustering methods M3C [20] and k-means clustering using the Gap statistic to select the number of clusters [36]. Since these methods are influenced by the number of input genes, we tested a range of median absolute deviation (MAD) thresholds. The number of clusters was assumed to be the smallest statistically significant value.

Small blue round cell tumor analysis

We then compared the clustering patterns across MYCN-NA neuroblastoma, osteosarcoma, Ewing sarcoma, embryonal rhabdomyosarcoma, alveolar rhabdomyosarcoma, and synovial sarcoma. We applied the TumorMap dimensionality reduction method [5] to visualize clustering of the full small blue round cell tumor gene expression matrix. We then applied the hydra framework to explore expression variation within each disease. Each disease expression matrix had unique statistical properties including sample size and subtype variation. This required us to adapt the minimum probability threshold for each disease dataset using the scan routine. The Jupyter notebooks for exploring these datasets can be found on GitHub (www.github.com/jpfeil/hydra-paper/analysis). We used agglomerative clustering to investigate patterns in the top 10 enriched gene sets for each disease’s expression subtypes.

Statistical analysis

A Kruskal-Wallis test was used to identify statistically significant differences across two or more groups, and a Mann-Whitney U test was used for pairwise tests using a Holm-Sidak correction for multiple hypothesis testing [71, 72]. We used the scipy [73] stats implementation of the Kruskal-Wallis test and the scikit-learn post hoc processing [74] implementation of pairwise Mann-Whitney U tests. Spearman rank and Pearson correlation values were calculated using the scipy library [72]. Survival analysis was done using the survminer package [75].

H&E slide preparation and pathologist review

Pediatric tumor samples were flash frozen, embedded in OCT, and 5μm cryosections were collected. Slides were hematoxylin and eosin (H&E) stained and imaged on a Leica DMi8, equipped with a HC PL APO 40x/0.85 NA objective and DFC7000T camera. H&E slides were reviewed by a licensed pathologist. Morphologic analysis was performed and the degree and type of inflammation estimated from the histologic sections. Grading of inflammation was either minimal (<10% of total nuclei consist of inflammatory cells) or moderate (20-30% of total nuclei consist of inflammatory cells). The type of inflammation (predominantly small mature lymphocytes or mixed inflammation consisting of small mature lymphocytes along with plasma cells and/or eosinophils) was noted for each tumor sample.

Results

Performance assessment using synthetic gene expression data

To assess how well hydra detects differentially expressed pathways as compared to common pathway enrichment approaches, we applied the hydra framework to synthetically-generated cancer gene expression data. We generated synthetic cancer gene expression data based on the TCGA glioblastoma multiforme and the MSigDB Hallmark gene sets as described above. We tested a range of effect sizes and percent differentially expressed genes (%DEG) within the MSigDB gene sets. We generated receiver operator curves (ROC) and calculated the area under the receiver operator curve (AUC) for each analysis. Overall, the hydra pipeline outperformed the single-sample GSEA approaches with a mean AUC of 0.93 (95% CI: 0.91—0.95). ssGSEA had a mean AUC of 0.72 (95% CI: 0.71—0.74) and GSVA had a mean AUC of 0.67 (95% CI: 0.66—0.68) (Fig 2A).

Fig 2

Hydra sweep is more sensitive than existing gene set enrichment approaches for detecting differential pathway expression in synthetic data and scales well to large datasets.

A: Mean receiver operator curves across effect sizes, percent differentially expressed genes (%DEG), and MSigDB Hallmark gene sets. A larger area under the curve (AUC) indicates better performance. The average AUC and 95% confidence interval for each method are in the ROC plot figure legends. B: Line plots comparing the mean AUC across a range of effect sizes and %DEG values. C: Box plot showing mean runtimes for differential pathway analysis where the effect size is fixed but the sample size varies. D: Line plot comparing the mean runtimes for differential pathway analysis across a range of sample sizes.

We further investigated the performance of these methods by plotting the AUC against the effect size at 10 and 25%DEG (Fig 2B). The hydra method performed better across all effect sizes, achieving near perfect performance above an effect size of 2.0 and 0.75 at 10 and 25%DEG, respectively. ssGSEA and GSVA performed similarly at low effect sizes, but ssGSEA performed better than GSVA as the effect size increased. Overall, the hydra framework performed significantly better than these standard gene set enrichment approaches, particularly at low effect sizes. Therefore, the hydra approach is better suited for subtyping within a disease cohort when the effect sizes are smaller and fewer genes are differentially expressed.

We performed a runtime analysis comparing hydra sweep, ssGSEA, and GSVA for identifying a single differentially expressed gene set, since these methods are directly comparable. Training the hydra model was the most computationally expensive step, but the classification of new samples was very fast. The average runtime for the hydra sweep algorithm was similar to ssGSEA, but the hydra runtimes were more variable across effect-sizes and number of differentially expressed genes. The GSVA approach was faster than hydra sweep and ssGSEA, but GSVA performed worse on the synthetic data analysis than ssGSEA and hydra. We repeated the above analysis with an effect size of 1.0, a %DEG of 25%, and a range of sample sizes, including 50, 100, 200, 300, 400, 500, 1000 samples. The hydra sweep and GSVA methods scaled well to large sample sizes, but the ssGSEA runtime increased exponentially as the sample size increased (Fig 2C & 2D).

Hydra analysis of high-risk neuroblastoma

High-risk neuroblastoma is an aggressive disease and is resistant to intensive therapy. Further subtyping of high-risk neuroblastoma may identify novel therapeutic targets and improve risk stratification. We hypothesized that unsupervised clustering of multimodally expressed genes associated with enriched Gene Ontology terms would identify expression subtypes of high-risk neuroblastoma tumors. TumorMap analysis [5] showed that the MYCN-non-amplified (MYCN-NA) neuroblastoma samples clustered separately from MYCN-amplified (MYCN-A) and stage 4S neuroblastoma samples (S2 Fig). We focused on the MYCN-NA neuroblastoma tumor samples because this is the largest set of samples (N = 70) and variation within MYCN-NA tumors is not well understood [28].

We applied the hydra filter analysis to the TARGET high-risk neuroblastoma cohort as described above. This analysis identified 931 genes within the MYCN-NA neuroblastoma cohort with a multimodal expression distribution. Of the 931 multimodally expressed genes, 358 genes were found to be potentially druggable by the Drug Gene Interaction Database (S1 File) and 60 genes were associated with an FDA-approved, anti-neoplastic drug [29].

We next examined whether unsupervised clustering of multimodally expressed genes revealed coordinated expression of annotated gene sets within the MSigDB database. Applying the hydra sweep command to the MYCN-NA neuroblastoma cohort discovered 105 gene sets with multimodal expression patterns. Each gene set sheds light on biological themes that are differentially expressed within the MYCN-NA neuroblastoma cohort. We clustered the differentially expressed gene sets to reveal these biological themes (S4 Fig). We found 6 major themes, including annotated cancer functions, cell cycle regulation, cell signaling pathways, immune functions, extracellular matrix reorganization, and metabolic pathway gene sets.

We applied the hydra enrich analysis to the MYCN-NA cohort to identify how the most highly enriched gene sets interact to form expression subtypes. This analysis found 428 genes with a minor component probability greater than 20% (S1 File). Gene Ontology analysis found enrichment for the following GO terms (FDR: q < 0.01): adaptive immune response (24 genes), mesenchyme development (12 genes), steroid hormone secretion (4 genes), and response to corticosterone (4 genes). DP-GMM analysis of the 44 enriched GO term genes identified three MYCN-NA neuroblastoma clusters (Fig 3A). The posterior probability for belonging to each cluster was 42%, 34%, and 17% for clusters 1, 2, and 3, respectively. The posterior probability for a sample belonging to a new cluster was about 6% in our analysis.

Fig 3

Hydra analysis identifies three distinct tumor microenvironment expression subtypes in MYCN non-amplified neuroblastoma samples.

A: Gene expression heatmap displaying expression profiles of hydra clusters. Heatmap columns (samples) are ordered by hydra cluster membership. Ward hierarchical clustering applied to rows (genes) identified coordinated expression of GO term genes. These GO term genes were originally identified by the hydra enrich command. B: GSEA performed on each cluster identified enrichment of tumor microenvironment and proliferative signaling gene sets. C: xCell enrichment score distributions for B-cells, CD8+ naive T-cells, and Fibroblasts, and the ESTIMATE TumorPurity score distributions for each cluster; enrichments for all cell types are available in S1 File. Abbreviations: Normalized Enrichment Score (NES), Epithelial to Mesenchymal Transition (EMT), Extracellular Matrix (ECM), Gene Ontology Biological Process (GOBP).

We next investigated cluster-specific expression signatures using GSEA (see Hydra Method section). Cluster 1 was enriched for adaptive immune response gene sets, cluster 2 was enriched for proliferative signaling gene sets, and cluster 3 was enriched for cancer-associated fibroblast gene sets (Fig 3B). Cluster 3 shares several features of a wound healing response, including fibroblast recruitment, extracellular matrix organization, and infiltration of immune cells [30].

Clusters 1 and 3 were enriched for tumor microenvironment-associated gene expression. To further validate this signal, we correlated the hydra clusters with enrichment scores from the tumor microenvironment profiling tools xCell [31] and ESTIMATE [32]. Cluster 1 had high average xCell enrichment scores associated with adaptive immune cell types including B-cells, CD4+ naive T-cells, and CD8+ naive T-cells (Kruskal-Wallis: p < 0.001). Cluster 2 was characterized by the absence of immune and stromal expression and higher tumor purity scores than clusters 1 and 3. The average ESTIMATE tumor purity was 88%, 96% and 82% for clusters 1, 2, and 3, respectively. Cluster 3 was enriched for fibroblast-associated expression by xCell analysis (Kruskal-Wallis: p < 0.001). Clusters 1 and 3 had higher ESTIMATE immune-associated expression levels than cluster 2 (average ImmuneScore per cluster: 58, -612, 56), but cluster 3 had the highest stromal expression signature score (average StromalScore per cluster: -1027, -1310, -135). Comparing ESTIMATE enrichment scores across clusters reveals clear trends in broad immune and stromal expression signatures. Lastly, we found a correlation between the hydra-identified tumor microenvironment subtype and CD274 and CTLA4 expression (S6 Fig).

We next correlated clusters with clinical features. We found no difference in patient survival outcomes across clusters (log-rank test, p > 0.05). Notably, cluster 1, which had the highest adaptive immune expression signal in MYCN-NA neuroblastoma, over-expresses cell-cycle regulation genes, which was not observed in other small blue cell tumors. We investigated associations with clinical covariates, including mutation burden, age, and tumor content as assessed by a clinical pathologist, but found no statistically significant differences (Kruskal-Wallis: p > 0.05). We then investigated associations between the hydra clusters and neuroblastoma-associated molecular aberrations and clinical features (S1 File). ATRX gene deletions were enriched in cluster 1 (Fisher’s Exact Test: p < 0.05). MKI low tumors were enriched in cluster 2 and 3 (Fisher’s Exact Test: p < 0.01). Chromosome 17 wild-type tumors were enriched in clusters 2 and 3 (Fisher’s Exact Test: p < 0.01). Analysis on a larger dataset may reveal additional clusters and correlations with clinical features.

Consensus clustering is a widely used approach for identifying tumor subtypes using gene expression data. We applied the M3C consensus clustering method, which is a more sophisticated version of consensus clustering that uses a null distribution to assess the statistical significance of the clustering [20, 21]. We used the top 5000 genes with the largest median absolute deviation (MAD) because this threshold is routinely used in unsupervised clustering of cancer gene expression data [33-35].

The M3C analysis resulted in the identification of two statistically significant clusters. One M3C cluster correlated with hydra clusters 1 and 3 and the other M3C cluster correlated with hydra cluster 2. Therefore, M3C clustering detected the tumor purity signal in the expression data, but was not able to separate the adaptive immune cell and fibroblast infiltrated clusters (hydra clusters 1 and 3). We also applied k-means clustering using the gap statistic approach [36, 37] for estimating the number of clusters, but this approach grouped all samples into a single cluster. We tested a range of MAD thresholds based on the median absolute deviation, but found similar results across thresholds (S3 Fig). Overall, the hydra approach was more sensitive at detecting distinct tumor microenvironment states than these other popular clustering methods.

To further investigate expression patterns within the hydra-identified tumor microenvironment subtypes, we performed GSEA by z-score normalizing each tumor’s gene expression data to its tumor microenvironment cluster. This is a novel GSEA approach that uses the tumor microenvironment state discovered by the hydra method to identify additional gene expression signals for individual samples. This approach revealed signals not present at the cohort level analysis (Fig 4). For example, enrichment of immune expression signatures within cluster 2 predicted differences in overall survival such that patients with higher immune expression had a better overall survival rate. Similarly, an elevated cell cycle signal within cluster 3 predicted worse survival compared to other cluster 3 samples with lower cell cycle expression. A metastatic expression signal was identified in the analysis of cluster 1 samples, but this signature did not correlate with a difference in survival. This approach may therefore provide appropriate background distributions for revealing and evaluating the significance of gene expression patterns and survival statistics within tumor subtypes.

Fig 4

Gene set enrichment analysis (GSEA) of MYCN-NA neuroblastoma identifies overall survival differences within hydra cluster 2 and cluster 3.

Cluster-level GSEA separated cluster 2 into high and low immune expression subtypes and cluster 3 into high and low cell cycle expression subtypes. A: Kaplan-Meier plot for immune expression subtypes within cluster 2. B: Kaplan-Meier plot comparing cell cycle expression subtypes within cluster 3.

N-of-1 tumor analysis for pediatric neuroblastoma

The command-line interface of the hydra toolkit includes a predict function for labeling samples using a pre-fit model. The MYCN-NA neuroblastoma model described above was used to predict expression subtypes on a new set of samples. We obtained tumor gene expression data from six stage 4, MYCN-NA neuroblastoma samples from the UCSC Treehouse gene expression compendium [5, 6]. The age at diagnosis ranged from 2 to 6 years. Four out of six samples had a deletion in the ATRX gene.

Application of the hydra N-of-1 analysis framework clustered 4 out of the 6 samples into cluster 1, which is characterized by adaptive immune cell expression. Three of the ATRX-deleted samples clustered with the high adaptive immune cell expression cluster (cluster 1) and one clustered in the low immune, high proliferative signaling cluster (cluster 2). We showed earlier that tumors with ATRX deletions tend to have higher adaptive immune expression, and we found a similar pattern in an independent set of MYCN-NA neuroblastoma samples.

Two of the samples with loss of ATRX came from the same patient but at different timepoints. The first sample (diagnostic sample) clustered with high adaptive immune cell expression (cluster 1), but the resection sample clustered with the low immune expression, high proliferative signaling cluster (cluster 2). We investigated possible explanations for the change in tumor microenvironment state. We performed GSEA comparing the samples from different timepoints to investigate potential mechanisms leading to immune evasion in these samples. GSEA found downregulation of the MHC Class I Antigen Processing & Presentation GO term in the resection sample (adjusted p-value < 0.002). Loss of antigen processing functions is a common mechanism of immune evasion across cancer types [38].

We obtained H&E stained sections for each of the hydra-identified clusters (S5 Fig). The cluster 1 sample had moderate levels of inflammation (30-50%) consisting of mature mononuclear cells, plasma cells, and eosinophils. The cluster 2 sample had minimal levels of inflammation (<10%) with some scattered mature mononuclear cells throughout the tumor. The cluster 3 sample looked similar to the cluster 1 slide with moderate levels of inflammation (30-50%), but also had regions of apparent necrosis. The inflammation and necrosis in the cluster 3 sample may correlate with the tissue remodeling/wound healing signature identified in the expression data.

Hydra analysis discovers complex tissue signatures

While the MYCN-NA neuroblastoma analysis above focused on immune and wound healing expression signatures, the hydra enrich method is unsupervised and can therefore detect any type of expression signature. To illustrate this, we applied the hydra filter/enrich analysis to the TARGET osteosarcoma cohort (N = 74) and discovered enrichment of the GO striated muscle contraction term (FDR < 0.01, Fig 5). Multivariate clustering for the GO striated muscle contraction gene set using the sweep routine identified two clusters. xCell analysis of the osteosarcoma cohort found significant enrichment of skeletal muscle expression in the second cluster (Mann-Whitney U test, p < 0.001). Surprisingly, the M3C clustering approach was not able to detect the strong muscle signature using the 5000 genes with the largest MAD (p > 0.05). We used the muscle expression signature to identify osteosarcoma tumors in the UCSC Treehouse Compendium which also contained a similar expression signature. We subsequently confirmed with a licensed pathologist that one of the muscle-expression positive tumor samples did contain significant muscle tissue infiltration. The hydra enrich analysis revealed expression signatures not routinely investigated when analyzing osteosarcoma data. Nevertheless, these signals contribute significantly to the tumor expression profile, so explaining these sources of variation is necessary to derive clinically relevant conclusions from gene expression data.

Fig 5

Hydra analysis of TARGET osteosarcoma cohort reveals skeletal muscle signature.

Hydra enrichment analysis on the TARGET osteosarcoma cohort revealed a subset of patients with high skeletal muscle expression. A: Clustered heatmap shows the muscle signature genes identified by hydra unsupervised enrichment analysis (purple: enriched for muscle signature; yellow: not enriched for muscle signature). B: xCell tumor microenvironment profiling identified significant differences in skeletal muscle expression compared to background (p < 0.001). C: H&E stained tumor slide confirms presence of striated muscle tissue within the tumor sample.

We applied the filter method to Ewing sarcoma and discovered multimodal expression of an important druggable gene, JAK1. Applying the multimodal expression model allowed us to deconstruct the Ewing sarcoma distribution into three components (S7 Fig). We found that the expression component with the highest JAK1 expression was also enriched for mast cell expression. Therefore, overexpression of JAK1 may not correspond to activation of the JAK/STAT signaling pathway in cancer cells but rather to the presence of mast cells within the tumor microenvironment. Furthermore, targeted inhibition of JAK1 using ruxolitinib was shown to inhibit essential mast cell functions, including degranulation [39]. Therefore, therapeutic intervention intending to inhibit JAK1 expression in cancer cells may inadvertently inhibit the patient’s mast cell functions. Overexpression analysis using the Ewing sarcoma JAK1 expression distribution may identify JAK1 as an actionable lead, but further investigation into the effect of inhibiting off-target JAK1 expression in mast cells is needed. The hydra framework facilitates the identification of important expression signatures which can be used to deconstruct complex tumor expression subtypes and identify potentially confounding expression signals.

We next quantified the number of multimodal druggable genes from the MYCN-NA neuroblastoma dataset that correlated with at least one xCell cell type signature. Out of the 358 druggable genes, we found that 77 correlated with a non-cancer cell type (Kruskal-Wallis test: Holm-Sidak adjusted p-value < 0.05, S1 File). Some of the druggable genes were expected to correlate with non-cancer cells, including the cytokines IL6 and TGFB2, which correlated with epithelial cells and fibroblasts, respectively. Other druggable genes were surprising, like AURKA and AURKB, which correlated with higher Th2 cell expression. Aurora kinases play essential roles in spindle formation during mitosis and the overexpression of these genes is associated with evading spindle formation checkpoints in cancer [40], but little is known in how these genes correlate with infiltrating immune cells. Aurora kinase inhibitors show limited clinical activity in solid tumors, but have been shown to have a greater effect in leukemias [40, 41].

Hydra analysis reveals recurrent expression subtypes across small blue round cell tumors

We next investigated whether similar hydra clusters could be identified across other small blue round cell tumors. We chose to focus on extracranial solid tumors because they are among the most common pediatric cancers, making up 20% of all pediatric cancer diagnoses [42], and while survival rates have improved, there are few effective treatment options for the subset of patients with relapse or refractory disease [43]. Identifying expression subtypes for these diseases may improve risk stratification and discover opportunities for new therapies. These tumors also share similar histopathological features, so we hypothesized that these tumors may share similar gene expression subtypes, despite significant differences in the raw expression profiles (Fig 6A).

Fig 6

Hydra enrich analysis of small blue round cell tumors reveals similar expression subtypes across cancer types.

A: TumorMap visualization of 6 small blue round cell tumor types. B: Hierarchically clustered heatmap for the top 10 enriched gene sets across the 21 small blue round cell tumor expression subtypes. Each column corresponds to a cancer type and an expression subtype (x-axis). Each row corresponds to a gene set. The expression subtype was manually assigned after reviewing the most highly enriched gene sets for each cancer expression subtype.

We first performed TumorMap analysis, which is a dimensionality reduction approach for visualizing genomic data on a 2D surface [5]. We found that small blue round cell tumor types—MYCN-NA neuroblastoma, osteosarcoma, Ewing sarcoma, synovial sarcoma, alveolar rhabdomyosarcoma, and embryonal rhabdomyosarcoma—all form separate TumorMap clusters (Fig 6A). This suggests there is a strong cell-of-origin signal driving the clustering of these cancer types, which is an observation that was recently made in the larger TCGA dataset of adult cancers [44]. While pan-cancer analysis emphasized the differences across small blue round cell tumors, we hypothesized that expression subtypes within cancer types would participate in shared biological themes.

We next performed hydra enrich analysis within each small round blue cell cancer type and found shared biological themes across all six small blue round tumor types. Hierarchical clustering of the top 10 statistically significant gene sets for each cancer type resulted in clustering by expression subtype and not the cancer type (Fig 6B). Common themes emerged across diseases including translational regulation, cell cycle regulation, immune effector cell signaling, inflammation, extracellular matrix organization, and tissue-of-origin signals. Furthermore, these signals predicted differences in patient outcomes in osteosarcoma and synovial sarcoma (Fig 7). In both cases, the presence of immune-associated expression correlated with better patient outcomes compared to tumors with proliferative signaling pathways associated with translation initiation and cell cycle regulation. Other osteosarcoma clusters were not included in the survival analysis due to insufficient number of samples with survival data (n < 5). Survival data were not available for the rhabdomyosarcoma and Ewing sarcoma expression datasets.

Fig 7

Hydra analysis identifies tumor microenvironment expression subtypes that correlate with patient outcomes in osteosarcoma and synovial sarcoma.

A: Kaplan-Meier plot showing overall survival curves for osteosarcoma wound healing and translation clusters. B: Kaplan-Meier plot showing metastasis survival curves for synovial sarcoma clusters.

Discussion

The hydra framework uses model-based clustering to facilitate the discovery of recurrent expression patterns within cancer gene expression cohorts. We leveraged recent improvements in model-based clustering algorithms to identify differentially expressed genes without a matched normal distribution. We modeled differential expression as a multimodal Gaussian distribution using nonparametric Bayesian statistics. We then enriched for biologically-annotated Gene Ontology terms and performed multivariate clustering to reveal expression subtypes. The hydra framework can be used for both identifying expression subtypes within large cohorts and classifying new tumor gene expression profiles using the trained models. The hydra framework outperformed standard gene set enrichment tools for identifying overexpression of the MSigDB Hallmark cancer gene sets in synthetic data. Application of this framework to small blue round cell tumors identified shared biological themes associated with the tumor microenvironment.

Multivariate gene expression analysis is typically underpowered because the number of genes greatly exceeds the number of samples. To address this limitation, we propose selecting for multimodally expressed genes before performing multivariate analysis. The hydra filter method reduces the number of genes and enriches for genes that participate in known biological processes, including those curated in the Gene Ontology and MSigDB databases. Selecting for multimodally expressed genes improves separation of known clinical subtypes better than the standard approach of using all expressed genes according to TumorMap analysis (S2 Fig). We also showed that the hydra approach of subsetting to multimodal genes improves detection of differential pathway expression, including the identification of expression subtypes associated with the TME.

Significant progress has been made in subtyping neuroblastomas and adapting therapy for aggressive subtypes, but unexplained heterogeneity remains [28]. Failure to account for this heterogeneity decreases the power of standard methods to detect important expression patterns. Identifying biomarkers using genome-wide technology may lead to improved risk stratification and the discovery of novel drug targets. Hydra analysis of the TARGET MYCN-NA neuroblastoma cohort found differential expression of tumor microenvironment markers, including markers of the adaptive immune response. Pediatric cancers are generally thought to be less immunogenic because they have lower mutation burdens than adult cancers, but the immunogenicity of pediatric cancer has not been sufficiently investigated [11, 12].

Our analysis found significant variation in immune marker expression, including markers of response to checkpoint blockade therapy, and identified ATRX deletions as a potential biomarker of immune infiltrated tumors in MYCN-NA neuroblastoma. Analysis of other small blue round cell tumors revealed similar expression signatures across tumor types, despite samples clustering by their histology in a pan-cancer TumorMap analysis. Identification of shared expression signatures across cancer types may suggest that these patients would respond similarly to therapies that target these pathways. In particular, the identification of a cross-disease subtype associated with high expression of immune markers may warrant further investigation of immunotherapies in small blue round cell tumors using a basket clinical trial design [45].

Hydra analysis found significant differences in tumor immune and stromal expression that may inform precision medicine applications. The tumor microenvironment has become an important therapeutic consideration, but few methods account for the tumor microenvironment directly. Tumor purity has been identified as a confounding factor in cancer gene expression subtyping efforts [46]. For example, tumor purity and tumor microenvironment expression have been shown to correlate with pancreatic cancer subtypes [47]. Furthermore, Aran et al. (2018) found that tumor purity was correlated with the mesenchymal glioblastoma subtype and recommended a differential expression approach to computationally remove the tumor purity signal. However, standard approaches for subtracting the tumor purity effect may not be ideal because several mechanisms may influence tumor purity, and each mechanism may result in a different expression pattern. For instance, our analysis of MYCN-NA neuroblastoma identified two gene expression signatures that correlated with lower predicted tumor purity. Cluster 1 had an adaptive immune expression signature and cluster 3 had a cancer-associated fibroblast signature. Therefore, the estimated tumor purity signal should not be subtracted without first accounting for the different mechanisms influencing tumor purity.

We also found shared biological pathway enrichment across small blue round cell tumors. While these diseases are related and may derive from similar cell lineages, current expression methods often emphasize difference across these diseases (Fig 6A). Unsupervised clustering of adult cancer types found that cell-of-origin signals strongly influence clustering of cancer gene expression data [44]. Although these diseases have distinct expression patterns on the surface, we discovered common themes once we subset the data to the cell-of-origin signal and applied the hydra analysis tools.

We found at least three shared TME states: immune silent, immune infiltrated, and wound healing subtypes. The wound healing subtypes predicted better overall survival in osteosarcoma and delayed metastases in synovial sarcoma tumors, which suggests the involvement of the host immune response limits the progression of these tumors. Amplification of the host immune response may further limit tumor growth and lead to immune-mediated tumor cell death. Additional research into immune modulating therapies is warranted in small blue round cell tumors and may lead to improved outcomes for some patients.

Conclusion

Precision oncology aims to differentiate tumors of the same diagnosis in order to match patients with the best treatment. We have developed the hydra framework to discover subtle but recurrent expression patterns within a cohort of samples with the same diagnosis, which is a novel strategy for pediatric precision oncology research. Our approach may help to uncover the biology underlying tumor progression and response to therapy. We have shown that hydra is more sensitive than standard gene set enrichment approaches for detecting differential pathway expression. Additionally, our framework provides tools to conduct unsupervised clustering analysis to discover expression subtypes. We applied the unsupervised hydra analysis to small blue round cell tumors and discovered distinct tumor microenvironment (TME) states. This shows that one of the strongest signals in clinical gene expression data comes from the TME, so careful modeling of the TME is required to maximize the impact of clinical gene expression analysis. The hydra framework provides unbiased clustering tools to characterize these sources of variation in specific disease populations and identify shared biological themes that can potentially be targeted therapeutically.

Example of bnpy memoized online variational inference clustering on toy data.

We used the bnpy moVB algorithm to infer the number of clusters from synthetic data. The model first randomly assigns clusters. Then, the model iteratively improves the model fit, creating and destroying clusters until the model converges on the correct number of clusters at lap 16 [56].

(TIF)

Click here for additional data file.

Enriching for multimodally expressed genes improves clustering of established neuroblastoma subtypes.

Standard TumorMap analysis of the TARGET neuroblastoma dataset resulted in stage 4S samples clustering with stage 4 neuroblastoma samples (left). An alternative TumorMap based solely on 1,498 multimodally expressed genes separated the stage 4S samples into a distinct cluster (right).

(TIF)

Click here for additional data file.

Consensus and k-means clustering applied to TARGET MYCN-NA dataset.

We tested a range of gene expression variation thresholds based on the median absolute deviation, but found that the clusters identified by this approach could not resolve the same clusters as the hydra approach. The barplot shows the number of clusters and the lineplot tracks the Rand index comparing the M3C and k-means clusters and the hydra clusters.

(TIF)

Click here for additional data file.

Hydra sweep analysis reveals differential pathway expression within MYCN-NA neuroblastoma without a matched cohort of normal tissue.

Unsupervised clustering of multimodal gene sets revealed biological themes associated with hallmark cancer functions, including cell cycle, immune cell signaling, extracellular matrix organization, and metabolism.

(TIF)

Click here for additional data file.

Hydra method correlates with distinct tumor features as assessed by licensed pathologist review of tumor H&E slides.

A-B: H&E sections from fresh frozen tumor tissue from MYCN-NA neuroblastoma sample at A: 2X magnification and B: 20X magnification. Tumor cells are medium to large with moderate amounts of cytoplasm and areas of rhabdoid appearing undifferentiated cells. There is a moderate amount of mixed inflammation present (30-50%) consisting mostly of mature mononuclear cells with some plasma cells and scattered eosinophils. C-D: H&E sections from fresh frozen tumor tissue from MYCN-NA neuroblastoma at C: 2X magnification and D: 20X magnification. Tumor cells are moderate to large in size with moderate amounts of cytoplasm. There is a minimal amount (<10%) of apparent mononuclear inflammation scattered throughout the tumor. E-F: H&E sections from fresh frozen tumor tissue from MYCN-NA neuroblastoma sample at (E) 2X magnification and (F) 20X magnification. Tumor cells are medium to large with moderate amounts of cytoplasm and areas of rhabdoid appearing undifferentiated cells. There are also areas of apparent necrosis. There is a moderate amount of inflammation present (30-50%) consisting mostly of mature mononuclear cells with some plasma cells and scattered eosinophils.

(TIF)

Click here for additional data file.

Hydra enrich analysis identifies correlation between expression subtypes and checkpoint blockade markers in MYCN-NA neuroblastoma.

(TIF)

Click here for additional data file.

Hydra analysis identified JAK1 expression clusters that correlate with mast cell expression signature in Ewing sarcoma.

A: JAK1 expression distribution for Ewing sarcoma cohort (top) and the JAK1 expression distributions for cluster 1 (green), 2 (orange), and 3 (blue). B: Boxplot showing the xCell mast cell enrichment score for the three clusters associated with JAK1 expression.

(TIF)

Click here for additional data file.

TARGET MYCN-NA neuroblastoma supplementary data.

(XLSX)

Click here for additional data file.

Hydra method documentation.

(PDF)

Click here for additional data file.

11 Dec 2019

Dear Dr Pfeil,

Thank you very much for submitting your manuscript, 'Hydra: A mixture modeling framework for subtyping pediatric cancer cohorts using multimodal gene expression signatures', to PLOS Computational Biology. As with all papers submitted to the journal, yours was fully evaluated by the PLOS Computational Biology editorial team, and in this case, by independent peer reviewers. The reviewers appreciated the attention to an important topic but identified some aspects of the manuscript that should be improved.

We would therefore like to ask you to modify the manuscript according to the review recommendations before we can consider your manuscript for acceptance. Your revisions should address the specific points made by each reviewer and we encourage you to respond to particular issues Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.raised.

In addition, when you are ready to resubmit, please be prepared to provide the following:

(1) A detailed list of your responses to the review comments and the changes you have made in the manuscript. We require a file of this nature before your manuscript is passed back to the editors.

(2) A copy of your manuscript with the changes highlighted (encouraged). We encourage authors, if possible to show clearly where changes have been made to their manuscript e.g. by highlighting text.

(3) A striking still image to accompany your article (optional). If the image is judged to be suitable by the editors, it may be featured on our website and might be chosen as the issue image for that month. These square, high-quality images should be accompanied by a short caption. Please note as well that there should be no copyright restrictions on the use of the image, so that it can be published under the Open-Access license and be subject only to appropriate attribution.

Before you resubmit your manuscript, please consult our Submission Checklist to ensure your manuscript is formatted correctly for PLOS Computational Biology: http://www.ploscompbiol.org/static/checklist.action. Some key points to remember are:

- Figures uploaded separately as TIFF or EPS files (if you wish, your figures may remain in your main manuscript file in addition).

- Supporting Information uploaded as separate files, titled 'Dataset', 'Figure', 'Table', 'Text', 'Protocol', 'Audio', or 'Video'.

- Funding information in the 'Financial Disclosure' box in the online system.

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com  PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

We hope to receive your revised manuscript within the next 30 days. If you anticipate any delay in its return, we ask that you let us know the expected resubmission date by email at ploscompbiol@plos.org.

If you have any questions or concerns while you make these revisions, please let us know.

Sincerely,

Florian Markowetz

Associate Editor

PLOS Computational Biology

Feilim Mac Gabhann

Editor-in-Chief

PLOS Computational Biology

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact ploscompbiol@plos.org immediately:

[LINK]

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: Review: Hydra: A mixture modelling framework for subtyping pediatric cancer cohorts using multimodal gene expression by Pfeil J et al.

Purpose of the research: Develop an improved computational platform to molecularly classify rare pediatric cancers with high accuracy by transcriptomic data that can be implemented in clinical settings.

Problem: Current modalities are insufficient to accurately identify novel, subtle cancer specific precision oncogenic targets and networks above thresholds relating to tumour heterogeneity from the microenvironment and the cell-of-origin.

What are Pfeil J et al trying to accomplish: Create a highly accurate computational analysis package that can be applied to single samples and easily utilized in clinical settings to facilitate identifying better treatment options for rare pediatric cancers.

How is it done today: The authors highlight several of the numerous current approaches for analysing bulk RNAseq data from patient samples and their current limitations. Importantly, they compared their new package to two of the current standards: single-sample gene set enrichment anlaysis (ssGSEA) and gene set variation analysis (GSVA). Based on the defined metrics, their Hydra approach outperformed the other packages overall.

What’s new in their approach: The authors elicit a Dirichlet process mixture model to detect multimodally distributed genes and signatures. By subsetting the data via this method, much of the background noise from the tumour microenvironment heterogeneity is reduced. The authors highlight that their approach was previously not feasible computationally for clinical setting but with new developments in Bayesian variational inference it is now possible. Their approach specifically performs better when effect sizes are smaller and fewer genes are differentially expressed which is an ideal platform for rare tumour types such as pediatric cancers.

Overall impression/suggestions:

- The authors were very thorough in the explanations and rationale for the methodology which is extremely helpful for individuals with limited knowledge in computational biology. The rationale for the study is sound and the approach is very reasonable.

- The abstract and author summary are concise, clear and accurately reflect the data presented in the manuscript.

- Would the authors elaborate as to why they chose TCGA GBM data to test hydra on initially? This a tumour predominantly found in adults and that has higher mutational frequency than pediatric brain tumours.

- The authors analysed 5 pediatric tumour datasets from UCSC treehouse compendium. It is very good to see that hydra and its different components can be applied to all the tumour types. However, it was confusing when performing some analyses on only 1 type of cancer or a few and not others. For continuity, would the authors consider focusing on neuroblastoma for each of the methods presented and use the other cancer types as a more supporting role? E.g. hydra filter/enrich also on neuroblastoma data.

- Have the authors considered applying hydra on the more common pediatric tumours such as lymphoma, leukemia or brain tumours? For example medulloblastoma has very well defined molecular subtypes based on transcriptomic data (methylomic, proteomic, et..)…would hydra better refine these currently well-defined subgroups by Paul Northcott, Michael Taylor, Stefan Pfister, Richard Gilbertson, etc…

- The histological data (images and analysis) are underwhelming to support their models. For example, figure 3C would be much better supported if the authors had access to human patient samples to compare the findings of the xCell score based on hydra clusters. Immunohistochemistry for B-cells, CD8+ T-cells and fibroblasts are easy to do if samples are available. For the immune system, the spatial context is especially important when determining the ‘hot’ and ‘cold’ nature of the tumours.

- The N-of-1 tumour analysis is very attractive. Rates of pediatric cancers are low which is a hurdle in facilitating better treatments for those affected. Having a more rigid histological analysis of the samples used (figure 4) is needed. General lymphocyte content is in not sufficient for inflammation. The pathologist could supply additional information on neutrophil, macrophages, etc… to better describe the type of inflammation occurring. This level of pathological resolution is needed to better support the findings from Hydra to be utilized on n-of-1 findings.

- The hydra filter/enrich analysis of the osteosarcoma data is very interesting. I think these data would be better suited in the main figures not supplemental.

- For actionable drug targets, the authors highlight JAK1 signaling in Ewing Sarcoma. The authors highlight how important interpretation of the data are (Mast cell infiltration vs. the tumour) for selecting drug targets regarding the complexities of the tumour microenvironment. Can the authors supply a more concrete cancer specific actionable target to support the filter method of Hydra…something less nuanced?

- For S4B and C, these data are impressive and should be moved to main figures.

- Why is S7 in the supplemental data? Did the authors assess all cluster combinations and identify only wound healing and translational regulation as important in osteosarcoma…what about cell cycle,etc..

General comments:

- Line 22: Though many accept that pediatric tumours have fewer mutations relative to adult cancers, this is not true for all pediatric cancers. Additionally, the spectrum of tumours in paediatrics and adults is remarkably different. Would the authors supply some citation for their statement or provide a comparison on mutational frequencies…even focusing on brain tumours comparing peidatrics to adults. Stefan Pfister’s group in Heidelberg have published some data related to this.

- It may be my version, but the figures were very blurry and hard to make out details.

- Line 473: these findings are not too surprisingly based since the 5 cancers assessed have links back to neural crest cells.

- Line 482: should reference S7 not S6.

- The method section describing the histological analysis is not sufficient. Please see above for suggestions on what to highlight.

- Several graphs in the figures have no labelled axes.

Reviewer #2: In the manuscript “Hydra: A mixture modeling framework for subtyping pediatric cancer cohorts using multimodal gene expression signatures“, the authors proposed a novel method for classifying tumor samples based on their expression profiles. The method focuses on multimodally expressed genes (identified through a Dirichlet process mixture model) and its effects are demonstrated in cohorts of pediatric cancers. The method is compared to a standard method in the field, gene set enrichment analysis (GSEA) and its equivalent for N-of-One type of analysis, single-sample GSEA. The method is a non-parametric mixture model and contain three tools named filter, enrich and sweep. The first result is obtained from synthetic gene expression data and demonstrate that hydra outperformed ssGSEA and GSVA (Gene set variation analysis). The second one identified 3 subclusters of the MYCN non-amplified neuroblastoma samples and relates these 3 subclusters to features of the tumor microenvironment. The third result illustrate the use of hydra to classify individual samples based on an annotated compendium of expression profiles, with five MYCN non-amplified neuroblastoma samples and the classification provided in the previous section. The fourth result section demonstrate the use of hydra’s functions for unsupervised pathway enrichment on the TARGET osteosarcoma dataset and on a Ewing sarcoma dataset. The final result shows how the hydra analysis on small blue round cell tumors clustered samples based on shared expression profiles and regardless of histology, while TumorMap clustering separates these different tumor types.

The method is well described and the results nicely illustrate what it can achieve, on cohorts that are both accessible and of interest to the authors of the manuscript. I have a set of minor comments and recommendation to try to improve the manuscript.

Comments:

1) I think a discussion of runtime is warranted, especially in comparison with the other methods, and an estimation of how runtime scales with the number of samples in the analysis.

2) In the first result with synthetic data, the data is modelled “as a multivariate Gaussian distribution” (l.204), which is the running assumption in Hydra (l. 94). It would be more powerful to model the synthetic data with different assumptions and measure the performance of the method.

3) l.125: The URL www.dockerhub.com/jpfeil/hydra failed for me, but the correct URL seems to be https://hub.docker.com/r/jpfeil/hydra. I note that the title of the manuscript indicate that hydra is a method for “subtyping pediatric cancer cohorts”, but the short description of the Docker package describes it as a “clustering pipeline to identify differentially expressed pathways”. I get how those two things are related, but I believe it would be good to clarify in the text .

4) Some sections are very technical. For instance, the lines 111 to 116 are referring to concepts and methods that may not be common knowledge for part of the audience. I leave it to the authors to decide if they can explain it more or prefer to leave it there, assuming that the reader will consult the references cited.

5) I would appreciate a guide, on the Github page or even better as part of the manuscript, describing how to use the method on a completely novel cancer type and novel datasets. How to process the data, build the model and run the different tools of Hydra? Some documentation already exists, but could be expanded a bit more if the authors want more users to adopt their method.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

28 Feb 2020

Dear Mr. Pfeil,

We are pleased to inform you that your manuscript 'Hydra: A mixture modeling framework for subtyping pediatric cancer cohorts using multimodal gene expression signatures' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Florian Markowetz

Deputy Editor

PLOS Computational Biology

Feilim Mac Gabhann

Editor-in-Chief

PLOS Computational Biology

***********************************************************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: Thank you for addressing all the comments/concerns raised during the initial review.

The comments to the questions and alterations to the manuscript have satisfied this reviewer.

Reviewer #2: Thanks for taking into consideration my comments and for improving your manuscript. I recommend its publication and believe it will be well received by the scientific community

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the

PLOS Computational Biology

data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review?

For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: Yes: Eric P Rahrmann

Reviewer #2: No

1 Apr 2020

PCOMPBIOL-D-19-01924R1

Hydra: A mixture modeling framework for subtyping pediatric cancer cohorts using multimodal gene expression signatures

Dear Dr Pfeil,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Laura Mallard

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol