dataesr/scientific-paragraphs-categorization huggingface dataset:
55-0. 56). Likewise, the WGCNA approach identified 322 genes in informative modules for docetaxel. Here we also found that LOOCV random forest models built from these genes had low predictive power (AUC 0. 6-0. 79). As previously mentioned, the WGCNA approach may overlook modules with informative genes whose signal is overwhelmed by less-informative genes, because it utilizes module eigengenes to determine if a specific module is related to the phenotype of interest (in this case response to chemotherapy). Once module eigengenes are used to identify informative module, the standard WGCNA workflow then identifies hub genes that may be particularly informative. Following this approach (see STAR Methods), we identified 14 hub genes for carboplatin (AUC with a GLM LOOCV 0. 62-0. 77) and 27 hub genes for docetaxel (AUC with a GLM LOOCV 0.
58-0. 68). Although these genes are representative of each module's eigengenes, they are not necessarily the genes that carry the predictive power. To address these issues, we applied a novel adaptation of CTD to determine which sets of genes in these large WGCNA modules were significantly connected in the pruned networks. CTD is an information-theoretic-based network method that identifies patterns of connectedness between analytes and assigns p values to highly connected subsets of analytes within large modules without computationally costly permutation testing (Figure  4A ). 65 With respect to gene expression data, CTD is used to ''connect the dots'' and identify these highly connected gene sets by finding subsets that are more connected in a weighted graph than by chance. Because response has not been added as a node in the expression graph, we hypothesize that significantly connected sets within the large modules may be connected due to their latent connection to response (Figure  4B ).
For each of the pruned networks, we identified small CTD submodules, which then allowed for the generation of generalized linear models (GLMs), and the identification of specific genes that are predictive for differential response to the two agent classes chosen. The p value calculated by CTD for these modules, as well as the significance of their linear regression with response, and whether the larger module was identified as informative by WGCNA can be found in Table  S1  related to Figure  4. In Table  S2  (related to Figure  4 ), the genes used for the multigene biomarkers and their modules are identified. An example of an informative submodule for the carboplatin epithelial network is shown in Figure  4C , whereas an example of an informative docetaxel epithelial submodule is shown in Figure  4D.
Comparison of the CTD/WGCNA approach with other feature selection methods We then investigated how the combined CTD/WGCNA approach was compared with six widely used feature selection methods (WGCNA hub genes,  52  correlation,  49  recursive feature elimination with a random forest model,  49  Boruta,  50  F-test,  51  and Mutual Information Regression). 51 To identify the WGCNA hub genes, we identified modules whose eigengenes were correlated with response and then determined which genes had the highest significance. For the correlation and recursive feature selection approaches, we identified the top 10 genes (as our CTD/WGCNA approach identified 9 genes for carboplatin and 6 genes for docetaxel) that were related to response in both the stromal and epithelial compartments. Boruta 50 is a random forest wrapper method that was applied using the Boruta function in R to  identify informative genes.
We also used the SelectKBest function from the python scikit-learn package to perform the F-test and mutual information regression analysis. The UpSetR package was then used to visualize the number of unique genes identified by each method as well as the number of genes multiple methods identified (Figures  5C  and 5D ). Although a few genes were identified with multiple methods, most potentially informative genes were unique to the feature selection method used. Genes that were identified as potentially informative by at least three feature selection methods for carboplatin included SV2A (5 methods), MAGED4 (4 methods), MAGED4B (4 methods), TRHDE (4 methods), and MED12L (3 methods). Genes that were identified by at least three feature selection methods for docetaxel included only ASGR1 (4 methods) and CRMP1 (3 methods).
CTD/WGCNA outperforms other methods of feature selection To evaluate CTD as an informative feature selection method, we tested each of these feature selection methods for their performance predicting treatment response using ''pseudo-bulk'' RNAseq data reconstructed from the PDX epithelial and stromal RNAseq gene expression datasets. We built these pseudobulk mRNA expression profiles by summing the epithelial and stromal mRNA expression values for common genes to approximate human clinical bulk RNAseq. To evaluate the performance of each feature selection method with a generalized linear model within our dataset, we compared the Mean Absolute Error (MAE) (Figure  5E ) and root-mean-square error (RMSE) (Figure  5F ) from the quantitative response predictions using a 70/30 training/test set splitting approach of the source data, which were reshuffled 10 times.
Preclinical validation of differential chemotherapy response predictors For the purposes of validation in other PDX datasets, and in clinical trial response assessments, we converted our quantitative response assessments to a qualitative modified RECIST 1. 1 classification for each PDX. We then determined how predictive our models were for both complete response (CR) versus all else (progressive disease [PD], stable disease (SD), a partial response [PR]) and CR/PR versus SD/PD using the full 45 PDX dataset, and how they were compared with models built using the informative module and hub genes identified with the standard WGCNA approach. ROC curves were constructed using LOOCV models for each of the treatments for two response cutoffs: complete response (CR) versus all else and complete response plus partial response (CR/PR) versus all else.
For both the CTD biomarkers and the hub genes identified by WGCNA, GLM LOOCV models were used to predict the qualitative response. However, due to the large number of genes identified as informative in the WGCA models, Random Forest (RF) models with default caret  49  parameters were used instead of GLMS to predict qualitative response to therapy. For carboplatin CR versus all else, the AUCs of these models are as follows: CTD = 0. 85, hub genes = 0. 77, and WGCNA modules = 0. 56 (Figure  6A ). Likewise, the AUCs for carboplatin CR/PR versus all else are CTD = 0. 80, hub genes = 0. 62, and WGCNA modules = 0. 55 Figure  6B ). All these methods had low predictive power for docetaxel CR (Figure  6C ), which we hypothesize may be associated with the unbalanced nature of responses in this collection (9 CR, 36 nCR).
For docetaxel CR/PR the AUCS were CTD = 0. 86, hub genes = 0. 69, and WGCNA modules = 0. 79 (Figure  6D ). Overall, we found that the CTD/WGCNA method was more predictive than the standard WGCNA approach and that both the carboplatin CTD/WGCNA multigene biomarker panel and the docetaxel biomarker panel were informative for both quantitative and qualitative response. To be maximally useful, biomarker panels predictive of differential drug response should be informative for choosing which drug might be more effective in which tumors. To evaluate our ability to choose which chemotherapy agent, docetaxel, or carboplatin, would have the best response in each PDX, we used a leave one out cross-validation (LOOCV) approach with the CTD multigene biomarkers to predict the quantitative response of each PDX to both carboplatin and docetaxel.
The ROC for the CTD genes is in red, the ROC for WGCNA hub genes is in green, and the ROC for the WGCNA modules is in blue. (B) ROC for carboplatin complete and partial response (CR/PR) versus all else. The ROC for the CTD genes is in red, the ROC for WGCNA hub genes is in green, and the ROC for the WGCNA modules is in blue. (C) ROC for docetaxel complete response (CR) versus all else. The ROC for the CTD genes is in red, the ROC for WGCNA hub genes is in green, and the ROC for the WGCNA modules is in blue. (D) ROC for docetaxel complete and partial response (CR/PR) versus all else. The ROC for the CTD genes is in red, the ROC for WGCNA hub genes is in green, and the ROC for the WGCNA modules is in blue. (E) The best response predictions for all PDXs. If both the actual best therapy and the predicted best therapy is carboplatin, the dot is blue.
Preclinical validation: Multigene biomarker panels are predictive for chemotherapy response in other PDX datasets To test the predictive power of our gene panels, we tested them using on other PDX datasets (Table  1 ) with either associated RNAseq (within platform) or Affymetrix array (across platform) gene expression data. With respect to within platform validation, the carboplatin and docetaxel were tested in a PDX dataset with RNAseq and qualitative drug response from the Rosalind & Morris Goodman Cancer Research Center (GSE142767). 42 Of the PDXs in this dataset (n = 37), 30 had mRECIST-like single-agent response data for cisplatin and paclitaxel. An LOOCV approach was used with the XGBoost tree model over the nine informative genes for carboplatin to predict cisplatin response, and a similar LOOCV approach was used over the six docetaxel genes to predict response to paclitaxel.
Despite the differences in the PDX composition and the specific therapies used (cisplatin versus carboplatin and paclitaxel versus docetaxel), we found that our multigene biomarker panels were predictive in this dataset. The AUCs for cisplatin response in this dataset were 0. 72 for CR versus all else and 0. 54 for CR/PR versus all else (Figures  7A  and 7B ). The AUCs for paclitaxel response were 0. 72 for CR versus all else and 0. 73 for CR/PR versus all else (Figures  7A  and 7B ). Thus, our biomarker panels had predictive power for therapies that were in the same class (taxane versus platinum-based). To test the cross-platform performance of our biomarker panel, the predictive power of our carboplatin biomarker was tested in a gene expression Affymetrix array dataset from the Curie Institute (n = 37) evaluating qualitative responses following treatment with cisplatin.
Clinical validation: Predictive power of the platinum multigene biomarker panel To investigate the predictive power of our PDX-based platinum response signature in a clinical dataset, we tested our carboplatin gene biomarker on an array dataset with 24 TNBC patient samples (GSE18864). 97 hese patients were treated with a cisplatin monotherapy and had response to therapy measured qualitatively with the RECIST 1. 1 classification system. 69 For this analysis, we also used an XGBoost tree model (LOOCV), this time using the nine informative genes for carboplatin response. The AUC for this model was 0. 73 (Figure  7D ), demonstrating that we have predictive power in a clinical dataset for platinum response, again across expression platforms. To determine if the expression of stromal genes indeed contributed to the predictive power of the clinical model, we determined the variable importance rankings for this model using the varImp function from caret.
Clinical validation: Predictive power of the taxane multigene biomarker panel Unfortunately, a single-agent taxane dataset could not be identified. However, we also tested the predictive power of our taxane biomarker panel for ''benefit from taxane'' by evaluating RNAseq data from a subset of TNBC patient samples from the BrighTNess clinical trial who were treated with paclitaxel followed by doxorubicin and cyclophosphamide using a LOOCV XGBoost approach. In this trial, response was evaluated only after the AC treatment was completed. Although full RECIST 1. 1 69 calls were not available, patients were evaluated for pCR versus not (GSE164458) (Loibl et al. , 2018). Although a different taxane agent was used (paclitaxel) and response was not evaluated immediately after the paclitaxel treatment, we still had some predictive power in this clinical dataset (AUC = 0.
CTD allows for the identification of biologically tractable multigene biomarker To identify our gene biomarker panels, we first had to address a methodological issue. The initial gene modules identified by WGCNA were large, and it was difficult to determine which specific genes were contributing to the modules being associated with response. Through the application of CTD, we identified highly connected sets of genes and a specific set of genes that are predictive of response to therapy. The biomarker panels identified by this approach outperformed several methods that are commonly used for feature selection including correlation analysis, a recursive feature extraction approach, and the hub genes identified with the standard WGCNA approach at predicting quantitative response. We also compared the performance of CTD and the standard WGCNA approach of identifying informative modules and hubs genes for predicting qualitative response to therapy.
Multigene biomarkers are predictive across expression platforms and within chemotherapy classes We then used our multigene biomarkers to predict response to therapy in other data collections. Because of biological differences between PDX models and patients, differences in data type (array vs RNAseq), and differences in chemotherapy agents (e. g. , carboplatin vs cisplatin), we built new models over our multigene biomarkers for each dataset. Several publications have demonstrated phenotypic and gene expression stability over multiple transplant generations, and in some cases, our response data were derived from PDX models at multiple passages (to be described in full elsewhere). Thus, a given PDX generally represents the tumor of origin with high fidelity. 40, 46, 47, 103, 104 However, in patients, several biological factors can lead to heterogeneity in breast cancer tumors and affect a tumor's response to therapy including tumor stage, location, grade/stage, and race/ethnicity (Mermut et al.
, 2021). Treatment of tumors also affects a tumor's expression pattern  105, 106  and its subsequent response to therapy. Therefore, the patient's prior treatment regimen may impact a PDX's response to therapy. Because these biological covariates may affect each of our datasets differently, we fix our multigene biomarkers, thereby allowing the models to be flexible, which, in turn, allows us to predict chemotherapy response in diverse datasets. By testing our multigene biomarkers across datasets from different cohorts, analyzed with different gene expression platforms, treated with different agents, we determined that (1) our multigene panels are predictive in both RNAseq and microarray datasets and (2) the multigene panels developed for carboplatin and docetaxel are predictive for other chemotherapy agents of the same class (platinums and taxanes, respectively).
Limitations of the study Although this study has shown that PDXs can serve as preclinical cohorts, there are elements of the human immune system that are not replicated in immunodeficient PDX models. Creating humanized PDX mice via immune system reconstitution may allow for further insight into how the immune system contributes to resistance to chemotherapy in TNBC. The use of humanized PDX models may also improve the concordance between patient and PDXs response to therapy. In addition, in the clinical paclitaxel dataset (Loibl  et al. , 2018), clinical response was not evaluated immediately after treatment with paclitaxel. Testing our taxane biomarker on a clinical dataset where response was measured after immediately following taxane treatment would allow us to determine better the predictive power of this biomarker for patient response.
Gene expression quantification RNAseq FASTQ files were processed using Xenome to separate murine stromal reads and human epithelial reads, which were both then quantified with rsem-calculate-expression. A reference index was created for RSEM (v1. 3. 0) 120 using rsem-prepare-reference on the hg38 and mm10 genome assembly (FASTA) and transcript feature (GTF) files. Expected counts from RSEM were collated to create a genes-by-samples matrix. This matrix was Upper-Quantile normalized per-sample to account for differences in sequencing depth.
Clustering of PDXs with primary tumors For the RNA-seq clustering of PDX with TCGA primary tumors, the log2(upper quantile normalized) RSEM counts that were classified as cancer reads in the PDXs by Xenome were clustered with log2(upper quantile normalized) RSEM counts from TCGA BRCA tumors. Clustering was done over the PAM50 genes with the hclust function from the stats v 4. 2. 1 package with Pearson's correlation coefficient. For the protein clustering, the z-scores of proteins were computed (TNBC vs non TNBC) for both the primary CPTAC tumors and the PDX tumors. These z-scores were then used to cluster the PDXs with the primary tumors with hclust over the 1,000 most variable proteins with Pearson's correlation coefficient using the stats v 4. 2. 1 package.
Network building To identify genes that could be predictive of response to chemotherapy, we used the network approach described below. For both the human and mouse reads, a network was built over all the PDXs and the responsive PDXs for both carboplatin and docetaxel independently. The gene expression data was obtained from deep RNAseq ($200M reads/sample) of the PDXs. Preprocessing included removing genes that were expressed in less than 20% of the samples or had a low expression across all the samples. The five thousand most variable genes were then selected for the human and mouse reads using the mad function from the stats v 4. 2. 1 package. These genes were used to build the species-specific networks using the WGNCA package v 1. 71. The following graphs were built: an epithelial graph with all the PDXs, an epithelial graph with the samples that are responsive to carboplatin, an epithelial graph with the samples that are responsive to docetaxel, a stromal graph from all the samples, a stromal graph with the samples that are responsive to carboplatin, and a stromal graph with the samples that are responsive to docetaxel.
WGCNA module discovery and the application of CTD Coarse modules were defined in the pruned graphs by applying WGCNA v 1. 71. 's cutreeDynamic function. After the initial module detection with the cutreeDynamic function, related modules were merged with WGCNA v 1. 71's mergeCloseModules function. The module WGCNA v 1. 71. 's moduleEigengenes function was then used to calculate the module eigengenes for each of these modules. To identify significantly connected submodules within these modules with CTD v 1. 1. 0 (https://github. com/ BRL-BCM/CTD),  65  we used the single-node diffusion walker (singleNode. getNodeRanksN function) CTD v 1. 1. 0 to get the ranks of each node. The mle. getPtBSbyK CTD v 1. 10 function was then used to encode each set of connected nodes, and the mle. getEncodingLength CTD v 1.
, 2021). To empower other users to utilize CTD, we have deployed CTD as a tool on the Seven Bridge's Cancer Genomic Cloud as a public application: (https://cgc. sbgenomics. com/public/apps/pedjaoetf/ctd-commit-project/ctd) Users will only need to upload their graphs as an adjacency matrix, as well as a list of genes in each module. CTD's code has been deposited in public GitHub repository at https://github. com/BRL-BCM/CTD where users can either download and use the package locally in R or run it inside Docker Desktop.
Comparison of our CTD/WGCNA approach to other feature selection methods To compare the MAE and RMSE of the CTD/WGCNA approach to other commonly used feature selection methods, we also choose potentially informative genes with three other commonly used feature selection methods (correlation, recursive feature extraction, hub genes from informative WGCNA modules). Correlation -for each treatment the top 10 most correlated epithelial and stromal genes were determined (base R v 4. 2. 1). Recursive feature extraction -for each treatment the top 10 most informative epithelial and stromal genes were determined (caret v 6. 0-93). Hub genes -to identify the hub genes for each treatment, WGCNA modules whose eigengene was correlation with response >0. 35 were identified (WGCNA v 1. 71). Then genes with a >0.
75 correlation with the module eigengene and a >0. 2 correlation with response were identified as the hub genes. For carboplatin there were 14 epithelial WGCNA hub genes and no stromal WGCNA hub genes. For docetaxel there were 26 epithelial WGCNA genes and 1 stromal WGCNA gene. Boruta -the Boruta R function was used to identify genes that were tentative and confirmed for importance in the response models (Boruta v 7. 0. 0). For carboplatin Boruta identified 22 genes, and for docetaxel Boruta identified 28 genes. Reticulate v 1. 26  121  was used to implement the python scikit-learn package in R. The scikit-learn v 1. 1 SelectKBest function was used to perform the F-test and mutual information regression feature selection methods with k = 10 (stromal and epithelial). The F-test approach identified 18 genes that are informative for carboplatin and 15 genes that are informative for docetaxel.
Rosalind & Morris Goodman Cancer Research Center RNAseq PDX dataset prediction An RNAseq dataset and qualitative response were obtained from GEO (GSE142767) (n = 37). The majority of these PDX models were Basal-like and triple negative, however there were some PDX models that were classified as HER2 (7) or positive for ER (5) or HER2 enrichment (3). A LOOCV XGBoost v 1. 6. 0 tree model was built over both the carboplatin genes (to predict response to cisplatin) and the docetaxel genes (to predict response to paclitaxel). A ROC curve was then calculated to determine the predictive power of predicting if the PDX has CR vs. all else or CR and PR vs. all else.
Curie dataset prediction The Curie dataset was obtained from the Curie institute. This dataset consists of the Affymetrix array expression of 37 TNBC PDXs and the human equivalent response to cisplatin. Due to the inconsistencies between the array and RNAseq platforms, an XGBoost v 1. 6. 0 tree model (LOOCV) was built over the carboplatin multigene biomarker with the xgbTree option in caret for this dataset to determine the AUC and predictive nature of these genes.
Dana-Farber cancer institute TNBC patient array dataset prediction The Dana-Farber (GSE18864) dataset was obtained from GEO. 97  This dataset included array data from 24 patient biopsies as well as qualitative response of these patients to neoadjuvant cisplatin. These patients were part of a cohort of 28 women with TNBC that were treated with neoadjuvant cisplatin every 21 days. An XGBoost v 1. 6. 0 tree model (LOOCV) was built over the carboplatin multigene biomarker with the xgbTree option in caret for this dataset to determine the AUC and predictive nature of these genes.
BrighTNess TNBC patient RNAseq dataset prediction RNAseq and patient response to paclitaxel followed by AC was obtained for the BrighTNess (GSE164458) dataset from GEO  38  Data from arm C, single agent paclitaxel followed by AC in 123 patients were used for this study. A LOOCV XGBoost v 1. 2. 0. 1 approach was used over the informative genes for to predict the response (pCR or not) to paclitaxel.