dataesr/scientific-paragraphs-categorization huggingface dataset:
Simple Summary: Gene expression monotonicity is an important feature in the evolution and progression of cancer, yet has not been investigated in bladder cancer (BLCA). Most of the transcriptomic stage investigations of BLCA are limited either to a subset of stages, or to a small number of samples. Here, we leverage publicly available data to create a meta-dataset of 1135 primary BLCA transcriptomes, to identify genes and processes with a monotonal change related to higher clinical or pathologic stages. Our analysis aims to deepen the current understanding of the disease's molecular pathogenesis, as well as to propose a prognostic gene signature based on the trait of monotonicity. Results demonstrate tumor dependencies on specific cell-cycle and metabolic microprocesses, and highlight an eight-gene signature capable of prognosing 5-year outcomes in both the discovery and validation sets.
Materials and Methods A comprehensive data mining strategy was employed to retrieve studies applying -omics technologies in BLCA. The overall workflow is summarized in Figure  1 , and described in detail in Appendix A. This analysis focused on primary, treatment-naïve tumor transcriptomes. The selection of datasets and samples eligible for the analysis, as well as the processing and integration of raw microarray data, is detailed in Appendix A. Briefly, we compiled a microarray discovery cohort from 12 well-characterized GEO datasets (GSE121711, GSE93527, E-MTAB-1940, GSE31684, GSE104922, GSE128959, GSE83586, GSE48276, GSE52219, GSE69795, GSE13507, and GSE48075). These data (summarized in Figure  1 ), comprised 1054 primary bladder cancer tumor transcriptomes of treatmentnaïve patients without any prior cancer history, along with profiles from 81 non-malignant urothelium tissues adjacent to the tumor site (NAU), corresponding to a total of 1135 gene expression profiles.
Stage distribution among the utilized datasets is shown in Table  A1. Two RNA-Seq datasets were used for validation (TCGA2017 and IMvigor210), as per availability. Table  1  shows sample allocation to the clinical variables for both the discovery and validation sets. Figure  1  illustrates the overall study design and workflow. In the discovery set, the ratio of men:women was 3. 5:1, with equal distribution between NMI and MI disease (p = 0. 99), and similar mean age at baseline diagnosis (68 years, p = 0. 81). Percentages of NAU, NMI, and MI in the dataset were 7. 1%, 43. 5%, and 49. 4%, respectively, with the grade distribution being as follows: 16. 5% low-grade and 48. 8% high-grade disease, with the remaining samples lacking available grade information. Detailed histological records were missing for 71.
5% of the cohort, with the most frequently reported histology among the available records being urothelial/papillary (23. 3%), and squamous differentiation being the most frequent variant (1. 3%). Gene expression across the 12 GEO datasets was harmonized using the empirical Bayes approach of the ComBat algorithm. Batch effect removal was assessed with the algorithm BatchQC (Figures  S1  and S2 )  [21] , with relative log expression and principal component analysis plots (Figure  S3 ), with comparative analysis of expression levels between housekeeping and non-housekeeping genes (Figure  S4 ), or using a set of 11 positive BLCA markers with known regulation across non-malignant-NMI-MI or non-malignantlow-grade-high-grade disease (Figure  S5 ). ComBat successfully eliminated batch effects (Figures  S1-S3 ), while maintaining the biological variability (Figures  S3  and S4 ), BatchQC reports indicated that the variability in the corrected data was explained by stage rather than by the batch variable (Figure  S2 ), allowing for an in-depth downstream analysis.
The  To identify genes that form a continuum of changes across BLCA stages, each of the 5 disease stages was initially compared against non-malignant adjacent urothelium (NAU). A total of 3018 genes differed significantly (Mann-Whitney p < 0. 05), while having the same orientation of change in all comparisons. We refer to this set of genes (n = 3108) as concordantly differentially expressed genes (CDEGs). CDEGs were utilized to infer pathway activation scores, and to create stage co-expression networks. Pathway activation scores per sample were calculated with the ssGSEA-GSVA method [22], using the molecular signature database libraries of Hallmark, Canonical Pathways (Reactome subset), C3 (GO biological processes subset), and C5 (GTRD subset of transcription factor targets). Dorothea (https://github.
Out of the 3108 × 3108 = 9,659,664 gene-pair weights calculated individually per condition (i. e. , NAU and 5 BLCA stages), gene-pairs with the highest GENIE3 weights, which were also positively correlated based on the Spearman's coefficient, were used to construct networks. The cutoff for this selection was determined based on the gene size of the resulting networks; to avoid network saturation, we opted to keep their gene sizes close to half the number of CDEGS (= 1554 genes), which resulted in setting the cutoff to the top 5600 gene-pairs. Networks were constructed with igraph and were analyzed with Louvain clustering  [26]  to identify local modules of coexpression relationships (communities) per condition. The top five in size (= number of genes) of the communities of co-expressed genes per condition were analyzed for Gene Ontology Biological Processes with clusterProfiler  [27].
Potential drug targets in the coexpression networks were defined based on the betweenness centrality metric  [28] , using default cutoffs (computed with igraph). We utilized CDEGs to extract genes whose expression levels was monotonically increasing or decreasing with higher stage. Monotonicity for a gene was defined as being a CDEG and additionally having a continuously larger/smaller fold change with increasing stage (fold change, as defined in each disease stage versus NAU). Functional annotation and enrichment were performed using PubMed and the online tool GeneCards (https: //ga. genecards. org/, accessed on 22 March 2021), respectively. Out of the monotonal subset, 43 genes were found to be of prognostic value (Cox univariate association with 5-year outcome). Eight of these genes, validated in the TCGA-BLCA dataset, were further utilized to construct a sample-wise scoring system-the 8-gene prognostic signature-by summing the expression values of the upregulated genes (n = 4) while subtracting the downregulated genes (n = 4).
Before calculations, each of the gene expression values per sample was divided by the gene's variation across the dataset, in order to minimize the effect of individual gene variability on the final signature score: S i = A i /VarA + B i /VarB + C i /VarC + D i /VarD -E i /VarE -F i /VarF -G i /VarG -H i /VarH (1) where S i denotes the sample-wise derived signature score, A i B i , C i , D i , and E i , F i , G i , H i , denote the sample-wise gene expression of the 4 upregulated and 4 downregulated genes, respectively, and Var denotes the gene expression variance across the entire set of samples. The 8-gene signature and the disease stage were further used as inputs in a multivariate Cox regression model, to identify whether they had independent prognostic value. This procedure was applied on both the discovery meta-cohort and the TCGA validation data.
Since our procedure for the detection of monotonal traits in genes and pathway activities involved multiple filtering criteria, to avoid over-elimination due to Type-II error, significance was defined at an unadjusted Mann-Whitney p < 0. 05. In contrast, significance for pathway over-representation (clusterProfiler output) was determined by FDR correction at p < 0. 05. Categorical variables were investigated for significance with Pearson's chi-squared test, and were adjusted for multiple hypotheses (with the package RVAide-Memoire). All reported correlation scores corresponded to Spearman's rank coefficient. Cox proportional hazards regression was performed with the packages survminer and survival, and statistical significance was determined via the log-rank method. CIBERSORT analysis was conducted on the web platform https://cibersort.
The results highlighted a set of 12 genes commonly co-expressed across stages for these processes, including COX7B, DLD, NDUFS4, UQCRFS1, PAIP2, RPL15, RPL30, RPL7, RPS23, RPS27, RPS27A, and RPS4X (Table  S10 ). with increasing bladder cancer stage. Further functional enrichment using GeneCards for the 157 genes verified their involvement in cell-cycle pathways, with the top hits being related to the regulation of the anaphase-promoting (APC) complex (score = 31. 53), to PLK1 (score = 24. 47) and Aurora B (score = 20. 95) signaling, and to TP53 (score = 19. 06) and RB1 (score = 17. 85) cell-cycle checkpoint control (Figure  4C , Table  S13 ). Univariate Cox regression analysis indicated 48 genes with potentially prognostic impact at p < 0. 01 (Table  S14 ). Due to the lack of an RNA-Seq dataset comprising the entire disease stage spectrum of BLCA incidents, the observed stage alterations in the discovery set were investigated for their reproducibility in the TCGA-BLCA RNA-Seq data  [41].
Among the 157 genes, 21 were of poorly described or unknown function (Table  S12 ), including the oncogenic factor TRIM65  [40] , which was found to be upregulated with increasing bladder cancer stage. Further functional enrichment using GeneCards for the 157 genes verified their involvement in cell-cycle pathways, with the top hits being related to the regulation of the anaphase-promoting (APC) complex (score = 31. 53), to PLK1 (score = 24. 47) and Aurora B (score = 20. 95) signaling, and to TP53 (score = 19. 06) and RB1 (score = 17. 85) cell-cycle checkpoint control (Figure  4C , Table  S13 ). Univariate Cox regression analysis indicated 48 genes with potentially prognostic impact at p < 0. 01 (Table  S14 ). Due to the lack of an RNA-Seq dataset comprising the entire disease stage spectrum of BLCA incidents, the observed stage alterations in the discovery set were investigated for their reproducibility in the TCGA-BLCA RNA-Seq data  [41].
Along with the disease stage, the eight-gene signature had independent prognostic value, both in the discovery and in the validation set (Figure  5B ). Specifically, we constructed a survival model to compare 5-year survival rates between those with high and low eight-gene signature scores (defined by a median cutoff). Patients with a high eight-gene signature score had a worse 5-year survival probability in the discovery set, and this finding was validated in the TCGA data (Figure  5C ). The eight-gene signature did not differ significantly between males and females (p = 0. 36), and was weakly associated with age and stage (Figure  S7 ), suggesting its independent value with respect to other clinical variables. In an attempt to verify the co-expression analysis findings, stage-specific co-expression networks were also created using the TCGA data, and were clustered using the Louvain algorithm.
GO-Biological Process analysis of the communities validated the differential segregation of the cell cycle, extracellular matrix, and immune activation processes to distinct communities (Figure  S8 ). In order to validate the value of AIF1 as a candidate biomarker of response to immunotherapy, we analyzed RNA-Seq data from the IMvigor210 study-a trial investigating responses to atezolizumab immunotherapy in patients with metastatic BLCA. High AIF1 expression in the IMvigor data was associated with a complete response to atezolizumab (Figure  5D ).
Discussion Integration of BLCA molecular data has been previously performed in the context of characterizing molecular subtypes  [42]  or validating results of either (single-cell) scRNA-Seq  [43]  or RNA-Seq re-analysis  [44]. In this study, we performed an integrated meta-analysis of datasets from non-tumor-bearing adjacent urothelium and BLCA stages, aiming to identify continuous as well as concerted gene expression alterations with increasing malignancy. To our knowledge, this is the first attempt to associate molecular alterations with clinical classification based on the analysis of more than 1000 well-characterized primary tumor datasets. Instead of focusing on molecular subtypes, we increased the power and addressed the disease as a continuum under the assumption that individual samples reflect different snapshots of the whole process.
Our novel design, based on the hypothesis of continuous evolution through the stages, was successfully applied here, and resulted in novel findings on gene regulation associated with cancer progression. Starting from a normalized expression dataset comprising 12 microarray cohorts, we identified genes that were differentially expressed between disease stages and NAU (CDEGs), and further analyzed them for pathways/processes that were progressively altered with higher stage, for changes in the co-expression profiles of stages, and for genes showing a monotonal change in expression with higher stage. Our analysis highlighted expected landmark pathways, such as the mTORC1 pathway [45] and MYC targets  [46] , which were upregulated, but also novel downregulated pathways, such as the circadian clock and the metabolism of heme.
Conclusions The progression of BLCA is associated with monotonal alterations both in individual genes and in molecular pathways, and these alterations can be leveraged as prognosticators, as well as for the development of novel combinatorial drug therapies. The suggested eight-gene signature may help in decision making with regards to therapeutic options and surveillance schedules. However, further validation of the signature in clinical trials is required before clinical implementation. Figure  S8 : Functional annotation of the top 5 largest molecular communities of co-expressed genes detected in the T2, T3, and T4 stages of the TCGA cohort, demonstrating a significant overlap with the cell-cycle, extracellular matrix, and immunity communities of the discovery data; Table  S1 : Processes, pathways and regulons that follow a monotonically increasing/decreasing activation with increasing Bladder Cancer stage (NAU = Non-malignant Adjacent Urothelium); Table  S2 : Gene Ontology Biological Processes for the Non-malignant Adjacent Urothelium (NAU) communities of coexpressed genes; Table  S3 : Gene Ontology Biological Processes for the Ta stage communities of coexpressed genes; Table  S4 : Gene Ontology Biological Processes for the T1 stage communities of coexpressed genes; Table  S5 : Gene Ontology Biological Processes for the T2 stage communities of coexpressed genes; Table  S6 : Gene Ontology Biological Processes for the T3 stage communities of coexpressed genes; Table  S7 : Gene Ontology Biological Processes for the T4 stage communities of coexpressed genes; Table  S8 : Genes of the cell-cycle community and their betweeness centrality scores; Table  S9 : Genes of the ECM and developmental community and their betweeness centrality scores; Table  S10 : Genes of the metabolism and translation community and their betweeness centrality scores; Table  S11 : Genes of the immune community and their betweeness centrality scores; Table  S12 : Functionanl analysis of the 157 monotonically changing genes based on searches in Pubmed and GeneCards; Table  S13 : Functional analysis results for implicated pathways of the 157 monotonically changing genes using the GeneCards tool and database; Table  S14 : Univariate cox regression analysis of the 157 monotonal genes for the discovery and the TCGA validation cohorts.
In green are highlighted genes who were validated with cox in the TCGA data; Table  S15 : Results of differential expression analysis (Mann-Whitney) between bladder cancer stages using only primary samples (70 T2, 90 T3, and 28 T4) from the TCGA Bladder Cancer 2017 cohort. Blue and red color indicate downregulation and upregulation, respectively; File S1: ComBat_expression_data; File S2: Clinical_data. Appendix A. 4 Integration and Processing of the Transcriptomes Pretreatment, raw microarray profiles from the 12 GEO datasets were summarized to the gene level with the package oligo  [74]. Affymetrix arrays were normalized via the RMA method, while Illumina arrays were filtered based on detection p values < 0. 05, followed by quantile normalization, the addition of 1, and transformation to the natural logarithmic scale.
Data were annotated using biomaRt (v2. 42. 1)  [75] , and microarray probes matching to multiple genes were excluded from downstream analysis. The probe with the highest mean across arrays was selected as representative in cases where multiple probes matched to the same gene. To completely avoid missing values, merging of expression matrices was based on the intersection of genes between datasets, using the HUGO Gene Symbol as a key. Batch effects across different studies in the discovery data were removed with ComBat  [76] , and the quality of the corrected data was assessed with BatchQC  [21]  (Figures  S1  and S2 ), with boxplots of expression distribution per sample and principal component analysis plots of sample relationships (Figure  S3 ), with gene expression comparisons of housekeeping genes to other genes (Figure  S4 )  [77] , and with a set of 12 positive control BLCA markers with known regulation across normal-NMI-MI or across normal-low-grade-high-grade disease (Figure  S5 ).
[70]3586[70]: Sjödahl G, et al. , 2017. The Journal of Pathology: This study utilizes tumor samples from 307 advanced bladder cancer patients in order to compare molecular classification to tumor cell phenotype (the latter assessed via immunohistochemistry). We excluded six samples due to the absence of stage information. • GSE48276 [71]: Choi W, et al. , 2014. Cancer Cell: The MDA molecular classification of muscle-invasive bladder cancer. The GEO dataset contains 116 expression profiles, 16 of which were collected post-treatment and, thus, excluded. An additional set of 20 profiles had no information on the timing of sample collection, and were also excluded. • GSE52219 [71]: Choi W, et al. , 2014. Cancer Cell: The MDA molecular classification of muscle-invasive bladder cancer. The GEO dataset contains 23 FFPE pretreatment tumors from a phase III clinical trial, which were used to validate the MDA classification.
• GSE69795 [72]: McConkey, et al. , 2016. European Urology: Investigation of gene expression signatures predicting response to methotrexate, vinblastine, doxorubicin, and cisplatin with bevacizumab. The GEO dataset contains 61 profiles, of which 23 and 1 were excluded due to post-treatment collection and the absence of stage information, respectively. • GSE13507 [73]: Kim WJ, et al. , 2010. Journal of Clinical Oncology: Analysis of primary bladder cancer for the identification of a prognostic gene expression classifier. The GEO dataset contains 256 profiles, of which 9 are normal bladders, 58 are NAU, 1 profile is taken from a C57BL/6 mouse sub-strain, 165 are primary tumor samples, and 23 correspond to recurrences. Due to the low number of normal bladder samples represented in the compiled discovery cohort, we grouped them together with the NAU samples.
Samples collected after recurrence (n = 23), as well as the expression profile of the murine strain, were excluded from analysis. • GSE48075 [71]: Choi W, et al. , 2014. Cancer Cell: The MDA molecular classification of muscle-invasive bladder cancer. It includes 142 primary tumors of both non-muscle-invasive and muscle-invasive bladder cancer. Validation (RNA-Seq, Datasets 2, Samples = 318) • TCGA-BLCA-2017 [42]: Robertson AG, et al. , 2017. Cell: The TCGA molecular subtyping of muscle-invasive bladder cancer. Molecular and clinical data were downloaded from cBioPortal, containing 412 tumor samples. Of these, 109 had history of other malignancies, and were excluded as potentially being secondary bladder cancers. Additionally, 105 samples had unknown history of prior treatment for non-muscle-invasive bladder cancer, 2 samples were collected post-NAC, and 8 more samples had no stage information, so they were all excluded from further analysis (n analyzed = 186 samples).
Funding:  Marika Mokou  is supported by the  European Union 's  Horizon-2020 Marie Sklodowska-Curie Research and Innovation Program   H2020-MSCA-IF-2019 ,  ReDrugBC  (Grant agreement ID:  898260 ). Emmanouil Mavrogeorgis is supported by the  European Union's Horizon-2020 Marie Sklodowska-Curie Research and Innovation Program, STRATEGY-CKD  (Grant agreement ID:  860329 ).
Data Availability Statement: This study used publicly available data. The raw microarrays used here are publicly available in the Gene Expression Omnibus under the accession codes GSE121711, GSE93527, GSE31684, GSE104922, GSE128959, GSE83586, GSE48276, GSE52219, GSE69795, GSE13507, and GSE48075, and in ArrayExpress under the accession code E-MTAB-1940. TCGA-BLCA-2017 molecular and clinical data were downloaded from cBioPortal. IMvigor molecular and clinical data were acquired through the IMvigor210CoreBiologies R package. The processed, ComBat-normalized matrix (discovery set) and the sample clinical data are available as Supplements. The codes used to generate the results and figures are available upon reasonable request.
• [70]3586[70]: Sjödahl G, et al. , 2017. The Journal of Pathology: This study utilizes tumor samples from 307 advanced bladder cancer patients in order to compare molecular classification to tumor cell phenotype (the latter assessed via immunohistochemistry). We excluded six samples due to the absence of stage information. • GSE48276 [71]: Choi W, et al. , 2014. Cancer Cell: The MDA molecular classification of muscle-invasive bladder cancer. The GEO dataset contains 116 expression profiles, 16 of which were collected post-treatment and, thus, excluded. An additional set of 20 profiles had no information on the timing of sample collection, and were also excluded. • GSE52219 [71]: Choi W, et al. , 2014. Cancer Cell: The MDA molecular classification of muscle-invasive bladder cancer. The GEO dataset contains 23 FFPE pretreatment tumors from a phase III clinical trial, which were used to validate the MDA classification.
Randomized reconstruction ranking To quickly review and rank thousands of neuronal reconstructions by the experimenter in a double-blinded, randomized manner across all samples in an experimental dataset, we built the last App in the PopRec pipeline-Reconstruction Validation App. This tool presents the neuron reconstruction with the original image in a one-by-one manner, both in the interactive 3D format. This review process can be run multiple times on the same neurons and when completed, it ranks the accuracy of each reconstruction by the number of times that neuron was reviewed  (Fig 1E, see Appendix Fig S4 for more information). We used this ranking to train a MATLAB classifier based on the morphological features and ranks for all of the neuron reconstructions (including multiple versions of each neuron), which allowed to both "clean out" poorly reconstructed neurons and also select the best versions of high-accuracy reconstructions  (Fig 1F).
We found full-length accuracy in 71. 7% (228/318), fine structure errors in 22. 3% (71/318), and significant errors in 6% (19/318) of the analyzed cells. We defined SD as the full dendritic trees that reside from the soma: a bipolar neuron would have two SD and multipolar more than two. Overall, the PopRec approach hence yielded 94% accuracy in terms of basic neuron morphology and 72% accuracy in terms of the fine dendritic details (Dataset EV3, Movie EV2).