Integrative analysis of genomic and epigenomic regulation reveals microRNA regulatory network mediated tumor heterogeneity and immune evasion in lower grade glioma

Integrative analysis of genomic and epigenomic regulation reveals microRNA regulatory network mediated tumor heterogeneity and immune evasion in lower grade glioma | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Integrative analysis of genomic and epigenomic regulation reveals microRNA regulatory network mediated tumor heterogeneity and immune evasion in lower grade glioma Zhen Yang, Xiaocen Liu, Hao Xu, Andrew E. Teschendorff, Lingjie Xu, and 9 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-3935802/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Background Lower grade glioma (LGG) is the most frequent primary tumors of the central nervous system and has been a major healthcare burden, however, the specific molecular mechanism underlying its initiation and progression remains to be elucidated. Although it is known that microRNAs (miRNAs) are widely involved in the regulation of malignant phenotypes of glioma, the underling mechanism for miRNA dysregulation remains largely unanswered. Methods In the present work, we developed a novel strategy to obtain the genome wide copy number variation (CNV) and promoter DNA methylation (DNAm) data of miRNAs and performed a systematic integrative study for the multi-omics data to identify mechanisms underlying miRNA dysregulation molecular subtyping in LGG. The relationship between LGG subtypes, prognosis, molecular features, tumor immune microenvironment and response to immune therapy was further analyzed. We also developed a prognostic model based on immune-related miRNAs that were differentially expressed between LGG samples. Then, the influence of the prognostic model on the immune microenvironment in LGG was comprehensively analyzed. Results We identify 719 miRNAs whose expression was associated with alterations of copy number variation or promoter methylation. Integrative multi-omics analysis of the copy number and methylation related miRNAs revealed four subtypes with differing prognoses, which were validated with independent cohort data. These glioma subtypes exhibited distinct immune-related characteristics as well as clinical and genetic features. We further screened immune-related miRNAs through investigation of their correlation with immune cell infiltrations and immune microenvironment. By construction of a miRNA regulatory network, we identified candidate miRNAs associated with immune evasion and response to glioma immunotherapy. We finally evaluated the associations between prognosis related miRNAs and immune regulation. Among them, miR-155-5p, miR-196a-5p, miR-196b-5p, miR-200a-3p, miR-503-5p, and miR-15b-5p were validated as immunoevasive biomarkers and to promote cell migration, invasion and proliferation for glioma through in vitro experiments. Conclusions Our study systematically reveals the crosstalk among DNA methylation, copy number variation and miRNA expression for immune regulation in glioma, and could have important implications for patient stratification and development of novel biomarkers for immunotherapy approaches. Lower grade glioma microRNA DNA methylation Copy number variation Omics-data integration Tumor immune microenvironment Regulatory network Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Figure 9 Figure 10 Background Lower grade glioma (LGG) is the most common and aggressive central nervous system (CNS) tumor, which accounting for about 40% of all brain malignancies ( 1 ). Most LGGs can be further classified according to specific type of cell with which they share histological features or classic molecular markers, such as isocitrate dehydrogenase (IDH) mutation or O 6 -methylguanine-DNA methyltransferase (MGMT) promoter methylation ( 2 , 3 ). Despite tremendous progress has been made over the past decades for its early detection, surgical paradigm and multidisciplinary treatment, such as neoadjuvant chemotherapy and radiotherapy, the postoperative overall survival of LGG patients remains extremely low, particularly for high-grade glioma with worse prognoses for their malignant aggressivity ( 4 ). In general, the prognosis of patients with glioma varies dramatically, which is dependent on different tumor grades, key gene alternation status, tumor microenvironment (TME), and combination of different efficacious method ( 5 ). Recent genomic studies have brought a comprehensive view of glioma based on molecular profiling and identified different molecular subtypes within the glioma that appear to correlate with disease etiology and therapy response ( 6 ). Hence, there is an urgent need to characterize more specific and precise molecular signatures of LGG for its accurate diagnosis, individualized treatment and the prognosis assessment. MicroRNAs (miRNAs) are ~ 22 nt small noncoding RNAs that play an important role in post-transcriptional regulation, which precisely adjust the gene expression level by targeting mRNAs for its degradation or translational repression ( 7 , 8 ). MiRNAs are critical for normal development of organisms and are involved in a variety of biological processes including cell cycle, cell proliferation, differentiation, apoptosis and cellular signaling ( 9 , 10 ). Aberrant expression of miRNAs is associated with many human diseases ( 11 , 12 ). For glioma, miRNA expression has been associated with the clinical and molecular classification ( 13 ), progression ( 14 , 15 ), prognosis ( 16 , 17 ) and response to targeted therapies ( 18 – 20 ). Altered miRNA regulation is widely involved in glioma pathogenesis via the modulation of oncogenes and the associated downstream signaling pathways. Till now, most studies were proposed concerning aberrant expression and functional curation of miRNAs in cancer and disease, but there still lacks systematic investigation of the mechanisms underlying miRNA regulatory activities, especially the genome-wide studies for genetic and epigenetic regulation of miRNAs in development and progression of glioma are still scarce. In recent years, large multi-omics studies greatly enhanced our understanding of disease dysregulation at genomic and epigenetic levels. Genome wide alternations such as copy number variations (CNVs) and DNA methylations (DNAm) are widely observed during tumorigenesis, which promote cancer progression ( 21 , 22 ). Copy number variation refers to the number of copies of a specific DNA segment in the genome varies among individuals ( 23 ). As a key genomic regulator that contribute to gene expression dysregulation via modulating mRNA levels and by influencing transcriptional regulation, CNV is significantly correlated with certain pathological processes such as cancer ( 24 ). Transcriptional disorders caused by CNV changes have been identified as driving events in LGG progression. In addition, gene expression regulation via DNA methylation at the epigenetic level also exerts a crucial part in the behaviors of various cancers ( 25 , 26 ). DNAm imbalances is able to affect the occurrence and development of cancer by defining different types of driver events, such as cell growth, proliferation, differentiation, and apoptosis processes ( 27 , 28 ). Genomic profiling studies indicated that these genomic and epigenomic dysregulations are highly heterogeneous and jointly regulate the occurrence and evolution of tumor ( 29 ). Till now, the influence of DNAm and CNV on noncoding RNA, especially miRNA expression remains largely unexplored. Studies indicated various regulators are required for precise control the expressions of miRNAs ( 30 ). Therefore, miRNAs with aberrant expression under genetic and epigenetic contexts may serve as valuable marker for investigation of the complex regulatory mechanisms in tumor progression, and a systematic analysis of DNAm and CNV patterns of miRNAs in large scale cancer samples is critically needed. In this study, we collected copy number variation, DNA methylation and miRNA expression profiles from the cohort of The Cancer Genome Atlas (TCGA) LGG patients, and identified miRNAs whose expression levels are regulated by genomic and/or epigenomic deregulation. In addition, by using the multi-omics integration study, we identified different molecular subtypes of LGG that are significantly associated with prognosis. Furthermore, specific miRNA biomarkers were proposed to drive the classification of these subtypes, which were validated in independent cohort data base on different platforms. By assessing the associations between immune cell infiltrations and immune pathways, we identified a subset of immune-related miRNAs, and several miRNAs are identified as immune evasion markers with high confidence. Our study developed novel strategies of multi-omics integration and characterized the genetic and epigenetic landscape of genes encoding miRNAs, and also provides a framework for identifying novel biomarkers to guide the clinical treatment for glioma patients. Material and Methods Data collection and preprocessing The genome wide profiling data for LGG was download from TCGA ( https://portal.gdc.cancer.gov/ ) ( 31 ). Here we obtained 512 samples with miRNA expression, 505 samples with DNA methylation and 515 samples with CNV detection, a total of 500 samples with all three different types of high-throughput profiling data were retained for further integrative study. The genome-wide miRNA profiling data which based on miRNA-seq was quantified as RPM (reads per million mapped reads). The expression value was log2(RPM + 1) transformed in order to regularize the data. For DNA methylation data based on Illumina Infinium HumanMethylation450 BeadChip array, methylation level of each probe was measured as the beta-value. We removed the probes whose beta-values are missing in more than 30% of the samples, the remaining probes with missing values were imputed by using the k-nearest neighbors (KNN) method ( 32 ). Finally, the BMIQ method was used to correct for the type II probe bias ( 33 ). For CNV profile, segmentation data generated by Affymetrix SNP 6.0 platform were used. The ‘nocnv.seg’ file contains chromosome location and segment mean value (log2(copy-number/2)) for each sample and was used to capture the focal somatic CNV. In addition, we also collected RNA-seq data for 500 LGG samples which was quantified as logarithm transformed (base 2) FPKM (Fragments Per Kilobase per Million), as well as the somatic mutation data for 508 samples generated by MuTect2 pipeline. Furthermore, clinical information of all samples including age, gender, survival status, histological type, tumor stage and overall survival time was also retrieved from the TCGA data portal. An external dataset of 198 glioma patients with miRNA expression profiling data was downloaded from the Chinese Glioma Genome Atlas (CGGA) ( 34 ). This miRNA profiling data is based on the Illumina Human v2 MicroRNA Expression BeadChip. We transfer the miRNA ID from the microarray data to that from miRBase v21 by using the “ID convertor” tool in dbDEMC v3.0 ( 35 ). The expression values were logarithmically transformed (base 2) and quantile normalized. The associated clinical data which includes survival status and overall survival time of the patients was also obtained. Patients in the in-house cohort and sample collection Human tissue samples of 30 pairs of adjacent non-tumor and glioma samples were collected from the Department of Neuro-surgery of the First Affiliated Hospital, Wannan Medical College (Wuhu, Anhui, P. R. China) and Department of Neuro-surgery of Huashan Hospital, Fudan University (Shanghai, P. R. China) from February 2023 to September 2023. Patients received no systemic treatment before the surgery. All samples obtained at the surgery were directly preserved in liquid nitrogen. Two pathologists evaluated all specimens according to the WHO guidelines. Patient consent was obtained. All surveys and experiments were approved by the Ethic Committee for Clinical Research of the First Affiliated Hospital of Wannan Medical College and Huashan Hospital, Fudan University, respectively. DNA copy numbers profiling construction for miRNAs To construct copy number profiling for miRNAs in LGG, we investigated whether the positions of miRNA loci can be mapped to genomic regions with copy number in a certain sample. We first obtained the genomic coordinates for mature miRNAs via miRBase version 21 ( 36 ). Then segment-mean values in distinct genomic loci (genome assembly hg38) were extracted from TCGA SNP array data. We assigned the locus of each miRNA to the copy number values by using the GISTIC2 ( 37 ). In this way, the copy number profiling data for individual estimates were obtained for a total of 2265 miRNAs. miRNA promoter identification and Infinium 450K Array probe reannotation to construct miRNA methylation profiles To characterize DNA methylation patterns for miRNAs, we employed a novel strategy to reannotate the probes of Infinium 450K arrays into miRNA-associated promoter regions. We first obtained genome coordinates for miRNA promoters from the fifth edition of the Functional Annotation of Mammalian Genome database (FANTOM5) ( 38 ). FANTOM5 incorporates miRNA transcription start site (TSS) information generated by Cap Analysis Gene Expression (CAGE) methods, which was based on genome assembly hg19. The coordinates of miRNA TSS were converted to genome assembly hg38 by using the liftOver tool from the UCSC genome browser. As a result, a total of 2050 mature miRNAs are identified for unique physical genomic location of TSS among 2587 miRNAs from miRBase v21 (Table S1 ). Here we define the promoter region for individual miRNA as between 2000 bp upstream to 500 bp downstream of TSS, which is akin to protein coding genes. By matching the genomic coordinates of the miRNA TSS with the 450K Array probe locations, we filtered 8825 cpg sites that are located within miRNA promoters (Table S2 ). Then we assign the average methylation value of the probes within promoter region for each miRNA to obtain its final methylation value. Finally, we obtained methylation profile for 1977 miRNAs for LGG. Identification of CNV and DNAm related miRNA sets The Pearson correlation coefficients (r) were calculated for paired profiles of CNV and miRNA-seq and between DNAm and miRNA-seq for each miRNA, respectively, the correlation coefficient was then transformed to the Fisher Z-scale according to the formula: $$FisherZ=\frac{1}{2}ln\frac{1+r}{1-r}$$ The miRNA of FDR adjusted p-value < 0.05 with correlation test constitutes the miRNA set significantly positively related to CNV (CNV-miRs) and the miRNA set negatively related to methylation (DNAm-miRs), respectively. Clustering analysis of multi-layered genomic profiles Nonnegative matrix factorization (NMF) is an unsupervised clustering method which was widely used in tumor subtypes identification based on omics data. Here we employ NMF to identify stable sample clusters by using expression of CNV-miRs and DNAm-miRs, respectively. Specifically, cluster analysis with standard “brunet” method and 50 iterations was used, we set the number of cluster k as 2 to 7, the optimal cluster number k was then determined based on the observed consensus map together with cophenetic correlation between clusters. At the same time, the mean silhouette width was computed by R package “NMF” to examine consensus membership matrix ( 39 ). For integrative cluster analysis of the multi-omics data for cancer subtyping, we applied “MOVICS” R package to integrate the CNV data of the CNV-miRs, methylation data of the DNAm-miRs, as well as the expression profile data of both CNV-miRs and DNAm-miRs ( 40 ). MOVICS provides an interface for 10 state-of-the-art multi-omics data clustering algorithms (SNF, PINSPlus, NEMO, COCA, LRAcluster, CIMLR, Consensus clustering, IntNMF, MoCluster, iClusterBayes) and combines the output of each algorithm to obtain more robust classification. The performance of these algorithms had been validated by previous researches ( 41 ). We calculated the clustering prediction index to find the optimal number of clusters and the consensus matrix was used to validate robust clustering of the samples. Measurement for immune cell infiltration and antitumor immunoactivity In this study, we identified and evaluated the levels of immune infiltrates by the Tumor immune estimation resource (TIMER) algorithm, which is a resource for systematic estimation of the abundances of six immune cell types (B cell, CD4 + T cell, CD8 + T cell, neutrophil, macrophage and dendritic cell) by using transcriptome profiles ( 42 ). The abundances of these cells were compared across different molecular subtypes for quantitative characterization of the tumor immune microenvironment. We used ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) to assess immune activity of the tumor samples. ESTIMATE yields the immune score and stromal score for each tumor sample that quantifies the levels of infiltration levels of immune and stromal cell based on its gene expression profiles ( 43 ). We collected 17 immune relevant pathways and the associated 1,811 genes from ImmPort ( 44 ). In addition, three different indices were constructed to measure the immune activities for tumor samples: Cytolytic activity (CYT) score, Major Histocompatibility Complex (MHC) score and Cytotoxic T Lymphocyte (CTL) score. The CYT score was defined as the geometric mean of GZMA and PRF1 : $${Score}_{CYT}=\sqrt{{Exp}_{GZMA}{*Exp}_{PRF1}}$$ These two genes are found upregulated in activated CD8 + T cells and play important roles in cytolysis and response to PD-1 and CTLA4 related immunotherapy ( 45 , 46 ). The MHC score of each sample was calculated as the average expression as nine genes of HLA-A , PSMB9 , HLA-B , PSMB8 , HLA-C , B2M , TAP2 , NLRC5 and TAP1 : $${Score}_{MHC}=\sum {Exp}_{n}/9$$ These genes were identified as the core gene set of MHC-I and are associated with activity of MHC-I antigen presentation ( 47 ). Finally, the CTL score for each sample was calculated as average expression from five genes of CD8A , CD8B , GZMA , GAMB and PRF1 : $${Score}_{CTL}=({{Exp}_{CD8A}+{Exp}_{CD8B}+Exp}_{GZMA}+{Exp}_{GZMB}+{Exp}_{PRF1})/5$$ These five genes were important factors to measure cytotoxicity of tumor-infiltrating T cells and immune cell effector function ( 48 ). In addition, we also obtained the gene sets that represents different immune signatures from several publications, including the immune check point, human leukocyte antigen (HLA), interferon (IFN) response and tumor-infiltrating lymphocytes (TILs), which are collected from publication of Ju et al. ( 49 ), and the gene sets for Macrophage/Monocyte, TGF-β response, IFN-γ response, and Wound healing collected from Vesteinn et al . ( 50 ). We use the single-sample Gene Set Enrichment Analysis (ssGSEA) algorithm to estimate the activity of tumor samples by the immune signatures obtained ( 51 ). Finally, another 63 immune cell type specific gene markers for 21 immune cell populations were collected from TISIDB ( 52 ). Immune evasion and immunotherapy response evaluation To evaluate complex interactions among cancer cells and immune microenvironment, we constructed the scoring system following Shuichi et al . to describe factors that relevant to the tumor proliferation (T-factor), antitumor immunity (I-factor) and immunosuppression (S-factor) ( 53 ). Briefly, nine representative gene sets related to those factors were collected from the molecular signature database (MsigDB V.7.1) ( 54 ), which include Gene Ontology (GO) terms of “SPINDLE_LOCALIZATION” and “POSITIVE_REGULATION_OF_GLUCOSE_METABOLIC_PROCESS” for T-factor, “T cell”, “B cell” Signatures and GO terms of “MHC_PROTEIN_COMPLEX” for “I-factor”, and “NEGATIVE_REGULATION_OF_CYTOKINE_PRODUCTION_INVOLVED_IN_INFLAMMATORY_RESPONSE”, “POSITIVE_REGULATION_OF_CELL_MIGRATION_INVOLVED_IN_SPROUTING_ANGIOGENESIS”, “NEGATIVE_REGULATION_OF_HYPOXIA_INDUCED_INTRINSIC_APOPTOTIC_SIGNALING_PATHWAY” as well as “EPITHELIAL_TO_MESENCHYMAL_TRANSITION” for “S-factor”, the individual ssGSEA enrichment score was calculated and z-score normalized, the mean values of z-scores for the gene sets allocated to T-factors, I-factors or S-factors and scored was calculated to obtain T-scores, I-scores or S-scores. Finally, we apply the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm to estimate the immunotherapy response of LGG patients ( 55 ). Higher TIDE scores implicated a higher potential of tumor immune evasion, thus less likely to benefit from immunotherapy. Immune-related miRNAs identification and prognosis model development To identify immune regulatory miRNAs involved in LGG subtyping, Pearson correlation analysis was performed between miRNAs expression and levels of immune cell infiltrating. MiRNAs with FDR less than 0.05 were considered as significant. We further collected the immune related miRNAs from the RNA2Immune database for functional validation ( 56 ). This database provides a high-quality experimentally supported resource that links ncRNA regulatory mechanisms to four sections of “immune cell function”, “immune diseases”, “cancer immunity”, and vaccines. We selected miRNAs from the “immune cell function”, and “cancer immunity” for validation propose. To find the miRNAs with the height prognostic value, we applied Cox-proportional hazards analysis based on the Least absolute shrinkage and selection operator (LASSO) estimation, which have been identified that suitable for model construction when dealing with collinearity problem of covariates. Among the immune miRNAs that were significant correlated with immune cell infiltrating, key immune miRNAs were selected by performing the LASSO-penalized Cox regression, which was implemented by the glmnet R package (version 4.1-7). Finally, an immune miRNA-related prognostic model was constructed by utilizing the regression coefficients derived from Cox regression analysis to multiply the normalized expression level of each immune miRNA. The “surv_cutpoint” function within the “survminer” package was applied to determine the best cutoff to classify patients into high-risk and low-risk groups. The Kaplan-Meier survival analysis and log-rank test were used to assess the predictive ability of the prognostic model. Validation of the prognostic model from independent cohorts To validate whether the predictions of the prognostic model were independent of other clinical features, we collected miRNA profiling data form CGGA and calculated the risk scores for the patients by using the model constructed. In addition, the predictive accuracies were compared using the time-dependent receiver operating characteristic (ROC) analyses and the associated Area Under Curve (AUC) Score. Construction and evaluation of the nomogram We assessed the correlation between prognostic model and other clinicopathological features. The Cox proportional hazards model was used to perform standard univariate and multivariate analysis. Prediction error curves were used to compare the accuracy of survival models. The Cox regression coefficients were used to construct the nomogram. Calibration plots were generated to explore the performance characteristics of the nomogram. Time dependent ROC analysis was performed to assess the predictive accuracy of the nomogram. Decision curve analysis was used to assess the clinical practicability of the nomogram. The statistical significance was set at 0.05. miRNA target prediction and functional enrichment analysis For both miRNAs and protein coding genes that differentially expressed compared between groups, the well-known empirical Bayes framework of moderated t-statistics were employed for assessing differential expression in microarray experiments ( 57 ). miRNA target prediction was obtained from databases of TargetScan ( 58 ) and miRDB ( 59 ). The miRNA-target interactions predicted by both of the two algorithms were identified as valid ones. In addition, we also collected experimentally validated miRNA targets form miRTarBase ( 60 ). For functional annotation of the miRNA targets involved in the regulatory network, The GO enrichment analysis was performed by using clusterProfiler ( 61 ), the enriched GO terms were obtained at the threshold of FDR < 0.05. In addition, hallmark gene sets were obtained from MsigDB, and Gene Set Variation Analysis (GSVA) analysis was used to infer the degree of activity of gene sets/pathways in individual sample Distinct genomic features across different groups The tumor mutation burden (TMB) and oncoplot of somatic mutation data were calculated and generated by using R package “maftools” ( 62 ). The different mutational frequencies for individual genes between groups were calculated by chi-square test, the p-value that less than 0.05 were identified as having achieved statistical significance. In addition, the druggable pathways and categories are also obtained by enrichment analysis. Quantitative real-time polymerase chain reaction (qRT–PCR) Total RNA was extracted from tissues with TRIzol reagent from 30 pairs of LGG and normal tissues according to the manufacturer’s protocol. Reverse transcription was performed using the PrimeScript RT Reagent Kit. The quantitative PCR assay was performed using the SYBR Green PCR kit. The expression of eight miRNAs of the in-house cohort was normalized by GAPDH, acting as an internal control. Relative miRNA expression levels were calculated by the 2 −ΔΔCT method. Each sample was repeated in triplicate and analyzed using relative quantification method. Cell culture and transfection The human glioma cell line U251 used in this study was purchased from the American Type Culture Collection (ATCC, Manassas, VA, USA). DNA fingerprinting, cell vitality detection, isozyme detection, and mycoplasma detection were used to characterize all cell lines. The cell line was regularly cultured in Dulbecco's modified Eagle media (DMEM; Gibco; Thermo Fisher Scientific; Waltham, MA, USA) supplemented with 10% fetal bovine serum (FBS; Gibco; Thermo Fisher Scientific) at 37°C in a humid environment with 5% CO 2 . Transfection of miRNA inhibitors were performed using Hieff TransTM in vitro siRNA/miRNA Transfection Reagent (Yeasen, China) according to the manufacturer’s instructions. U251 cells were plated in T-25 cells culture flasks and 96 wells plates, 24 wells plates, and transfected with miRNA inhibitors. At 48–72 h after transfection, the cells were collected and used for experiments. Western blotting The membrane was incubated with the following primary antibodies (1:1000) at 4°C overnight after being blocked with 5% non-fat milk at room temperature for 1 hour: anti-PD-L1 (Cat#13684), anti-β-actin (Cat#4970), anti-CTLA4 (Cat#53560), and anti-FOXP3 (Cat#19713) antibodies. The membrane was washed three times, followed by treatment with horseradish peroxidase-conjugated anti-mouse or anti-rabbit secondary antibody. After that, a Tanon 5200 system (Tanon, Shanghai, PR China) was used to photograph and wash the membrane. Using Quantity One software from Tanon (Shanghai, PR China), the intensity of protein bands was assessed using densitometric analysis. Target gene protein expression was calibrated to β-actin protein expression. A triplicate of the experiment was run. CCK-8 assay Cells were seeded, cultivated for an overnight period, and then given the CCK8 solution in 96-well plates. The cells were then treated with the CCK-8 reagent at 37°C in a humid environment that contained 5% CO 2 . The absorbance at 450 nm was measured on an enzyme analyzer (Tecan infinite M2009PR, Tecan, Männedorf, Switzerland). The experiment was carried out three times. Transwell assay Cell migration and invasion were performed with Transwell chambers (BD Biosciences, USA). After 48 h of corresponding treatment, U251 cells in each group were inoculated into the upper chamber of the Transwell chamber at the density of 4×10 4 cells/well (100 µL culture medium containing 5% fetal bovine serum). Besides, 500 µL culture medium containing 10% fetal bovine serum was added to the 24-well culture plate of the lower chamber. After 24 h of routine culture, the chamber was removed, and cells from the upper layer of microporous membranes were wiped with cotton swabs. After that, cells were immobilized in 4% paraformaldehyde solution for 10 minutes at room temperature and stained with 0.5% crystal violet solution (Sigma-Aldrich; Merck KGaA) for 15 minutes. In the last step, 5 visual fields were randomly selected and observed under an optical microscope (Nikon, Japan) to count the number of cells invading the sub-layer of the microporous membranes of the chamber. Wound healing assay Cell migration was observed using a wound-healing assay. When the transfected U251 cells were maintained in a 6-well plate, achieving 90–95% confluence, scratches were generated using the micropipette tips. The wound state was observed at 0 h and 24 h after scratching with an X71 inverted microscope (Olympus, Tokyo, Japan). Statistical analysis Unless stated otherwise, all statistical tests were performed by R software (version:4.1.3) The significance of differences between two groups was determined by Wilcoxon rank sum test. Hypergeometric test was used to determine the significance of overlap between two groups of genes. For patient survival analyses, Kaplan-Meier plots were created, we used the Cox proportional hazard model and the log-rank test to determine the difference of survival between different groups of the patients. Unless stated otherwise, a p-value < 0.05 was considered significant. If necessary, p-values were corrected for multiple tests with the Benjamini-Hochberg procedure to calculate false discovery rate (FDR). Results MiRNA deregulation mediated by copy number variation and DNA methylation Genomic and epigenomic profiling data of copy number variation, DNA methylation, and miRNA expression were obtained for 500 LGG patients from TCGA. We first explored the global CNV characteristics of LGG on the whole genome. Based on the copy number amplitude calculated by GISTIC2, we found that there is a concentrated copy number amplification on q32.1 of chromosome 1, p11.2 of chromosome 7 and 12q.14 of chromosome 11 and copy number deletion on q37.3 of chromosome 2, p21.3 of chromosome 3, q23.31 of chromosome 10, and q13.2 of chromosome 9 (Fig. S1 ). The copy number amplitude data was then mapped to miRNA loci to build copy number profiles for 2265 miRNAs. Next, by mapping the 8825 450K array probes to miRNA promoter regions, we obtained DNA methylation profiling data for 1977 miRNAs (Table S1 and Table S2 ). The metaplot was used to depict the methylation level around miRNA transcription start sites. In general, the methylation level around transcription start sites was significantly lower than that from gene body and intergenic regions for miRNAs (Fig. S2 ). To further investigate the coordination between CNV and DNAm abnormalities, we also compared the frequencies of the miRNA associated aberrant CNV and DNAm in the genome. We define the CNV value > 0.3 as CNV amplification (CNVamp); CNV value 0.8 was defined as hypermethylation (HyperM); β-value < 0.2 was defined as hypomethylation (HypoM). The frequencies of CNVamp, CNVdel, HyperM and HypoM for each sample were calculated. We found that the frequencies of CNVamp are mildly correlated with the frequencies of CNVdel (Fig. S3A, r = 0.1, p = 0.03), whereas the frequencies of HyperM are significantly negatively correlated with the frequencies of HypoM (Fig. S3B, r = -0.7, p < 2.2e − 16). In addition, the frequencies of all the aberrant copy number are significantly correlated with the that of all aberrant DNA methylation (Fig. S3C, r = 0.11, p = 0.015). Abnormalities in directional CNVamp, CNVdel, HyperM and HypoM are also tightly correlated (Fig. S3D-G). In summary, our observation suggests that the LGG patients with frequent aberration of miRNA copy number are more likely to have frequent aberration of promoter methylation. These correlated frequencies of aberrant CNV and DNAm highlights the widespread involvement of genetic and epigenetic regulations for miRNAs in the pathogenesis of LGG. To assess the global effects of genomic and/or epigenomic aberrations on miRNA expression, the correlation coefficients between CNV or DNAm with corresponding miRNA expression were calculated and then normalized by using Fisher’s Z-transformation. It was found that the overall correlation coefficients between CNV and the corresponding miRNA expression shifted to the right (skewness = 2.43, p < 2.2e − 16, D’Agostino test), whereas that between DNAm and the miRNA expression significantly shifted to the left (skewness = -1.11, p < 2.2e − 16, D’Agostino test) (Fig. 1 A). This indicated the overall positive and negative regulatory effect on miRNA expression by CNV and DNAm, respectively, which is similar to that of protein coding genes. We next focus on the positively correlated miRNAs for CNV (CNV-miR) and the negatively correlated miRNAs for DNAm (DNAm-miR). A total of 351 CNV-miRs (Table S3) and 541 DNAm-miRs (Table S4) were identified at the thresholds of FDR < 0.05. Genome-wide landscape for the associations between CNV and DNAm for the two different groups of miRNAs are depicted as Fig. S4. Some top ranked CNV-miRs and DNAm-miRs with significant associations between the expression with CNV and DNAm are presented as Fig. S5 and Fig. S6, respectively. Many of these miRNAs have been identified that play a key role in different cancer types. One of the examples is the hsa-miR-26a, whose expression is highly correlated with CNV value (r = 0.62, p < 2.2e − 16). This miRNA has been found as a key driver gene for glioblastoma, its copy number amplification associated overexpression was observed in glioblastoma to transform cells and to promotes tumor growth ( 63 , 64 ). Another example is miR-155-5p for its negative correlated pattern between expression and promoter methylation (r = -0.74, p < 2.2e − 16). This negative correlation has been observed in several other studies with regard to glioma, which was indicated that significantly correlated with patients’ survival ( 65 , 66 ). Similar phenomenon has been observed in many other cancer types including breast cancer ( 67 ), gastric cancer ( 68 ) and T-cell lymphoma ( 69 ), which suggests miR-155-5p act as a pan-cancer regulator with functions regardless of tumor origin. We next inspected the genome distribution of CNV-miRs and DNAm-miRs and found that both groups of miRNAs are mainly distributed on chromosome 14, followed by chromosome 1 and chromosome 19 (Fig. 1 B). The CNV-miRs and DNAm-miRs present an overlap of 173 miRNAs, which suggests that a significant number of miRNAs are coordinated regulated by CNV and DNAm, whereas others are mutually exclusive regulated (Fig. 1 C). Additionally, we also checked the positional distribution of the probes within miRNA genes and found that DNAm-miRs probes are more frequently located in promoter regions of 2000 bp upstream from the TSS to 200 bp upstream from the TSS (TSS2000) rather than within 200 bp upstream of the transcription start site (TSS200) and gene body (Fig. 1 D). Moreover, DNAm-miRs probes are more frequently located in N-shores rather than in S-Shores, OpenSea and CpG island regions (Fig. 1 E). We further investigated the functions of the miRNA by performing the Gene Ontology enrichment analysis based on the target genes of the two different groups of miRNAs. To do this, we first dissect the influence of DNAm and CNV on miRNA expression, respectively. We use DNAm and CNV as covariates in a regularized linear regression model and subsequently inferring the regulatory connectives. Specifically, we propose expression level of specific miRNA as Y ( Y miR ), and its methylation level as M ( M miR ), as well as the copy number level as N ( N miR ), we infer its association coefficients with the DNAm ( \({\beta }_{m}\) ) and CNV ( \({\beta }_{n}\) ) for each miRNA across all samples as: $${Y}_{miR}\approx {\beta }_{0}+{\beta }_{m}*{M}_{miR}+{\beta }_{n}*{N}_{miR}$$ We filtered those miRNAs with negative coefficients for DNAm and positive coefficients for CNV under FDR < 0.05. In this way, we identified a subset of 425 DNAm-miRs and 286 CNV-miRs that uniquely influenced by these two different mechanisms. Enrichment analysis indicated these DNAm-miRs are mainly involved in nervous system development and brain morphogenesis as well as some tumor related pathways, such as Wnt signaling and Notch signaling (Fig. 1 F). Whereas for CNV-miRs, functions are mainly about neuron development, cell proliferation and immune response pathways (Fig. 1 G). These observations indicated that CNV and DNAm mediated miRNA dysregulation may fulfil diverse functions in LGG, and thus it could help to define patients as more appropriate subtypes and to find suitable treatment methods. Molecular subtypes defined by CNV-miRs and DNAm-miRs Next, we investigated whether the expression profiles of CNV-miRs and DNAm-miRs can be used for patients’ stratification that associated with prognosis. We applied the non-negative matrix factorization (NMF) clustering analysis to miRNA expression profile with the cluster number k set as 2–10, the k values were then determined for both CNV-miRs and DNAm-miRs, respectively. The optimal number of k was selected as 4 based on profiles of DNAm-miRs (Fig. 2 A and Fig. S7), and as 3 based on profiles of CNV-miRs (Fig. 2 B and Fig. S8). For both clustering results that identified by DNAm-miRs and CNV-miRs, Kaplan–Meier plot analyses demonstrated that the molecular subtypes among each classification were significant associated with overall survival of the patients (Fig. 2 C-D). In addition, we compared the molecular subtypes obtained and the histological types and found that these classifications are distinct from each other (Fig. 2 E), whereas the comparison between the classifications generated from miRNA profiles indicated that high proportions of samples from CNV-miRs based subtypes were remarkably overlapped with the DNAm-miRs based ones (Fig. 2 F). These results indicate the heterogeneous for LGG samples at molecular level. We next sought to identify molecular subtypes that reflected the multi-layered pattern of the CNV-miRs and DNAm-miRs. We employed an integrated clustering method MOVICS for the combined analysis of the genomic profiling data regarding miRNA expression, DNAm and CNV. We identified four integrative clusters (ICs) from the 10 multi-omics ensemble clustering algorithms. The number of ICs was determined by referring to the cluster prediction index, gap statistical analysis and silhouette score, and was comprehensively evaluated by the consensus clustering matrix (Fig. S9A-C). Then, the clustering results were further combined by the consensus ensemble approach with distinctive molecular patterns of the multi-omics data (Fig. 3 A). Our classification system was closely related to overall survival of the patients (Log − rank p = 3.15e − 34; Fig. 3 B). Notably, cancer subtype 2 (IC2) exhibited the worst prognosis among all ICs evaluated (Fig. 3 C). For purpose of comparison, we also examined the overall survival among histological types, and only mild differences for their prognosis can be observed (Log − rank p = 0.005; Fig. 3 D). There are higher proportions of overlapping samples between histological type of astrocytoma and IC1 cluster, and between oligodendroglioma and IC3 cluster (Fig. 3 E). These results indicated that the integrated analysis of multi-omics data associated miRNAs identified clinically-relevant molecular subtypes in which miRNA dysregulation caused by epigenomic and genomic aberrance influences prognostic outcomes. To further identify the clinical characteristics of the four integrative clusters, we compared the tumor grade, 1p/19q co-deletion, MGMT promoter methylation and TERT expression status among the clusters identified. We found that patients in these groups exhibited asymmetric distributions for these characteristics. For example, the major proportion of IC2 samples are grade III tumor, which is significantly higher than that of other groups (p < 0.0001; Fig. S10A). In addition, all the IC1 and IC2 cluster samples present 1p/19q co-deletion, in contrast, nearly none of the samples in IC3 present co-deletion, whereas 64% of the samples in IC4 present co-deletion (p < 0.0001; Fig. S10B). There are also significant more patients in IC2 present MGMT promoter methylation (p < 0.0001; Fig. S10C), and fewer patients in IC2 and IC3 present TERT expression (p < 0.0001; Fig. S10D). We next evaluated the prognostic predictive value of the tumor subtypes in different grades. We performed survival analyses for the patients with grade II and III glioma respectively. Stratified survival analyses showed that patients in the IC2 group have a worse prognosis than patients in the other groups in both grade II and III samples (Fig. S10E-F), which indicated the prognostic independence of our glioma classifications method for patient survival. Distinct genomic features among different molecular subtypes Finally, we examined the single nucleotide variation (SNV) in different subtypes obtained to inspect the molecular heterogeneity of LGG samples. We screened a set of genes that altered significantly among subtypes, surprisingly, the somatic mutational spectrum of these genes agrees fantastically with the clusters identified. As shown in Fig. 4 A, we found that IDH1 mutates predominantly in clusters of IC1, IC3 and IC4, and the mutation frequency of this gene was further showed statistically significant differences among clusters (chi-square test p = 3.53e − 64, Fig. 4 B, left panel), whereas the TP53 predominates in IC1 (chi-square test p = 1.85e − 66, Fig. 4 B, middle panel). In contract, IC2 cluster bears nearly all the mutations for EGFR (chi-square test p = 3.33e − 29, Fig. 4 B, right panel). A mutually exclusive mutation pattern for the IDH1 and EGFR can be observed. In addition, some distinctly different mutant genes including ATRX , CIC and FUBP1 were also identified in our analysis. We next checked the associations between mutational spectrums of these genes and the patients’ survival. In general, only mild or no survival differences can be observed between groups defined by the mutational status of these genes (Fig. S11). This analysis indicates that our suggested classification system is superior to that based on mutations of individual genes for prognosis of LGG patients. We further extended our analysis to all mutated genes and calculated tumor mutation burden (TMB) for all the samples, as a result, IC2 presents highest mutation load across all clusters (Fig. 4 C), which indicated genetic heterogeneity are highly correlated with the subtypes defined by DNA methylation and copy number associated miRNAs. Inspired by the correlations between the mutational spectrum and the multi-omics clustering results for LGG, we further checked the pathways affected in order to obtain the potential targets for precise treatment, particularly for IC2. By using maftools, we identified the most significant pathways including “RTK-RAS”, “WNT”, “NOTCH” etc., which are differentially affected between IC2 and other groups for the genes involved (Fig. 4 D). The potential druggable gene categories for IC2 and other groups are shown in Fig. 4 E-F. We could see that the top list genes for IC2 includes the COL6A3 , EGFR , MUC16 , OBSCN and PIK3CA , whereas ARID1A , ATRX , BCOR , CIC and FUBP1 are included for other groups. These intriguing results strongly suggest that the distinct mutational events might synergistically involve in miRNA based LGG subtyping and thus provides novel insights to development of potential therapeutic strategies for clinical utility in the future. Immune landscape of LGG molecular subtypes The composition of the tumor immune microenvironment (TIME) has been identified as the critical factor of tumor–immune interactions to influence clinical outcome of patients and also the response to immune treatment ( 70 ). Many studies have been conducted on the TIME to uncover the underlying mechanisms of cancer ( 71 ). To evaluate whether the difference in multi-omics pattern imply immune heterogeneity and microenvironment of LGG, we thus investigated immune cell infiltrations and immune functions discrepancy among different clusters. We first focus on the 17 immunologically relevant pathways derived from ImmPort ( 44 ). We calculated the gene expression fold changes between IC2 and other clusters and performed the GSEA for the collected immune pathways. Analysis results showed that more than half of the pathways are enriched by the genes present dysregulation between glioma clusters, such as interleukins, chemokines, cytokines and cytokine receptors (Fig. 5 A). Next, we calculated the immune score by ESTIMATE for each tumor sample. We observe the immune scores were significantly increased in IC2 cluster, which suggests patients in IC2 cluster present higher infiltration of immune cells in tumor sample (Fig. 5 B). In addition, IC2 cluster also yields higher stromal scores than other groups, but lower tumor purity scores obtained by ESTIMATE (Wilcoxon p < 0.05, Fig. S12), which indicate that IC2 cluster samples are more heterogeneous than that of other groups. Moreover, we also calculated the antitumor immunoactivity for glioma patients, including CYT, MHC and CTL scores. All the three immunoactivity scores were significantly increased in IC2 cluster (Fig. 5 C, left panel), indicating patients in IC2 cluster displayed greater tumor cell cytotoxic activity and stronger antigen recognition capacities. In addition, we collected eight different immune signature gene sets and the ssGSEA was uses to estimate the enrichment score of tumor samples (Table S5), also, IC2 cluster also presents higher enrichment of different immune signatures (Fig. 5 C, right panel). For instance, higher IFN-γ score were observed in IC2 cluster, which indicated that IC2 may have stronger anti-tumor immunity and inhibited angiogenesis activities in tumor, thus to influence the efficiency of cancer immunotherapies ( 72 ). We next evaluated the relative proportion of immune cells of glioma patients using TIMER. Similarly, we found the infiltrate levels of CD4 + T, neutrophil, macrophage and dendritic cells were highest in IC2 among all cluster (Fig. 5 D), which presents the typical phenotypes for immune “hot” tumor. Finally, we collected 63 well-known marker genes of 21 immune cells from public resources (Table S6), the clustering heatmap shows that these genes exhibited diverse expression patterns among clusters, where samples in IC2 showed higher expression in general. (Fig. 5 E). Collectively, these results suggested the DNAm and CNV mediated the heterogeneity of miRNA expression is highly correlated with tumor immune status and could determine the immunophenotypes of glioma. Immune phenotype associated miRNA identification and prognosis model construction Having demonstrated the relationship between immune phenotype heterogeneity and the DNAm as well as CNV associated miRNA expression dysregulation, we next sought to comprehensively identify the miRNAs that determine the immune cell infiltrations and immune functions. By calculating the Pearson correlation coefficient between cell fractions and expression of DNAm-miRs and CNV-miRs, we identified 517 and 325 miRNAs significantly correlated with infiltrating immune cells respectively. In general, the number of positively correlated miRNAs was higher than that of negatively correlated ones for both groups (Fig. 6 A). Interestingly, we observed that the number of miRNAs with immune cell infiltration associations increases with that of cell types, and there are 195 miRNAs were identified that correlated with all the six immune cell infiltrations (Fig. 6 B, Table S7). In addition, we observe that miRNA expression correlated with diverse immune cell type infiltration in majority of the clusters identified. Most miRNAs were positively associated with immune cells, whereas the negative correlations are mainly occurring in CD4 + T cells (Fig. 6 C). In order to further explore effects of the immune-related key miRNAs on LGG, we collected experimentally verified immune related miRNAs from the RNA2Immune database, which is a public repository to catalogue ncRNA-immune associations ( 56 ). Here we selected 464 and 574 miRNAs from the sections of “immune cell function” and “cancer immunity”, respectively (Table S8). By comparing with those miRNAs that correlated with immune cell infiltrations, a total of 64 miRNAs were identified that overlap with the external validations (Fig. 6 D), in which 47 miRNAs are from DNAm-miRs and 17 are from CNV-miRs (Fig. 6 E). Studies have extensively indicated that miRNAs present strong potential for predicting cancer prognosis and exhibited great potential for clinical applications as therapeutic targets ( 73 ). Thus, we integrated the survival information to evaluate the potential ability of the identified miRNA as candidate prognostic marker for LGG. We performed a Cox proportional hazard regression analyses to identify the survival related miRNAs. Under a stringent threshold of p < 0.0005, a total of 34 miRNAs were identified that correlated with patients’ survival. Interestingly, we found that almost all of these miRNAs are positively correlated with immune cell infiltration levels, and only miR-9-5p exhibited negative correlations (Fig. 6 F). To construct the prognostic prediction model, a LASSO proportional hazards regression analysis was then performed to further determine the relationship between patients’ survival and the expressions of the identified 34 miRNAs. We found a significant correlation with the overall survival of LGG patients with the following 8 miRNAs: (hsa-miR-10b-5p, hsa-miR-155-5p, hsa-miR-196a-5p, hsa-miR-196b-5p, hsa-miR-200a-3p, hsa-miR-204-5p, hsa-miR-503-5p, hsa-miR-15b-5p). These markers were considered to function as prognostically significant immune-related miRNAs. We then implemented these 8 miRNAs in development of a miRNA based prognostic risk score that determine the degree of tumor immune infiltration. These 8 miRNAs were weighted based on the LASSO regression coefficient as follows: Risk Score = (0.1617* Exp miR−10b−5p ) + (0.3419* Exp miR−155−5p ) + (0.0695* Exp miR−196a−5p ) + (0.1391* Exp miR−196b−5p ) + (0.1474* Exp miR−200a−3p ) + (0.0840* Exp miR−204v5p ) + (0.0996* Exp miR−503−5p ) + (0.2083* Exp miR−15b−5p ) We observed that the IC2 cluster samples present height risk scores among all groups as expected (Wilcoxon p < 2.2e − 16, Fig. S13). Besides, all the samples in this training set were separated into high-risk (N = 84) and low-risk groups (N = 411) according to the level of the risk scores, the distribution of risk score among high-risk and low-risk groups and expression heatmap of the miRNAs are indicated in Fig. 6 G. We can observe a clear distinction between for the signature miRNA expression from the two groups of patients. In order to determine whether the prognostic model was robust, a validation set obtained from CGGA was utilized to verify the validity and accuracy of the prognostic model based on the 8 miRNAs. We collected miRNA microarray data for 198 samples from this cohort, 188 of which have the survival data available, each sample of the validation set also acquired a risk score according to the same formula, and was then stratified as being either high-risk (N = 121) or low-risk groups (N = 67). The risk score distribution and gene expression data for the validation set are shown in Fig. S14. Patients with higher risk scores were noted to have worse survival in contrast to those with lower risk scores, the overall survival analysis revealed significant difference between the two groups (Log − rank p = 4.24e − 08, Fig. 6 H). The predictive potential of the prognostic model using time-dependent ROC curves was shows as Fig. 6 I. The area under the ROC curve (AUC) of the prognostic model for OS was 0.65 at 1 year, 0.75 at 2 years, 0.78 at 3 years and 0.76 at 5 years, suggesting that our prediction model has favorable efficacy for predicting both short- and long-term prognosis. qRT-PCR validation of signature miRNAs and clinical utility of risk model We next validated the miRNAs as novel prognosis markers for glioma in our in-house cohort. We detected the expression level of the eight miRNAs form tumor and normal samples of 30 LGG patients by qRT-PCR (Table S9). The normalized expression levels from qRT-PCR showed that all these miRNAs present elevated expression in LGG tumor sample, although with different significance obtained. Two miRNAs of miR-10b-5p and miR-204-5p present mild upregulation for about 1.5 times of fold change, whereas others are significantly up-regulated in about 20–50 times, Notably, miR-155-5p showed highest level of fold change in 56.4 times (Fig. 7 A). These cross platform and cross-racial data indicated the robustness of these signature miRNAs presented here for patient’s prognosis. Furthermore, we conducted univariate and multivariate Cox regression analyses to explore whether the prognostic value of our risk model was independent of other clinical features in the TCGA LGG cohort. For the univariate Cox regression model, our prediction classifier was the strongest variable correlated with prognosis among clinical features including histological type, IDH mutation, 1p/19q codeletion, MGMT promoter methylation and TERT expression status. After multivariate adjustment by the clinical features, our classifier remains the strongest and independent prognostic factor in survival analysis (Fig. 7 B). This results indicated that our model performed better than conventional clinicopathological features for clinical outcome prediction. To provide clinically correlated quantitative approach to predicting the prognosis of LGL patients, a nomogram that integrated the risk model and independent clinical risk factors (tumor grade and 1p/19q codeletion) was constructed (Fig. 7 C). Our risk model was found to contribute the most risk points (ranging from 0 to 100) compared with the other clinical features, which was consistent with the Cox multivariate regression results. The calibration plot indicated that the nomogram performed well compared with an ideal model (Fig. 7 D). The AUC score at 1-, 2-, 3- and 5- years were 0.86, 0.84, 0.83 and 0.80 for the nomogram, respectively (Fig. 7 E). In summary, these findings suggest that the nomogram was a better model for predicting short-term or long-term survival of LGG patients than individual prognostic factors. miRNAs associated with immune evasion of LGG Since the IC2 cluster exhibited higher levels of immune cell infiltration, whereas patients in this cluster had poorer prognosis than others, we conducted in-depth functional investigations to explore possible reasons for this contradiction. We first performed functional enrichment analysis of GSVA on Hallmark gene sets obtained from MSigDB. Specially, besides the higher immune activities, the IC2 cluster is also endowed with higher energy expenditure states, such as glycolysis, hypoxia and angiogenesis, whereas IC4 is more closely related to metabolic activities, such as the bile acid and fatty acid metabolism, as well as protein secretion. In contract, IC3 and IC1 clusters are more distinguished by proliferative-related pathways, such as mitotic spindle, E2F targets, and G2M checkpoints, and also by classical tumor pathways including the MYC and Wnt-signaling (Fig. 8 A). Next, we further investigated whether the clusters obtained are correlated with immunosuppressive effects. We used the methodology from Shuichi et al . to calculate three different scores in order to measure tumor proliferation and glycolysis rates (T-score), immune cell and antitumor immunity level (I-score), as well as immunosuppressive level (S-score), respectively. Here, the IC2 cluster conveys the highest for all three different scores, particularly for the S-score, which is much higher in IC2 than other groups (Fig. 8 B), which indicated an immune evasion state may occurred in the patients of IC2 cluster. We next inspected PD-L1 expression and found its expression is much higher in IC2 than other groups (Fig. 8 C), which verify the hypothesis of an immune evasion status for the IC2 cluster. To identify key miRNAs that involved in LGG immune evasion, we collected experimentally validated regulatory associations between PD-L1 and miRNAs from miRTarBase, and further screened those by performing a Pearson correlation analysis between the expression of PD-L1 and DNAm-miRs as well as CNV-miRs. Base on negative correlation of PCC < 0 and the threshold of FDR < 0.05, we discovered many miRNA modulators that play key roles in PD-L1 mediated immune evasion. One of the examples is the miR-338-5p, our analysis indicated this miRNA presents significant negative correlation with PD-L1 for their expressions (r = − 0.1, p = 0.0309) (Fig. 8 D), this is possibly due to the DNA methylation dysregulation in its promoter region, a negative correlation between expression of this miRNA and the promoter methylation was observed (r = − 0.15, p = 0.0011) (Fig. 8 E). This regulatory relationship between miR-338-5p and PD-L1 has been identified for regulatory T cells (Tregs) mediated immune evasion ( 74 ). Another example is the miR-17-5p, whose expression is significantly negatively correlated with that of PD-L1 (r = − 0.12, p = 0.0056) (Fig. S15A), and its expression dysregulation was identified that caused by its copy number variation (r = 0.24, p = 7.14e − 08) (Fig. S15B). Inverse correlation between PD-L1 and miR-17-5p has been confirmed in many cancer types. For instance, up-regulation of PD-L1 due to post-transcriptional control by miR-17-5p confers resistance to BRAF or MEK inhibitors (BRAFi or MEKi) for metastatic melanoma ( 75 ). In addition, ceRNA regulatory loops involving miR-17-5p/PD-L1 have been observed in lung cancer and breast cancer to promote tumor progression and immune checkpoint blockades (ICB) resistance ( 76 , 77 ). Besides these, similar regulatory interactions were also identified for eight miRNAs including let-7a-3p, miR-30b-3p, miR-9-3p, miR-548f-3p, miR-139-3p miR-9-5p, miR-128-3p and miR-4797-3p. These miRNAs present significant negative correlations for their expression with PD-L1 (Fig. S16), which indicated this miRNA regulatory module may play critical roles in PD-L1 associated immune evasion for LGG (Fig. 8 F). To provide a more comprehensive view of miRNA related tumor immune evasion mechanism, we used TIDE method to evaluate the efficacy of ICB based immunotherapy for LGG. We found patients in the IC2 and IC3 clusters presents significant higher TIDE scores than other groups, which results in a lower efficacy of immunotherapy (Fig. 9 A). According to TIDE analysis, only about one quarter of the patients in IC2 and IC3 may response to PD-1 and CTLA-4 based treatment, in contract, about 64%-79% of the patients in IC1 and IC4 could benefit from the treatment (Fig. 9 B). By Pearson correlation analysis between TIDE scores and miRNA expression, we identified 225 miRNAs with significant correlation under FDR < 0.05, which may involve in immunotherapy response. One such example is miR-10b-5p, a significant positive correlation was observed between TIDE score and miR-10b-5p expression (Fig. 9 C). This miRNA has been identified that impairs TET2-mediated PD-L1 transcription inhibition, by which to promote immune evasion and tumor progression in glioblastoma ( 78 ). By further integration of the TIDE score correlated miRNAs and the putative targets obtained from public resources, we constructed a miRNA immune evasive network to describe the mechanism explaining roles of miRNA involved in immune envision regulation (Fig. 9 D). Within this network, 29 unique targets of 82 immune envision related miRNAs that under the CNV and DNAm burden in LGG were identified. We screened genes that simultaneously targeted by two or more miRNAs and found PLAGL2 , SCD1 , PIK3CB , ELOVL6 , NR2E1 and MET demonstrated the highest connectivity, which suggests the importance of these genes in immune envision and the associated immunotherapy response for LGG. Indeed, many of these genes have been validated for their functions in immune regulation and cancer development. For instance, the PLAGL2 and the PIK3CB were identified as key hub gene for the regulatory network that significantly associated with tumor immune infiltration and sensitivity of chemotherapeutic drugs for cancer ( 79 , 80 ). SCD1 was found that expressed in cancer cells and immune cells causes immune-resistant conditions, inhibition of SCD1 enhances the recruitment of dendritic cells into tumors and the subsequent induction and tumor accumulation of antitumor CD8 + T cells, thus this gene was identified as a potential target to enhance the antitumor effects of an anti-PD-1 therapy ( 81 ). For miRNAs, miR-214-3p, miR-25-3p and miR-27b-3p were identified as key hub regulators with highest connectivity, and also have been validated for their roles in immune escape ( 82 , 83 ). In summary, our analysis revealed miRNAs act as key regulators for genes that involved in cancer immune evasion, and thus could act as an adjuvant to several major immunotherapies by synergistically increasing their efficacy potential. Validation of the association between prognosis related miRNAs and immune evasion in vitro experiments We next evaluated the capacities of prognosis related miRNAs as novel immune evasion regulators for glioma in vitro experiments. miRNA inhibitor technology was used to knock down the miRNAs in the U251 cell line and measure the expression of miRNAs. We found that miRNA inhibitor significantly reduced the expression of miRNAs compared with transfected NC inhibitor cells (Fig. 10 A). We further performed CCK8 to investigate the effect of miRNAs on the viability of glioma cells. Silencing of miR-155-5p, miR-196a-5p, miR-196b-5p, miR-200a-3p, miR-503-5p, and miR-15b-5p delayed the proliferation of U251 cell lines, while silencing of miR-10b-5p or miR-204-5p has no effect (Fig. 10 B). Wound healing and transwell assays were employed to explore whether cell migration and invasion were influenced after silencing of miRNAs. We found that the knockdown of miR-155-5p, miR-196a-5p, miR-196b-5p, miR-200a-3p, miR-503-5p, and miR-15b-5p significantly decreased cell migration and invasion compared with the control, while silencing of miR-10b-5p or miR-204-5p has no effect (Fig. 10 C-F). Moreover, the knockdown of miR-155-5p, miR-196a-5p miR-196b-5p, miR-200a-3p, miR-503-5p, and miR-15b-5p significantly inhibited the well-known immune evasion markers of PD-L1 , CTLA4 , and FOXP3 at the protein level, while silencing of miR-10b-5p or miR-204-5p has no effect (Fig. 10 G). Taken together, these results suggested that miR-155-5p, miR-196a-5p, miR-196b-5p, miR-200a-3p, miR-503-5p, and miR-15b-5p significantly promote in vitro cell migration, invasion, proliferation, and are associated with immune evasion in glioma. Discussion MiRNAs expression deregulation has been widely observed in cancer development and progression. However, less attention has been paid towards understanding the underlying mechanisms of miRNA expression dysregulation. In recent years, more and more evidences have indicated the genetic and epigenetic machineries are extensively involved in the process of miRNA biogenesis and impair its expression in cancer cells ( 84 , 85 ). Here we focus on two widely-known genetic and epigenetic processes of DNA methylation and copy number variation, in order to get deeper insights into this intricate regulation as well as their roles in cancer development. We developed strategies to obtain promoter methylation and copy number variation for miRNA genes from high-throughput profiling, combined with miRNA expression, we obtained miRNAs whose expression affected by DNAm and CNV dysregulation in lower grade glioma, and also identified prognostic subtypes based on multidimensional omics analysis. Different with previous studies which only investigate individual mechanism ( 86 , 87 ), we sought to inspect and give a more comprehensive description for both effect of DNAm and CNV on miRNA expression simultaneously. Our analysis not only provides novel molecular subtyping system for LGG, but also reveals that the regulatory mechanism surrounding miRNA may be closely related to tumor microenvironment and immune evasion, as well as the resulting prognostic outcome differences. Furthermore, this tumor classification around miRNAs is also closely related to genome wide single nucleotide polymorphisms. This work indicated that integrative analysis of multi-omics profiling focus on miRNA regulation could identify meaningful molecular subtypes which are in connection with mechanisms and sample heterogeneity in tumors as well as potential biomarkers and therapy targets. The relationship between genetic and epigenetic regulatory mechanisms has been extensively investigated by using integrative multi-omics study for molecular subtyping ( 88 , 89 ). Nonetheless, illustrating the complex interactions for DNAm- or CNV- related miRNA dysregulation remains limited, mainly due to that obtaining high-throughput profiles around miRNA at genetic and epigenetic layers is still challenging. Toward this end, we developed novel strategy and successfully extracted miRNA promoter methylation information by using credible miRNA promoter annotation and repurposing the 450K microarray data, which were originally designed to detect the methylation of protein coding genes. Although many miRNA genes reside within intronic regions and are thought to co-transcribed with their host genes, increasing evidences have indicated that intronic miRNAs carry their own promoters, therefore have independent transcription ( 90 ). In this case, we decided to use experimentally verified promoter region for individual miRNAs but without that from host genes being used. In addition, we also obtained miRNA CNV profiling by examining miRNAs loci which localize in somatic copy number variation regions. Our results confirmed such associations between DNAm or CNV aberration and miRNA expression. Akin to protein coding genes, the methylation at promoter regions and the expression for corresponding miRNAs are often observed to have negative correlation, whereas for CNVs, the positive correlations dominant for all associations, although few opposite trends also can be observed. We found that about one third of the miRNAs in the human genome tend to be affected by either copy number or methylation alterations, and many of them were driven by a combination of both mechanisms. This phenomenon was also observed for protein coding genes from glioblastoma and ovarian cancer ( 91 ), suggesting the crosstalk between genetic and epigenetic layers are widespread in cancer cell, this synergistic regulation between different regulatory mechanisms of DNAm and CNV may function in maintaining the precise expression of genes. The concomitant regulation for miRNA expression by alterations in genomic CNV and DNAm was also indicted by the genomic distribution of miRNA genes affected. Surprisingly, both the CNV-miRs and DNA-miRs exhibited genomic preference for three different chromosomes, especially the chromosome 14. In addition, the frequencies of the aberrant DNA copy number and DNA methylation are also correlated. The concentration distribution of miRNA genes may raise from expansion of miRNA gene families in the genome. Several studies have indicated that some miRNA families may evolved from whole genome or tandem duplication events, their family members underwent fast amplification and rapid divergence, thus may have conserved regulatory elements ( 92 ). This could be another mechanism for elaborating the regulation of miRNA family members. This co-regulation in the genome may allow them to target the genes within key pathways that involved in cancer, and also provides us another perspective to understand the mechanism of miRNA evolution. Our classification analysis based on CNV and DNAm associated miRNAs revealed novel molecular key features that represent novel precision therapy targets and biomarkers for LGG management. These LGG samples were subdivided into four subtypes which shows significant prognostic difference. Survival analysis showed that patients’ survival was much lower in the IC2 subtype than that in other three subtypes. Results indicated that our classification is superior to those of protein-coding gene based mutational profiling or multi-omics classification methods for prediction of patients’ prognosis ( 93 , 94 ). As a kind of highly effective gene expression regulator, one miRNA may regulate dozens to hundreds of target genes. One could hypothesize that disturbance to only a few miRNAs may cause huge downstream runaway effect ( 95 ). In such case, miRNAs show a higher value for tumor classification, development of relevant tools around miRNA omics data may provide us a better alternative for early diagnosis and accurate prognosis predictions in cancers. Analysis of different biological characteristics among subtypes revealed intra-tumor heterogeneity and molecular interpretability of LGG. Significant differences in the tumor immune microenvironment of the molecular subtypes were observed. Studies have extensively shown that tumor infiltrating lymphocytes (TILs) are involved in tumor progression and invasiveness. TILs include various lymphocytes with different activities, particularly the CD8 + and CD4 + T cells, which are usually associate with favorable prognostic in patients ( 96 ). As found in our study, higher level of TILs but the poor survival of IC2 indicated this group as an immunosuppressive subtype. miRNAs are emerging as important regulators involved in the cancer immune process and play essential roles in cancer immunotherapy ( 97 ). Therefore, it is necessary to determine the relationship between miRNAs and immune regulation in LGG. By estimating the correlations with infiltrating immunocytes and immune pathway activities, we identified immune-related miRNAs that may have prognosis value for LGG patients with 8 essential immune miRNAs were prioritized, many of them have been reported to be associated with immune regulation and cancer therapy from previous studies. For instance, miR-15b-5p mediates the degradation of PD-L1 mRNA through interaction with its 3'-untranslated region, by which to induce CD8 + T and NK cell activation and cytotoxicity against neuroblastoma ( 98 ). Similar observations also apply to another miRNA of miR-155-5p, by targeting PD-L1 , this miRNA was found to activated CD8 + T cell function and suppress ovarian cancer progression ( 99 ). The prognosis model constructed based on the expression of these miRNAs was validated by an independent cohort. It is worth noting that despite the different platforms and technologies being used to assess expression status, these miRNAs were consistently found that associated with patient’s risk in the discovery and validation cohort, suggesting the robust of the proposed prediction model. These results indicated that the integration of DNAm and CNV is an effective method to study the immune regulation of miRNAs for cancer prognosis. Remarkably, our proposed sample classification has a good concordance with the mutation profiles among the subtypes, particularly for IDH1 , TP53 and EGFR . Somatic mutation is another mechanism that may impact the biogenesis of miRNAs. However, the detailed process remains unrevealed. One inspiring example is the IDH1 , somatic mutations of this gene are strongly associated with the DNA methylome and transcriptome reorganization in LGG ( 100 ). Whether DNA methylation mediates this indirect association between IDH1 mutation and miRNA expression remains further analysis. Besides, the association of miRNA expression and mutation frequencies of other genes within LGG remains unknown. Nonetheless, our present findings provide potential therapy targets for corresponding subtypes and brings novel enlightenment to the study of the intrinsic regulatory mechanisms for LGG. Conclusions In summary, we have designed and implemented a novel strategy to facilitate the integration of copy number variation, DNA methylation and miRNA expression data. Analysis results demonstrate the added value of for molecular subtyping of LGG by applying our method to the multi-omics data. Systematic study also characterized the immune landscape for glioma patients based on miRNAs related multi-omics data and depicted the crosstalk among DNA methylation, copy number variation and miRNA expression for immune regulation and evasion. The repertoire of immune-related miRNAs will facilitate the development of immunotherapeutic targets for glioma in future. Declarations Authors’ contributions Z.Y, K.L and H.Y conceived the study. X.L and H.X performed experimental validation analysis. A.E.T supervised data analysis. L.X and J.L performed data analysis. M.F, J.L, H.Z, Y.W, L.Z, G.Z and Y.H collected and processed samples. Z.Y wrote the paper. A.E.T revised the paper. All authors read and approved the paper. Funding The study was supported by the National Natural Science Foundation of China (91959106, 32071270, 82072370, 82103376); Outstanding Innovative Research Team for Molecular Enzymology and Detection in Anhui Provincial Universities (2022AH010012); Program for Excellent Sci-tech Innovation Teams of Universities in Anhui Province (2022AH010074); the Major Science and Technology Projects in Anhui Province (202003a06020009); Anhui Natural Science Foundation (2208085MC48); Key University Science Research Project of Anhui Province (2023AH051746, 2023AH051768); Climbing Peak Training Program for Innovative Technology team of Yijishan Hospital, Wannan Medical College (PF201904); Peak Training Program for Scientific Research of Yijishan Hospital, Wannan Medical College (GF2019T01, GF2019G15); the Open Project of Key Laboratory of Non-coding RNA Transformation Research of Anhui Higher Education Institution, Wannan Medical College (RNA202205); Science and Technology Application Basic Research Project of Wuhu (2022jc60); and Key Research and Development Foundation supported by Science and Technology Department of Sichuan Province (2023YFS0243). The sponsors have no role in study design, in the collection, analysis and interpretation of data, in the writing of the report, and in the decision to submit the article for publication. Availability of data and materials Data analyzed in this manuscript is already publicly available from The Cancer Genome Atlas (TCGA) data portal: https://portal.gdc.cancer.gov/ and the Chinese Glioma Genome Atlas (CGGA) data portal: http://www.cgga.org.cn/. Ethics approval and consent to participate Experimental methods using human glioma tissues were approved by the Medical Ethics Committee of Wannan Medical College and are in compliance with the Helsinki Declaration of the world medical association. All of the donors gave written informed consent for their sample donation could also be used for scientific investigations. Consent for publication Not applicable. Competing interests The authors declare that they have no competing interests. References Ostrom QT, Bauchet L, Davis FG, Deltour I, Fisher JL, Langer CE, et al. The epidemiology of glioma in adults: a state of the science review. Neuro Oncol. 2014;16(7):896–913. Goodenberger ML, Jenkins RB. Genetics of adult glioma. Cancer Genet. 2012;205(12):613–21. Weller M, Stupp R, Reifenberger G, Brandes AA, van den Bent MJ, Wick W, et al. MGMT promoter methylation in malignant gliomas: ready for personalized medicine? Nat Rev Neurol. 2010;6(1):39–51. van den Bent MJ. Chemotherapy for low-grade glioma: when, for whom, which regimen? Curr Opin Neurol. 2015;28(6):633–938. Yang K, Wu Z, Zhang H, Zhang N, Wu W, Wang Z, et al. Glioma targeted therapy: insight into future of molecular approaches. Mol Cancer. 2022;21(1):39. Louis DN, Perry A, Reifenberger G, von Deimling A, Figarella-Branger D, Cavenee WK, et al. The 2016 World Health Organization Classification of Tumors of the Central Nervous System: a summary. Acta Neuropathol. 2016;131(6):803–20. Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004;116(2):281–97. O'Brien J, Hayder H, Zayed Y, Peng C. Overview of MicroRNA Biogenesis, Mechanisms of Actions, and Circulation. Front Endocrinol (Lausanne). 2018;9:402. Shivdasani RA. MicroRNAs: regulators of gene expression and cell differentiation. Blood. 2006;108(12):3646–53. Hwang HW, Mendell JT. MicroRNAs in cell proliferation, cell death, and tumorigenesis. Br J Cancer. 2007;96 Suppl:R40-4. Tufekci KU, Oner MG, Meuwissen RL, Genc S. The role of microRNAs in human diseases. Methods Mol Biol. 2014;1107:33–50. Paul P, Chakraborty A, Sarkar D, Langthasa M, Rahman M, Bari M, et al. Interplay between miRNAs and human diseases. J Cell Physiol. 2018;233(3):2007–18. Ohno M, Matsuzaki J, Kawauchi J, Aoki Y, Miura J, Takizawa S, et al. Assessment of the Diagnostic Utility of Serum MicroRNA Classification in Patients With Diffuse Glioma. JAMA Netw Open. 2019;2(12):e1916953. Zhang Y, Dutta A, Abounader R. The role of microRNAs in glioma initiation and progression. Front Biosci (Landmark Ed). 2012;17(2):700–12. Li Y, Xu J, Chen H, Bai J, Li S, Zhao Z, et al. Comprehensive analysis of the functional microRNA-mRNA regulatory network identifies miRNA signatures associated with glioma malignant progression. Nucleic Acids Res. 2013;41(22):e203. Wang Q, Li P, Li A, Jiang W, Wang H, Wang J, et al. Plasma specific miRNAs as predictive biomarkers for diagnosis and prognosis of glioma. J Exp Clin Cancer Res. 2012;31(1):97. Zhang Y, Chen J, Xue Q, Wang J, Zhao L, Han K, et al. Prognostic Significance of MicroRNAs in Glioma: A Systematic Review and Meta-Analysis. Biomed Res Int. 2019;2019:4015969. Tumilson CA, Lea RW, Alder JE, Shaw L. Circulating microRNA biomarkers for glioma and predicting response to therapy. Mol Neurobiol. 2014;50(2):545–58. Mondal I, Kulshreshtha R. Potential of microRNA based diagnostics and therapeutics in glioma: a patent review. Expert Opin Ther Pat. 2021;31(1):91–106. Ahmadpour S, Taghavi T, Sheida A, Tamehri Zadeh SS, Hamblin MR, Mirzaei H. Effects of microRNAs and long non-coding RNAs on chemotherapy response in glioma. Epigenomics. 2022;14(9):549–63. Xiao Y, Bi M, Guo H, Li M. Multi-omics approaches for biomarker discovery in early ovarian cancer diagnosis. EBioMedicine. 2022;79:104001. Heo YJ, Hwa C, Lee GH, Park JM, An JY. Integrative Multi-Omics Approaches in Cancer Research: From Biological Networks to Clinical Subtypes. Mol Cells. 2021;44(7):433–43. Pos O, Radvanszky J, Buglyo G, Pos Z, Rusnakova D, Nagy B, et al. DNA copy number variation: Main characteristics, evolutionary significance, and pathological aspects. Biomed J. 2021;44(5):548–59. Krijgsman O, Carvalho B, Meijer GA, Steenbergen RD, Ylstra B. Focal chromosomal copy number aberrations in cancer-Needles in a genome haystack. Biochim Biophys Acta. 2014;1843(11):2698–704. De Carvalho DD, You JS, Jones PA. DNA methylation and cellular reprogramming. Trends Cell Biol. 2010;20(10):609–17. Blattler A, Farnham PJ. Cross-talk between site-specific transcription factors and DNA methylation states. J Biol Chem. 2013;288(48):34287–94. Borgel J, Guibert S, Li Y, Chiba H, Schubeler D, Sasaki H, et al. Targets and dynamics of promoter DNA methylation during early mouse development. Nat Genet. 2010;42(12):1093–100. Fialkova V, Vidomanova E, Balharek T, Marcinek J, Kudela E, Hanysova S, et al. DNA methylation as mechanism of apoptotic resistance development in endometrial cancer patients. Gen Physiol Biophys. 2017;36(5):521–9. Sun W, Bunn P, Jin C, Little P, Zhabotynsky V, Perou CM, et al. The association between copy number aberration, DNA methylation and gene expression in tumor samples. Nucleic Acids Res. 2018;46(6):3009–18. Ali Syeda Z, Langden SSS, Munkhzul C, Lee M, Song SJ. Regulatory Mechanism of MicroRNA Expression in Cancer. Int J Mol Sci. 2020;21(5). Cancer Genome Atlas, Research N, Brat DJ, Verhaak RG, Aldape KD, Yung WK, Salama SR, et al. Comprehensive, Integrative Genomic Analysis of Diffuse Lower-Grade Gliomas. N Engl J Med. 2015;372(26):2481–98. Di Lena P, Sala C, Prodi A, Nardini C. Missing value estimation methods for DNA methylation data. Bioinformatics. 2019;35(19):3786–93. Teschendorff AE, Marabita F, Lechner M, Bartlett T, Tegner J, Gomez-Cabrero D, et al. A beta-mixture quantile normalization method for correcting probe design bias in Illumina Infinium 450 k DNA methylation data. Bioinformatics. 2013;29(2):189–96. Zhao Z, Zhang KN, Wang Q, Li G, Zeng F, Zhang Y, et al. Chinese Glioma Genome Atlas (CGGA): A Comprehensive Resource with Functional Genomic Data from Chinese Glioma Patients. Genomics Proteom Bioinf. 2021;19(1):1–12. Xu F, Wang Y, Ling Y, Zhou C, Wang H, Teschendorff AE, et al. dbDEMC 3.0: Functional Exploration of Differentially Expressed miRNAs in Cancers of Human and Model Organisms. Genomics Proteom Bioinf. 2022;20(3):446–54. Kozomara A, Birgaoanu M, Griffiths-Jones S. miRBase: from microRNA sequences to function. Nucleic Acids Res. 2019;47(D1):D155–D62. Mermel CH, Schumacher SE, Hill B, Meyerson ML, Beroukhim R, Getz G. GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers. Genome Biol. 2011;12(4):R41. de Rie D, Abugessaisa I, Alam T, Arner E, Arner P, Ashoor H, et al. An integrated expression atlas of miRNAs and their promoters in human and mouse. Nat Biotechnol. 2017;35(9):872–8. Mirzal A. Nonparametric Tikhonov Regularized NMF and Its Application in Cancer Clustering. IEEE/ACM Trans Comput Biol Bioinform. 2014;11(6):1208–17. Lu X, Meng J, Zhou Y, Jiang L, Yan F. MOVICS: an R package for multi-omics integration and visualization in cancer subtyping. Bioinformatics. 2021;36(22–23):5539–41. Pierre-Jean M, Deleuze JF, Le Floch E, Mauger F. Clustering and variable selection evaluation of 13 unsupervised methods for multi-omics data integration. Brief Bioinform. 2020;21(6):2011–30. Li T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, et al. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res. 2020;48(W1):W509–W14. Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612. Bhattacharya S, Dunn P, Thomas CG, Smith B, Schaefer H, Chen J, et al. ImmPort, toward repurposing of open access immunological assay data for translational and clinical research. Sci Data. 2018;5:180015. Rooney MS, Shukla SA, Wu CJ, Getz G, Hacohen N. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell. 2015;160(1–2):48–61. Narayanan S, Kawaguchi T, Yan L, Peng X, Qi Q, Takabe K. Cytolytic Activity Score to Assess Anticancer Immunity in Colorectal Cancer. Ann Surg Oncol. 2018;25(8):2323–31. Lauss M, Donia M, Harbst K, Andersen R, Mitra S, Rosengren F, et al. Mutational and putative neoantigen load predict clinical benefit of adoptive T cell therapy in melanoma. Nat Commun. 2017;8(1):1738. Liu Y, Liang G, Xu H, Dong W, Dong Z, Qiu Z, et al. Tumors exploit FTO-mediated regulation of glycolytic metabolism to evade immune surveillance. Cell Metab. 2021;33(6):1221–33e11. Ju M, Bi J, Wei Q, Jiang L, Guan Q, Zhang M et al. Pan-cancer analysis of NLRP3 inflammasome with potential implications in prognosis and immunotherapy in human cancer. Brief Bioinform. 2021;22(4). Thorsson V, Gibbs DL, Brown SD, Wolf D, Bortone DS, Ou Yang TH, et al. The Immune Landscape of Cancer. Immunity. 2018;48(4):812–30. e14. Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7. Ru B, Wong CN, Tong Y, Zhong JY, Zhong SSW, Wu WC, et al. TISIDB: an integrated repository portal for tumor-immune system interactions. Bioinformatics. 2019;35(20):4200–2. Shinohara S, Takahashi Y, Komuro H, Matsui T, Sugita Y, Demachi-Okamura A et al. New evaluation of the tumor immune microenvironment of non-small cell lung cancer and its association with prognosis. J Immunother Cancer. 2022;10(4). Liberzon A, Birger C, Thorvaldsdottir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 2015;1(6):417–25. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med. 2018;24(10):1550–8. Wang J, Li S, Wang T, Xu S, Wang X, Kong X et al. RNA2Immune: A Database of Experimentally Supported Data Linking Non-coding RNA Regulation to the Immune System. Genomics Proteom Bioinf. 2022. Smyth GK. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. 2004;3:Article3. Nam JW, Rissland OS, Koppstein D, Abreu-Goodger C, Jan CH, Agarwal V, et al. Global analyses of the effect of different cellular contexts on microRNA targeting. Mol Cell. 2014;53(6):1031–43. Chen Y, Wang X. miRDB: an online database for prediction of functional microRNA targets. Nucleic Acids Res. 2020;48(D1):D127–D31. Huang HY, Lin YC, Cui S, Huang Y, Tang Y, Xu J, et al. miRTarBase update 2022: an informative resource for experimentally validated miRNA-target interactions. Nucleic Acids Res. 2022;50(D1):D222–D30. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innov (Camb). 2021;2(3):100141. Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018;28(11):1747–56. Kim H, Huang W, Jiang X, Pennicooke B, Park PJ, Johnson MD. Integrative genome analysis reveals an oncomir/oncogene cluster regulating glioblastoma survivorship. Proc Natl Acad Sci U S A. 2010;107(5):2183–8. Favero F, McGranahan N, Salm M, Birkbak NJ, Sanborn JZ, Benz SC, et al. Glioblastoma adaptation traced through decline of an IDH1 clonal driver and macro-evolution of a double-minute chromosome. Ann Oncol. 2015;26(5):880–7. Schliesser MG, Claus R, Hielscher T, Grimm C, Weichenhan D, Blaes J, et al. Prognostic relevance of miRNA-155 methylation in anaplastic glioma. Oncotarget. 2016;7(50):82028–45. Wu X, Wan Q, Wang J, Hou P, Zhang Q, Wang Q, et al. Epigenetic Activation of lncRNA MIR155HG Mediated by Promoter Hypomethylation and SP1 is Correlated with Immune Infiltration in Glioma. Onco Targets Ther. 2022;15:219–35. Vrba L, Munoz-Rodriguez JL, Stampfer MR, Futscher BW. miRNA gene promoters are frequent targets of aberrant DNA methylation in human breast cancer. PLoS ONE. 2013;8(1):e54398. Li CL, Nie H, Wang M, Su LP, Li JF, Yu YY, et al. microRNA-155 is downregulated in gastric cancer cells and involved in cell metastasis. Oncol Rep. 2012;27(6):1960–6. Sandoval J, Diaz-Lagares A, Salgado R, Servitje O, Climent F, Ortiz-Romero PL, et al. MicroRNA expression profiling and DNA methylation signature for deregulated microRNA in cutaneous T-cell lymphoma. J Invest Dermatol. 2015;135(4):1128–37. Binnewies M, Roberts EW, Kersten K, Chan V, Fearon DF, Merad M, et al. Understanding the tumor immune microenvironment (TIME) for effective therapy. Nat Med. 2018;24(5):541–50. Tang T, Huang X, Zhang G, Hong Z, Bai X, Liang T. Advantages of targeting the tumor immune microenvironment over blocking immune checkpoint in cancer immunotherapy. Signal Transduct Target Ther. 2021;6(1):72. Gocher AM, Workman CJ, Vignali DAA. Interferon-gamma: teammate or opponent in the tumour microenvironment? Nat Rev Immunol. 2022;22(3):158–72. Kim T, Croce CM. MicroRNA: trends in clinical trials of cancer diagnosis and therapy strategies. Exp Mol Med. 2023;55(7):1314–21. Holla S, Stephen-Victor E, Prakhar P, Sharma M, Saha C, Udupa V, et al. Mycobacteria-responsive sonic hedgehog signaling mediates programmed death-ligand 1- and prostaglandin E2-induced regulatory T cell expansion. Sci Rep. 2016;6:24193. Audrito V, Serra S, Stingi A, Orso F, Gaudino F, Bologna C, et al. PD-L1 up-regulation in melanoma increases disease aggressiveness and is mediated through miR-17-5p. Oncotarget. 2017;8(9):15894–911. Cheng G, Li Y, Liu Z, Song X. lncRNA PSMA3-AS1 promotes the progression of non-small cell lung cancer through targeting miR-17-5p/PD-L1. Adv Clin Exp Med. 2021;30(10):1043–50. Selem NA, Nafae H, Manie T, Youness RA, Gad MZ. Let-7a/cMyc/CCAT1/miR-17-5p Circuit Re-sensitizes Atezolizumab Resistance in Triple Negative Breast Cancer through Modulating PD-L1. Pathol Res Pract. 2023;248:154579. Du W, Chen D, Wei K, Yu D, Gan Z, Xu G, et al. MiR-10b-5p Impairs TET2-Mediated Inhibition of PD-L1 Transcription Thus Promoting Immune Evasion and Tumor Progression in Glioblastoma. Tohoku J Exp Med. 2023;260(3):205–14. Zhao P, Zhen H, Zhao H, Huang Y, Cao B. Identification of hub genes and potential molecular mechanisms related to radiotherapy sensitivity in rectal cancer based on multiple datasets. J Transl Med. 2023;21(1):176. Ye J, Zeng T. Mining database and verification of PIK3CB as a marker predicting prognosis and immune infiltration in renal clear cell carcinoma. Med (Baltim). 2022;101(22):e29254. Katoh Y, Yaguchi T, Kubo A, Iwata T, Morii K, Kato D et al. Inhibition of stearoyl-CoA desaturase 1 (SCD1) enhances the antitumor T cell response through regulating beta-catenin signaling in cancer cells and ER stress in T cells and synergizes with anti-PD-1 antibody. J Immunother Cancer. 2022;10(7). Li H, Yang Z, Yang X, Zhang F, Wang J, Wu Z, et al. LINC01123 promotes immune escape by sponging miR-214-3p to regulate B7-H3 in head and neck squamous-cell carcinoma. Cell Death Dis. 2022;13(2):109. You J, Wu W, Lu M, Xie Y, Miao R, Gu M, et al. Hepatic exosomes with declined MiR-27b-3p trigger RIG-I/TBK1 signal pathway in macrophages. Liver Int. 2022;42(7):1676–91. Marcinkowska M, Szymanski M, Krzyzosiak WJ, Kozlowski P. Copy number variation of microRNA genes in the human genome. BMC Genomics. 2011;12:183. Gulyaeva LF, Kushlinskiy NE. Regulatory mechanisms of microRNA expression. J Transl Med. 2016;14(1):143. Roy R, Chatterjee A, Das D, Ray A, Singh R, Chattopadhyay E, et al. Genome-wide miRNA methylome analysis in oral cancer: possible biomarkers associated with patient survival. Epigenomics. 2019;11(5):473–87. Veerappa AM, Murthy MN, Vishweswaraiah S, Lingaiah K, Suresh RV, Nachappa SA, et al. Copy number variations burden on miRNA genes reveals layers of complexities involved in the regulation of pathways and phenotypic expression. PLoS ONE. 2014;9(2):e90391. Chu G, Ji X, Wang Y, Niu H. Integrated multiomics analysis and machine learning refine molecular subtypes and prognosis for muscle-invasive urothelial cancer. Mol Ther Nucleic Acids. 2023;33:110–26. Wang B, Li M, Li R. Identification and verification of prognostic cancer subtype based on multi-omics analysis for kidney renal papillary cell carcinoma. Front Oncol. 2023;13:1169395. Liu B, Shyr Y, Cai J, Liu Q. Interplay between miRNAs and host genes and their role in cancer. Brief Funct Genomics. 2018;18(4):255–66. Louhimo R, Hautaniemi S. CNAmet: an R package for integrating copy number, methylation and expression data. Bioinformatics. 2011;27(6):887–8. Zhang F, Zhang Y, Lv X, Xu B, Zhang H, Yan J, et al. Evolution of an X-Linked miRNA Family Predominantly Expressed in Mammalian Male Germ Cells. Mol Biol Evol. 2019;36(4):663–78. Lin WW, Ou GY, Zhao WJ. Mutational profiling of low-grade gliomas identifies prognosis and immunotherapy-related biomarkers and tumour immune microenvironment characteristics. J Cell Mol Med. 2021;25(21):10111–25. Cheng Q, Duan W, He S, Li C, Cao H, Liu K, et al. Multi-Omics Data Integration Analysis of an Immune-Related Gene Signature in LGG Patients With Epilepsy. Front Cell Dev Biol. 2021;9:686909. Bartel DP. MicroRNAs: target recognition and regulatory functions. Cell. 2009;136(2):215–33. Quigley DA, Kristensen V. Predicting prognosis and therapeutic response from interactions between lymphocytes and tumor cells. Mol Oncol. 2015;9(10):2054–62. Paladini L, Fabris L, Bottai G, Raschioni C, Calin GA, Santarpia L. Targeting microRNAs as key modulators of tumor immune response. J Exp Clin Cancer Res. 2016;35:103. Pathania AS, Prathipati P, Olwenyi OA, Chava S, Smith OV, Gupta SC, et al. miR-15a and miR-15b modulate natural killer and CD8(+)T-cell activation and anti-tumor immune response by targeting PD-L1 in neuroblastoma. Mol Ther Oncolytics. 2022;25:308–29. Li X, Wang S, Mu W, Barry J, Han A, Carpenter RL, et al. Reactive oxygen species reprogram macrophages to suppress antitumor immune response through the exosomal miR-155-5p/PD-L1 pathway. J Exp Clin Cancer Res. 2022;41(1):41. Turcan S, Rohle D, Goenka A, Walsh LA, Fang F, Yilmaz E, et al. IDH1 mutation is sufficient to establish the glioma hypermethylator phenotype. Nature. 2012;483(7390):479–83. Additional Declarations No competing interests reported. Supplementary Files AdditionalFile1.pdf Additional file 1: A pdf document containing all Supplementary Figures. AdditionalFile2.xlsx Additional file 2: An Excel document containing all Supplementary Tables. Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-3935802","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":271584286,"identity":"2d8a5946-b891-45dd-a54b-0baa59d31fbe","order_by":0,"name":"Zhen Yang","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA80lEQVRIiWNgGAWjYDACCSjNz3wAyj5ArBbJtgRStRgcI1aL/OzmY495au7YbT7GY3jj5w4GOb4bCYyfC/BoYZxzLN2Y59iz5G3HeIwte88wGEveSGCWnoFHC7NEjpk0D9vhZLP7PWYSvG0MiRtuJLAx8+DRwiaR/02a59/hZOM2HjPJv20M9QS18EjksEnzth22M2DjMQMyGBIMCGmRkEgzk5zbdzhB4hhbsbVsm4ThzDMPm6XxaZGfkfxM4s23w/b8bcwbb75ts5HnO5588DM+LSDABFSQ2AC1FYgZGwhoACr5wcBgT1DVKBgFo2AUjFwAAPD0R4RA0Pk8AAAAAElFTkSuQmCC","orcid":"","institution":"Fudan University","correspondingAuthor":true,"prefix":"","firstName":"Zhen","middleName":"","lastName":"Yang","suffix":""},{"id":271584287,"identity":"aed9dd15-e7d8-47ca-afe6-9c2e0291ec0b","order_by":1,"name":"Xiaocen Liu","email":"","orcid":"","institution":"The First Affiliated Hospital of Wannan Medical College, Yijishan Hospital of Wannan Medical College)","correspondingAuthor":false,"prefix":"","firstName":"Xiaocen","middleName":"","lastName":"Liu","suffix":""},{"id":271584288,"identity":"e5fed56f-7c2e-4ab3-af25-7010c46de45c","order_by":2,"name":"Hao Xu","email":"","orcid":"","institution":"Fudan University","correspondingAuthor":false,"prefix":"","firstName":"Hao","middleName":"","lastName":"Xu","suffix":""},{"id":271584289,"identity":"25f1bfe5-95db-422f-b18d-d329076a0da3","order_by":3,"name":"Andrew E. Teschendorff","email":"","orcid":"","institution":"University of Chinese Academy of Sciences, Chinese Academy of Sciences","correspondingAuthor":false,"prefix":"","firstName":"Andrew","middleName":"E.","lastName":"Teschendorff","suffix":""},{"id":271584290,"identity":"76f63b60-931a-4c97-8d82-f55ffafeeedf","order_by":4,"name":"Lingjie Xu","email":"","orcid":"","institution":"West China hospital, Sichuan University","correspondingAuthor":false,"prefix":"","firstName":"Lingjie","middleName":"","lastName":"Xu","suffix":""},{"id":271584291,"identity":"86ee0e48-c3e6-41dc-a413-97bf89c21017","order_by":5,"name":"Jingyi Li","email":"","orcid":"","institution":"Capital Medical University","correspondingAuthor":false,"prefix":"","firstName":"Jingyi","middleName":"","lastName":"Li","suffix":""},{"id":271584292,"identity":"28e1188c-c968-421e-8ab6-c8ed39e4e615","order_by":6,"name":"Minjie Fu","email":"","orcid":"","institution":"Fudan University","correspondingAuthor":false,"prefix":"","firstName":"Minjie","middleName":"","lastName":"Fu","suffix":""},{"id":271584293,"identity":"f4f67e3c-b8cc-4fdb-abb7-77a9918a1777","order_by":7,"name":"Jun Liu","email":"","orcid":"","institution":"The First Affiliated Hospital of Wannan Medical College, Yijishan Hospital of Wannan Medical College)","correspondingAuthor":false,"prefix":"","firstName":"Jun","middleName":"","lastName":"Liu","suffix":""},{"id":271584294,"identity":"3890bc1f-e2ab-456e-9454-f82b5ac88277","order_by":8,"name":"Hanyu Zhou","email":"","orcid":"","institution":"Wannan Medical College","correspondingAuthor":false,"prefix":"","firstName":"Hanyu","middleName":"","lastName":"Zhou","suffix":""},{"id":271584295,"identity":"0e9f3b82-02b6-4531-8837-fc4b14b45976","order_by":9,"name":"Yingying Wang","email":"","orcid":"","institution":"The First Affiliated Hospital of Wannan Medical College, Yijishan Hospital of Wannan Medical College)","correspondingAuthor":false,"prefix":"","firstName":"Yingying","middleName":"","lastName":"Wang","suffix":""},{"id":271584296,"identity":"a87b5091-401f-441c-9b2f-538da90e65cb","order_by":10,"name":"Licheng Zhang","email":"","orcid":"","institution":"Fudan University","correspondingAuthor":false,"prefix":"","firstName":"Licheng","middleName":"","lastName":"Zhang","suffix":""},{"id":271584297,"identity":"f84b4529-abbd-4f47-a8cc-bae67050e96c","order_by":11,"name":"Yungang He","email":"","orcid":"","institution":"Shanghai Fifth People's Hospital, Fudan University","correspondingAuthor":false,"prefix":"","firstName":"Yungang","middleName":"","lastName":"He","suffix":""},{"id":271584298,"identity":"d2cdcbbd-1c06-446b-9cc9-6feb6fdcbd20","order_by":12,"name":"Kun Lv","email":"","orcid":"","institution":"Wannan Medical College","correspondingAuthor":false,"prefix":"","firstName":"Kun","middleName":"","lastName":"Lv","suffix":""},{"id":271584299,"identity":"d8aabcda-5c34-4445-a354-e65723f8c69e","order_by":13,"name":"Hui Yang","email":"","orcid":"","institution":"Wannan Medical College","correspondingAuthor":false,"prefix":"","firstName":"Hui","middleName":"","lastName":"Yang","suffix":""}],"badges":[],"createdAt":"2024-02-07 04:49:56","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-3935802/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-3935802/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":50937592,"identity":"21e03aa0-9ebf-4190-8693-80230c35b8f9","added_by":"auto","created_at":"2024-02-09 21:20:03","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":172820,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eCharacteristics of the copy number related miRNAs (CNV-miRs) and DNA methylation related miRNAs (DNAm-miRs) in LGG. \u003c/strong\u003eA). Distribution of Z-transformed correlation coefficients between miRNA expression level and copy number or DNA methylation across samples; B). Frequencies of CNV-miRs and DNAm-miRs across different chromosomes; C). Venn diagram shows the count of overlaps between CNV-miRs and DNAm-miRs; D-E). The proportional frequencies in each category of DNA methylation probes for DNAm-miRs and all miRNAs, respectively. Methylation probes are categorized based on the positional relations with transcription start site and CpG islands, respectively; F-G). GO terms for biological processes enriched by target genes of DNAm-miRs and CNV-miRs, respectively.\u003c/p\u003e","description":"","filename":"Figure1.png","url":"https://assets-eu.researchsquare.com/files/rs-3935802/v1/58260c6c0f5121f0a25d97c1.png"},{"id":50937590,"identity":"b343e066-a0c0-4f16-81f4-c97b2e170254","added_by":"auto","created_at":"2024-02-09 21:20:03","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":566330,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eIdentification of LGG molecular subtypes based on expression profiling of CNV-miRs and DNAm-miRs. \u003c/strong\u003eA). NMF clustering results for DNAm-miRs; B). NMF clustering results for CNV-miRs. Kaplan–Meier plot indicated differences for OS among subtypes identified by NMF clustering of DNAm-miRs (C) and CNV-miRs (D), respectively; E). Alluvial plot illustrates the overlap between the molecular subtypes identified from CNV-miRs and DNAm-miRs, as well as the histological subtypes; F). Sample number of the overlaps between the CNV-miR and the DNAm-miR subtypes.\u003c/p\u003e","description":"","filename":"Figure2.png","url":"https://assets-eu.researchsquare.com/files/rs-3935802/v1/8b957b592b77f9182951b495.png"},{"id":50937591,"identity":"dcc5dcba-24c1-4e8a-bf1f-71b52f11caf9","added_by":"auto","created_at":"2024-02-09 21:20:03","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":2053690,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eLGG sample classification by integrative analysis of multi-omics data.\u003c/strong\u003e A). Comprehensive clustering heatmaps for miRNA expression, DNA methylation and copy number variation identify consensus ensemble subtypes of LGG. Samples are labeled for the subtypes identified for CNV-miR and the DNAm-miR and also the histological subtypes; B). Kaplan–Meier plot indicated differences for OS among ensemble subtypes; C). P-values generated by log-rank test that among different integrative clusters; D). Kaplan–Meier plot indicated differences for OS among histological subtypes; E). Sample number of the overlaps between the integrative clusters and the histological subtypes.\u003c/p\u003e","description":"","filename":"Figure3.png","url":"https://assets-eu.researchsquare.com/files/rs-3935802/v1/398ba09351314ac68a40255e.png"},{"id":50937600,"identity":"2184a136-a4f1-4a66-ba14-9f4ec4d6d68d","added_by":"auto","created_at":"2024-02-09 21:20:04","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":315406,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eThe landscape of genomic features of LGG subtypes.\u003c/strong\u003eA). Oncoplot showing the differentially mutated genes among the subtypes of LGG; B). Alluvial plot illustrates the overlap between the LGG subtypes and mutational status of \u003cem\u003eIDH1\u003c/em\u003e, \u003cem\u003eTP53\u003c/em\u003e and \u003cem\u003eEGFR\u003c/em\u003e, respectively; C). Boxplot showing the distribution of tumor mutation burden among different clusters; D). Barplot of the pathway alteration frequency and the fraction of sample affected for IC2 and other clusters; E). Potentially druggable gene categories from mutation datasets for IC2; F). Potentially druggable gene categories from mutation datasets for other clusters.\u003c/p\u003e","description":"","filename":"Figure4.png","url":"https://assets-eu.researchsquare.com/files/rs-3935802/v1/6de726f8d8852fdec815bc1a.png"},{"id":50937597,"identity":"d68ac36b-b4c2-45ca-8f3d-0b946c8ee604","added_by":"auto","created_at":"2024-02-09 21:20:03","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":536590,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eCharacteristics of tumor immunophenotypes in LGG subgroups. \u003c/strong\u003eA). Bar plot showing immune related pathways enriched by genes differentially expressed IC2 and other clusters, significant threshold at p \u0026lt; 0.05 are indicated; B). Boxplot showing the distribution of the ESTIMATE immune score across different clusters; C). Radar plot showing the mean value of scaled antitumor immunoactivity and immune signature scores across different clusters; D). Violin plots showing the immune cell infiltration levels obtained from TIMER for different clusters; E). Heatmap showing the relative expression of immune cell type marker genes in different clusters.\u003c/p\u003e","description":"","filename":"Figure5.png","url":"https://assets-eu.researchsquare.com/files/rs-3935802/v1/1911ac8bb10983791da93e6d.png"},{"id":50938088,"identity":"d4e827d5-2ccf-4aa4-92c3-e142174f0edd","added_by":"auto","created_at":"2024-02-09 21:28:04","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":340036,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003emiRNAs correlated with immune cell infiltrations in LGG.\u003c/strong\u003eA). Bar plots showing the number of DNAm-miRs and CNV-miRs in which the expression correlated with immune cell infiltrations, respectively. The number of positive and negative correlated miRNAs are indicated by yellow and purple, respectively; B). the number of miRNAs that correlated with different numbers of immune cell types; C). Heatmaps showing the number of miRNAs in which expression correlated with different types of immune cell infiltrations within different clusters; D). Venn diagram indicated the number of specific and shared miRNAs that related to immune cell and that from sections of “Immune Function” and “Cancer Immunology” of RNA2Immune; E). Bar plots shows the number of immune miRNAs from DNAm-miR and CNV-miR, respectively; F). Bubble plots of correlation between expression and immune cell infiltration for prognosis related miRNAs; G). Risk score and expression heatmap of the eight signature miRNAs in the training set; H). Kaplan–Meier plot of survival for patients with different risk scores from an independent validation cohort of CGGA; I). Time dependent ROC curve and AUC indicated the prognostic performance for survival prediction in validation set.\u003c/p\u003e","description":"","filename":"Figure6.png","url":"https://assets-eu.researchsquare.com/files/rs-3935802/v1/d92dfee6b237490c39fc95b9.png"},{"id":50937593,"identity":"8f82b3ad-d9b4-42dd-88f4-264d467e98b8","added_by":"auto","created_at":"2024-02-09 21:20:03","extension":"png","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":290576,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eRelationship between the risk model and other clinical information.\u003c/strong\u003e A). qRT-PCR validation for the expression of eight signature miRNAs of our in-house cohort; B). Univariate and multivariate regression analysis of the relation between our prognostic model and other clinicopathological features regarding prognostic value; C). Nomogram for predicting the probability of 1-, 2-, 3- and 5-year survival for LGG patients; D). Calibration plot of the nomogram for predicting the probability of patient’s survival at 1-, 2-, 3- and 5 years; E). Time-dependent ROC curves based on the nomogram for 1-, 2-, 3- and 5 years survival probability.\u003c/p\u003e","description":"","filename":"Figure7.png","url":"https://assets-eu.researchsquare.com/files/rs-3935802/v1/208cf94337e01786c6b3ffa4.png"},{"id":50937601,"identity":"b368d8ce-1214-4b56-9ff1-c0205b191e56","added_by":"auto","created_at":"2024-02-09 21:20:04","extension":"png","order_by":8,"title":"Figure 8","display":"","copyAsset":false,"role":"figure","size":538442,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eBiological function and immune evasion of LGG subtypes. \u003c/strong\u003eA). Heatmap of the enrichment score obtained from ssGSEA on Hallmark genesets; B). Violin plots showing the distribution of “T−score”, “I−score”, and “S−score” among different clusters; C). Boxplot showing the distribution of PD-L1 expression among different clusters; D). Scatter plot showing significant correlation between expression of PD-L1 and miR-338-5p; E). Scatter plot showing significant correlation between expression and promoter DNA methylation of miR-338-5p; F). Regulatory module consists of PD-L1 and the associated miRNA regulators.\u003c/p\u003e","description":"","filename":"Figure8.png","url":"https://assets-eu.researchsquare.com/files/rs-3935802/v1/856d1f5e452ce40406f3b3bd.png"},{"id":50937596,"identity":"32baa2e5-3cd8-480d-9066-e87d637373c3","added_by":"auto","created_at":"2024-02-09 21:20:03","extension":"png","order_by":9,"title":"Figure 9","display":"","copyAsset":false,"role":"figure","size":503937,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eIdentification of miRNA immunoevasive biomarkers.\u003c/strong\u003eA). Violin plots showing the distribution of TIDE score among different clusters; B). Stacked bar chart of the fractions of immune therapy response for patients among different clusters; C). Scatter plot showing significant correlation between TIDE score and miR-10b-5p expression; D) miRNA immune evasive network for the putative targets and immune evasive miRNAs in LGG.\u003c/p\u003e","description":"","filename":"Figure9.png","url":"https://assets-eu.researchsquare.com/files/rs-3935802/v1/1df209fb5053314819d970c6.png"},{"id":50937598,"identity":"191c6d6d-9833-450a-8ffe-7dc69120f8c6","added_by":"auto","created_at":"2024-02-09 21:20:03","extension":"png","order_by":10,"title":"Figure 10","display":"","copyAsset":false,"role":"figure","size":2646003,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eRole of the prognosis related miRNAs in the proliferation, migration, invasion and immune evasion of glioma cell lines.\u003c/strong\u003eA). The expression of miRNAs (miR-10b-5p, miR-155-5p, miR-196a-5p, miR-196b-5p, miR-200a-3p, miR-204-5p, miR-503-5p, and miR-15b-5p) in U251 cells was determined by qRT-PCR; B). CCK-8 assays show that the inhibition of miRNAs (miR-155-5p, miR-196a-5p, miR-196b-5p, miR-200a-3p, miR-503-5p, and miR-15b-5p) decreased cell proliferation in U251 cell lines cells; C-D). Wound healing assays show that miRNAs (miR-155-5p, miR-196a-5p, miR-196b-5p, miR-200a-3p, miR-503-5p, and miR-15b-5p) knockdown significantly reduced the cell migration ability of U251 cells with the representative images and the quantitative analysis; E-F). Transwell invasion of miRNAs (miR-155-5p, miR-196a-5p, miR-196b-5p, miR-200a-3p, miR-503-5p, and miR-15b-5p) knockdown cells is significantly reduced compared with control cells; G). The expression of \u003cem\u003eFOXP3\u003c/em\u003e, \u003cem\u003eCTLA-4\u003c/em\u003e and \u003cem\u003ePD-L1\u003c/em\u003e in U251 cells transduced with miRNAs (miR-10b-5p, miR-155-5p, miR-196a-5p, miR-196b-5p, miR-200a-3p, miR-204-5p, miR-503-5p, and miR-15b-5p) inhibitor was examined by Western blot; β-actin was an internal control. Data are the mean ± SD from three independent experiments. * P \u0026lt; 0.05.\u003c/p\u003e","description":"","filename":"Figure10.png","url":"https://assets-eu.researchsquare.com/files/rs-3935802/v1/f5d76e102a7c6e1da01c0ebe.png"},{"id":50974808,"identity":"60e9f7c5-05df-4851-87a4-d4cfd640555a","added_by":"auto","created_at":"2024-02-12 00:40:16","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":6540279,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-3935802/v1/ab2c01db-40d9-45fb-9437-ea1536f067c8.pdf"},{"id":50937595,"identity":"4c50110f-ee01-4076-acc5-c1b022f22d28","added_by":"auto","created_at":"2024-02-09 21:20:03","extension":"pdf","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":5192664,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eAdditional file 1: \u003c/strong\u003eA pdf document containing all Supplementary Figures.\u003c/p\u003e","description":"","filename":"AdditionalFile1.pdf","url":"https://assets-eu.researchsquare.com/files/rs-3935802/v1/712ac8b91dce8135bd4dbb87.pdf"},{"id":50938087,"identity":"9d9d9770-8018-4775-a209-ca43d573c678","added_by":"auto","created_at":"2024-02-09 21:28:03","extension":"xlsx","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":467352,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eAdditional file 2: \u003c/strong\u003eAn Excel document containing all Supplementary Tables.\u003c/p\u003e","description":"","filename":"AdditionalFile2.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-3935802/v1/5c2817622e12d2dcac51227e.xlsx"}],"financialInterests":"No competing interests reported.","formattedTitle":"Integrative analysis of genomic and epigenomic regulation reveals microRNA regulatory network mediated tumor heterogeneity and immune evasion in lower grade glioma","fulltext":[{"header":"Background","content":"\u003cp\u003eLower grade glioma (LGG) is the most common and aggressive central nervous system (CNS) tumor, which accounting for about 40% of all brain malignancies (\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e). Most LGGs can be further classified according to specific type of cell with which they share histological features or classic molecular markers, such as isocitrate dehydrogenase (IDH) mutation or O\u003csup\u003e6\u003c/sup\u003e-methylguanine-DNA methyltransferase (MGMT) promoter methylation (\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e, \u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e). Despite tremendous progress has been made over the past decades for its early detection, surgical paradigm and multidisciplinary treatment, such as neoadjuvant chemotherapy and radiotherapy, the postoperative overall survival of LGG patients remains extremely low, particularly for high-grade glioma with worse prognoses for their malignant aggressivity (\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e). In general, the prognosis of patients with glioma varies dramatically, which is dependent on different tumor grades, key gene alternation status, tumor microenvironment (TME), and combination of different efficacious method (\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e). Recent genomic studies have brought a comprehensive view of glioma based on molecular profiling and identified different molecular subtypes within the glioma that appear to correlate with disease etiology and therapy response (\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e). Hence, there is an urgent need to characterize more specific and precise molecular signatures of LGG for its accurate diagnosis, individualized treatment and the prognosis assessment.\u003c/p\u003e \u003cp\u003eMicroRNAs (miRNAs) are ~\u0026thinsp;22 nt small noncoding RNAs that play an important role in post-transcriptional regulation, which precisely adjust the gene expression level by targeting mRNAs for its degradation or translational repression (\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e, \u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e). MiRNAs are critical for normal development of organisms and are involved in a variety of biological processes including cell cycle, cell proliferation, differentiation, apoptosis and cellular signaling (\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e, \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e). Aberrant expression of miRNAs is associated with many human diseases (\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e, \u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e). For glioma, miRNA expression has been associated with the clinical and molecular classification (\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e), progression (\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e, \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e), prognosis (\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e, \u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e) and response to targeted therapies (\u003cspan additionalcitationids=\"CR19\" citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e). Altered miRNA regulation is widely involved in glioma pathogenesis via the modulation of oncogenes and the associated downstream signaling pathways. Till now, most studies were proposed concerning aberrant expression and functional curation of miRNAs in cancer and disease, but there still lacks systematic investigation of the mechanisms underlying miRNA regulatory activities, especially the genome-wide studies for genetic and epigenetic regulation of miRNAs in development and progression of glioma are still scarce.\u003c/p\u003e \u003cp\u003eIn recent years, large multi-omics studies greatly enhanced our understanding of disease dysregulation at genomic and epigenetic levels. Genome wide alternations such as copy number variations (CNVs) and DNA methylations (DNAm) are widely observed during tumorigenesis, which promote cancer progression (\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e, \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e). Copy number variation refers to the number of copies of a specific DNA segment in the genome varies among individuals (\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e). As a key genomic regulator that contribute to gene expression dysregulation via modulating mRNA levels and by influencing transcriptional regulation, CNV is significantly correlated with certain pathological processes such as cancer (\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e). Transcriptional disorders caused by CNV changes have been identified as driving events in LGG progression. In addition, gene expression regulation via DNA methylation at the epigenetic level also exerts a crucial part in the behaviors of various cancers (\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e, \u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e). DNAm imbalances is able to affect the occurrence and development of cancer by defining different types of driver events, such as cell growth, proliferation, differentiation, and apoptosis processes (\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e, \u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e). Genomic profiling studies indicated that these genomic and epigenomic dysregulations are highly heterogeneous and jointly regulate the occurrence and evolution of tumor (\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e). Till now, the influence of DNAm and CNV on noncoding RNA, especially miRNA expression remains largely unexplored. Studies indicated various regulators are required for precise control the expressions of miRNAs (\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e). Therefore, miRNAs with aberrant expression under genetic and epigenetic contexts may serve as valuable marker for investigation of the complex regulatory mechanisms in tumor progression, and a systematic analysis of DNAm and CNV patterns of miRNAs in large scale cancer samples is critically needed.\u003c/p\u003e \u003cp\u003eIn this study, we collected copy number variation, DNA methylation and miRNA expression profiles from the cohort of The Cancer Genome Atlas (TCGA) LGG patients, and identified miRNAs whose expression levels are regulated by genomic and/or epigenomic deregulation. In addition, by using the multi-omics integration study, we identified different molecular subtypes of LGG that are significantly associated with prognosis. Furthermore, specific miRNA biomarkers were proposed to drive the classification of these subtypes, which were validated in independent cohort data base on different platforms. By assessing the associations between immune cell infiltrations and immune pathways, we identified a subset of immune-related miRNAs, and several miRNAs are identified as immune evasion markers with high confidence. Our study developed novel strategies of multi-omics integration and characterized the genetic and epigenetic landscape of genes encoding miRNAs, and also provides a framework for identifying novel biomarkers to guide the clinical treatment for glioma patients.\u003c/p\u003e"},{"header":"Material and Methods","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003eData collection and preprocessing\u003c/h2\u003e \u003cp\u003eThe genome wide profiling data for LGG was download from TCGA (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://portal.gdc.cancer.gov/\u003c/span\u003e\u003cspan address=\"https://portal.gdc.cancer.gov/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) (\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e). Here we obtained 512 samples with miRNA expression, 505 samples with DNA methylation and 515 samples with CNV detection, a total of 500 samples with all three different types of high-throughput profiling data were retained for further integrative study. The genome-wide miRNA profiling data which based on miRNA-seq was quantified as RPM (reads per million mapped reads). The expression value was log2(RPM\u0026thinsp;+\u0026thinsp;1) transformed in order to regularize the data. For DNA methylation data based on Illumina Infinium HumanMethylation450 BeadChip array, methylation level of each probe was measured as the beta-value. We removed the probes whose beta-values are missing in more than 30% of the samples, the remaining probes with missing values were imputed by using the k-nearest neighbors (KNN) method (\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e). Finally, the BMIQ method was used to correct for the type II probe bias (\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e). For CNV profile, segmentation data generated by Affymetrix SNP 6.0 platform were used. The \u0026lsquo;nocnv.seg\u0026rsquo; file contains chromosome location and segment mean value (log2(copy-number/2)) for each sample and was used to capture the focal somatic CNV. In addition, we also collected RNA-seq data for 500 LGG samples which was quantified as logarithm transformed (base 2) FPKM (Fragments Per Kilobase per Million), as well as the somatic mutation data for 508 samples generated by MuTect2 pipeline. Furthermore, clinical information of all samples including age, gender, survival status, histological type, tumor stage and overall survival time was also retrieved from the TCGA data portal.\u003c/p\u003e \u003cp\u003eAn external dataset of 198 glioma patients with miRNA expression profiling data was downloaded from the Chinese Glioma Genome Atlas (CGGA) (\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e). This miRNA profiling data is based on the Illumina Human v2 MicroRNA Expression BeadChip. We transfer the miRNA ID from the microarray data to that from miRBase v21 by using the \u0026ldquo;ID convertor\u0026rdquo; tool in dbDEMC v3.0 (\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e). The expression values were logarithmically transformed (base 2) and quantile normalized. The associated clinical data which includes survival status and overall survival time of the patients was also obtained.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec4\" class=\"Section2\"\u003e \u003ch2\u003ePatients in the in-house cohort and sample collection\u003c/h2\u003e \u003cp\u003eHuman tissue samples of 30 pairs of adjacent non-tumor and glioma samples were collected from the Department of Neuro-surgery of the First Affiliated Hospital, Wannan Medical College (Wuhu, Anhui, P. R. China) and Department of Neuro-surgery of Huashan Hospital, Fudan University (Shanghai, P. R. China) from February 2023 to September 2023. Patients received no systemic treatment before the surgery. All samples obtained at the surgery were directly preserved in liquid nitrogen. Two pathologists evaluated all specimens according to the WHO guidelines. Patient consent was obtained. All surveys and experiments were approved by the Ethic Committee for Clinical Research of the First Affiliated Hospital of Wannan Medical College and Huashan Hospital, Fudan University, respectively.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec5\" class=\"Section2\"\u003e \u003ch2\u003eDNA copy numbers profiling construction for miRNAs\u003c/h2\u003e \u003cp\u003eTo construct copy number profiling for miRNAs in LGG, we investigated whether the positions of miRNA loci can be mapped to genomic regions with copy number in a certain sample. We first obtained the genomic coordinates for mature miRNAs via miRBase version 21 (\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e). Then segment-mean values in distinct genomic loci (genome assembly hg38) were extracted from TCGA SNP array data. We assigned the locus of each miRNA to the copy number values by using the GISTIC2 (\u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e). In this way, the copy number profiling data for individual estimates were obtained for a total of 2265 miRNAs.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec6\" class=\"Section2\"\u003e \u003ch2\u003emiRNA promoter identification and Infinium 450K Array probe reannotation to construct miRNA methylation profiles\u003c/h2\u003e \u003cp\u003eTo characterize DNA methylation patterns for miRNAs, we employed a novel strategy to reannotate the probes of Infinium 450K arrays into miRNA-associated promoter regions. We first obtained genome coordinates for miRNA promoters from the fifth edition of the Functional Annotation of Mammalian Genome database (FANTOM5) (\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e). FANTOM5 incorporates miRNA transcription start site (TSS) information generated by Cap Analysis Gene Expression (CAGE) methods, which was based on genome assembly hg19. The coordinates of miRNA TSS were converted to genome assembly hg38 by using the liftOver tool from the UCSC genome browser. As a result, a total of 2050 mature miRNAs are identified for unique physical genomic location of TSS among 2587 miRNAs from miRBase v21 (Table \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e). Here we define the promoter region for individual miRNA as between 2000 bp upstream to 500 bp downstream of TSS, which is akin to protein coding genes. By matching the genomic coordinates of the miRNA TSS with the 450K Array probe locations, we filtered 8825 cpg sites that are located within miRNA promoters (Table \u003cspan refid=\"MOESM2\" class=\"InternalRef\"\u003eS2\u003c/span\u003e). Then we assign the average methylation value of the probes within promoter region for each miRNA to obtain its final methylation value. Finally, we obtained methylation profile for 1977 miRNAs for LGG.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec7\" class=\"Section2\"\u003e \u003ch2\u003eIdentification of CNV and DNAm related miRNA sets\u003c/h2\u003e \u003cp\u003eThe Pearson correlation coefficients (r) were calculated for paired profiles of CNV and miRNA-seq and between DNAm and miRNA-seq for each miRNA, respectively, the correlation coefficient was then transformed to the Fisher Z-scale according to the formula:\u003cdiv id=\"Equa\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equa\" name=\"EquationSource\"\u003e\n$$FisherZ=\\frac{1}{2}ln\\frac{1+r}{1-r}$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003eThe miRNA of FDR adjusted p-value\u0026thinsp;\u0026lt;\u0026thinsp;0.05 with correlation test constitutes the miRNA set significantly positively related to CNV (CNV-miRs) and the miRNA set negatively related to methylation (DNAm-miRs), respectively.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec8\" class=\"Section2\"\u003e \u003ch2\u003eClustering analysis of multi-layered genomic profiles\u003c/h2\u003e \u003cp\u003eNonnegative matrix factorization (NMF) is an unsupervised clustering method which was widely used in tumor subtypes identification based on omics data. Here we employ NMF to identify stable sample clusters by using expression of CNV-miRs and DNAm-miRs, respectively. Specifically, cluster analysis with standard \u0026ldquo;brunet\u0026rdquo; method and 50 iterations was used, we set the number of cluster k as 2 to 7, the optimal cluster number k was then determined based on the observed consensus map together with cophenetic correlation between clusters. At the same time, the mean silhouette width was computed by R package \u0026ldquo;NMF\u0026rdquo; to examine consensus membership matrix (\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e). For integrative cluster analysis of the multi-omics data for cancer subtyping, we applied \u0026ldquo;MOVICS\u0026rdquo; R package to integrate the CNV data of the CNV-miRs, methylation data of the DNAm-miRs, as well as the expression profile data of both CNV-miRs and DNAm-miRs (\u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e). MOVICS provides an interface for 10 state-of-the-art multi-omics data clustering algorithms (SNF, PINSPlus, NEMO, COCA, LRAcluster, CIMLR, Consensus clustering, IntNMF, MoCluster, iClusterBayes) and combines the output of each algorithm to obtain more robust classification. The performance of these algorithms had been validated by previous researches (\u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e). We calculated the clustering prediction index to find the optimal number of clusters and the consensus matrix was used to validate robust clustering of the samples.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec9\" class=\"Section2\"\u003e \u003ch2\u003eMeasurement for immune cell infiltration and antitumor immunoactivity\u003c/h2\u003e \u003cp\u003eIn this study, we identified and evaluated the levels of immune infiltrates by the Tumor immune estimation resource (TIMER) algorithm, which is a resource for systematic estimation of the abundances of six immune cell types (B cell, CD4\u0026thinsp;+\u0026thinsp;T cell, CD8\u0026thinsp;+\u0026thinsp;T cell, neutrophil, macrophage and dendritic cell) by using transcriptome profiles (\u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e). The abundances of these cells were compared across different molecular subtypes for quantitative characterization of the tumor immune microenvironment. We used ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) to assess immune activity of the tumor samples. ESTIMATE yields the immune score and stromal score for each tumor sample that quantifies the levels of infiltration levels of immune and stromal cell based on its gene expression profiles (\u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e43\u003c/span\u003e). We collected 17 immune relevant pathways and the associated 1,811 genes from ImmPort (\u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e). In addition, three different indices were constructed to measure the immune activities for tumor samples: Cytolytic activity (CYT) score, Major Histocompatibility Complex (MHC) score and Cytotoxic T Lymphocyte (CTL) score. The CYT score was defined as the geometric mean of \u003cem\u003eGZMA\u003c/em\u003e and \u003cem\u003ePRF1\u003c/em\u003e:\u003cdiv id=\"Equb\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equb\" name=\"EquationSource\"\u003e\n$${Score}_{CYT}=\\sqrt{{Exp}_{GZMA}{*Exp}_{PRF1}}$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003eThese two genes are found upregulated in activated CD8\u0026thinsp;+\u0026thinsp;T cells and play important roles in cytolysis and response to \u003cem\u003ePD-1\u003c/em\u003e and \u003cem\u003eCTLA4\u003c/em\u003e related immunotherapy (\u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e45\u003c/span\u003e, \u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e46\u003c/span\u003e). The MHC score of each sample was calculated as the average expression as nine genes of \u003cem\u003eHLA-A\u003c/em\u003e, \u003cem\u003ePSMB9\u003c/em\u003e, \u003cem\u003eHLA-B\u003c/em\u003e, \u003cem\u003ePSMB8\u003c/em\u003e, \u003cem\u003eHLA-C\u003c/em\u003e, \u003cem\u003eB2M\u003c/em\u003e, \u003cem\u003eTAP2\u003c/em\u003e, \u003cem\u003eNLRC5\u003c/em\u003e and \u003cem\u003eTAP1\u003c/em\u003e:\u003cdiv id=\"Equc\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equc\" name=\"EquationSource\"\u003e\n$${Score}_{MHC}=\\sum {Exp}_{n}/9$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003eThese genes were identified as the core gene set of MHC-I and are associated with activity of MHC-I antigen presentation (\u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e47\u003c/span\u003e). Finally, the CTL score for each sample was calculated as average expression from five genes of \u003cem\u003eCD8A\u003c/em\u003e, \u003cem\u003eCD8B\u003c/em\u003e, \u003cem\u003eGZMA\u003c/em\u003e, \u003cem\u003eGAMB\u003c/em\u003e and \u003cem\u003ePRF1\u003c/em\u003e:\u003cdiv id=\"Equd\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equd\" name=\"EquationSource\"\u003e\n$${Score}_{CTL}=({{Exp}_{CD8A}+{Exp}_{CD8B}+Exp}_{GZMA}+{Exp}_{GZMB}+{Exp}_{PRF1})/5$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003eThese five genes were important factors to measure cytotoxicity of tumor-infiltrating T cells and immune cell effector function (\u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e48\u003c/span\u003e). In addition, we also obtained the gene sets that represents different immune signatures from several publications, including the immune check point, human leukocyte antigen (HLA), interferon (IFN) response and tumor-infiltrating lymphocytes (TILs), which are collected from publication of Ju \u003cem\u003eet al.\u003c/em\u003e (\u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e49\u003c/span\u003e), and the gene sets for Macrophage/Monocyte, TGF-β response, IFN-γ response, and Wound healing collected from Vesteinn \u003cem\u003eet al\u003c/em\u003e. (\u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e50\u003c/span\u003e). We use the single-sample Gene Set Enrichment Analysis (ssGSEA) algorithm to estimate the activity of tumor samples by the immune signatures obtained (\u003cspan citationid=\"CR51\" class=\"CitationRef\"\u003e51\u003c/span\u003e). Finally, another 63 immune cell type specific gene markers for 21 immune cell populations were collected from TISIDB (\u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e52\u003c/span\u003e).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec10\" class=\"Section2\"\u003e \u003ch2\u003eImmune evasion and immunotherapy response evaluation\u003c/h2\u003e \u003cp\u003eTo evaluate complex interactions among cancer cells and immune microenvironment, we constructed the scoring system following Shuichi \u003cem\u003eet al\u003c/em\u003e. to describe factors that relevant to the tumor proliferation (T-factor), antitumor immunity (I-factor) and immunosuppression (S-factor) (\u003cspan citationid=\"CR53\" class=\"CitationRef\"\u003e53\u003c/span\u003e). Briefly, nine representative gene sets related to those factors were collected from the molecular signature database (MsigDB V.7.1) (\u003cspan citationid=\"CR54\" class=\"CitationRef\"\u003e54\u003c/span\u003e), which include Gene Ontology (GO) terms of \u0026ldquo;SPINDLE_LOCALIZATION\u0026rdquo; and \u0026ldquo;POSITIVE_REGULATION_OF_GLUCOSE_METABOLIC_PROCESS\u0026rdquo; for T-factor, \u0026ldquo;T cell\u0026rdquo;, \u0026ldquo;B cell\u0026rdquo; Signatures and GO terms of \u0026ldquo;MHC_PROTEIN_COMPLEX\u0026rdquo; for \u0026ldquo;I-factor\u0026rdquo;, and \u0026ldquo;NEGATIVE_REGULATION_OF_CYTOKINE_PRODUCTION_INVOLVED_IN_INFLAMMATORY_RESPONSE\u0026rdquo;, \u0026ldquo;POSITIVE_REGULATION_OF_CELL_MIGRATION_INVOLVED_IN_SPROUTING_ANGIOGENESIS\u0026rdquo;, \u0026ldquo;NEGATIVE_REGULATION_OF_HYPOXIA_INDUCED_INTRINSIC_APOPTOTIC_SIGNALING_PATHWAY\u0026rdquo; as well as \u0026ldquo;EPITHELIAL_TO_MESENCHYMAL_TRANSITION\u0026rdquo; for \u0026ldquo;S-factor\u0026rdquo;, the individual ssGSEA enrichment score was calculated and z-score normalized, the mean values of z-scores for the gene sets allocated to T-factors, I-factors or S-factors and scored was calculated to obtain T-scores, I-scores or S-scores. Finally, we apply the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm to estimate the immunotherapy response of LGG patients (\u003cspan citationid=\"CR55\" class=\"CitationRef\"\u003e55\u003c/span\u003e). Higher TIDE scores implicated a higher potential of tumor immune evasion, thus less likely to benefit from immunotherapy.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec11\" class=\"Section2\"\u003e \u003ch2\u003eImmune-related miRNAs identification and prognosis model development\u003c/h2\u003e \u003cp\u003eTo identify immune regulatory miRNAs involved in LGG subtyping, Pearson correlation analysis was performed between miRNAs expression and levels of immune cell infiltrating. MiRNAs with FDR less than 0.05 were considered as significant. We further collected the immune related miRNAs from the RNA2Immune database for functional validation (\u003cspan citationid=\"CR56\" class=\"CitationRef\"\u003e56\u003c/span\u003e). This database provides a high-quality experimentally supported resource that links ncRNA regulatory mechanisms to four sections of \u0026ldquo;immune cell function\u0026rdquo;, \u0026ldquo;immune diseases\u0026rdquo;, \u0026ldquo;cancer immunity\u0026rdquo;, and vaccines. We selected miRNAs from the \u0026ldquo;immune cell function\u0026rdquo;, and \u0026ldquo;cancer immunity\u0026rdquo; for validation propose. To find the miRNAs with the height prognostic value, we applied Cox-proportional hazards analysis based on the Least absolute shrinkage and selection operator (LASSO) estimation, which have been identified that suitable for model construction when dealing with collinearity problem of covariates. Among the immune miRNAs that were significant correlated with immune cell infiltrating, key immune miRNAs were selected by performing the LASSO-penalized Cox regression, which was implemented by the glmnet R package (version 4.1-7). Finally, an immune miRNA-related prognostic model was constructed by utilizing the regression coefficients derived from Cox regression analysis to multiply the normalized expression level of each immune miRNA. The \u0026ldquo;surv_cutpoint\u0026rdquo; function within the \u0026ldquo;survminer\u0026rdquo; package was applied to determine the best cutoff to classify patients into high-risk and low-risk groups. The Kaplan-Meier survival analysis and log-rank test were used to assess the predictive ability of the prognostic model.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec12\" class=\"Section2\"\u003e \u003ch2\u003eValidation of the prognostic model from independent cohorts\u003c/h2\u003e \u003cp\u003eTo validate whether the predictions of the prognostic model were independent of other clinical features, we collected miRNA profiling data form CGGA and calculated the risk scores for the patients by using the model constructed. In addition, the predictive accuracies were compared using the time-dependent receiver operating characteristic (ROC) analyses and the associated Area Under Curve (AUC) Score.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec13\" class=\"Section2\"\u003e \u003ch2\u003eConstruction and evaluation of the nomogram\u003c/h2\u003e \u003cp\u003eWe assessed the correlation between prognostic model and other clinicopathological features. The Cox proportional hazards model was used to perform standard univariate and multivariate analysis. Prediction error curves were used to compare the accuracy of survival models. The Cox regression coefficients were used to construct the nomogram. Calibration plots were generated to explore the performance characteristics of the nomogram. Time dependent ROC analysis was performed to assess the predictive accuracy of the nomogram. Decision curve analysis was used to assess the clinical practicability of the nomogram. The statistical significance was set at 0.05.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec14\" class=\"Section2\"\u003e \u003ch2\u003emiRNA target prediction and functional enrichment analysis\u003c/h2\u003e \u003cp\u003eFor both miRNAs and protein coding genes that differentially expressed compared between groups, the well-known empirical Bayes framework of moderated t-statistics were employed for assessing differential expression in microarray experiments (\u003cspan citationid=\"CR57\" class=\"CitationRef\"\u003e57\u003c/span\u003e). miRNA target prediction was obtained from databases of TargetScan (\u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e58\u003c/span\u003e) and miRDB (\u003cspan citationid=\"CR59\" class=\"CitationRef\"\u003e59\u003c/span\u003e). The miRNA-target interactions predicted by both of the two algorithms were identified as valid ones. In addition, we also collected experimentally validated miRNA targets form miRTarBase (\u003cspan citationid=\"CR60\" class=\"CitationRef\"\u003e60\u003c/span\u003e). For functional annotation of the miRNA targets involved in the regulatory network, The GO enrichment analysis was performed by using clusterProfiler (\u003cspan citationid=\"CR61\" class=\"CitationRef\"\u003e61\u003c/span\u003e), the enriched GO terms were obtained at the threshold of FDR\u0026thinsp;\u0026lt;\u0026thinsp;0.05. In addition, hallmark gene sets were obtained from MsigDB, and Gene Set Variation Analysis (GSVA) analysis was used to infer the degree of activity of gene sets/pathways in individual sample\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec15\" class=\"Section2\"\u003e \u003ch2\u003eDistinct genomic features across different groups\u003c/h2\u003e \u003cp\u003eThe tumor mutation burden (TMB) and oncoplot of somatic mutation data were calculated and generated by using R package \u0026ldquo;maftools\u0026rdquo; (\u003cspan citationid=\"CR62\" class=\"CitationRef\"\u003e62\u003c/span\u003e). The different mutational frequencies for individual genes between groups were calculated by chi-square test, the p-value that less than 0.05 were identified as having achieved statistical significance. In addition, the druggable pathways and categories are also obtained by enrichment analysis.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec16\" class=\"Section2\"\u003e \u003ch2\u003eQuantitative real-time polymerase chain reaction (qRT\u0026ndash;PCR)\u003c/h2\u003e \u003cp\u003eTotal RNA was extracted from tissues with TRIzol reagent from 30 pairs of LGG and normal tissues according to the manufacturer\u0026rsquo;s protocol. Reverse transcription was performed using the PrimeScript RT Reagent Kit. The quantitative PCR assay was performed using the SYBR Green PCR kit. The expression of eight miRNAs of the in-house cohort was normalized by GAPDH, acting as an internal control. Relative miRNA expression levels were calculated by the 2\u003csup\u003e\u0026minus;ΔΔCT\u003c/sup\u003e method. Each sample was repeated in triplicate and analyzed using relative quantification method.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec17\" class=\"Section2\"\u003e \u003ch2\u003eCell culture and transfection\u003c/h2\u003e \u003cp\u003eThe human glioma cell line U251 used in this study was purchased from the American Type Culture Collection (ATCC, Manassas, VA, USA). DNA fingerprinting, cell vitality detection, isozyme detection, and mycoplasma detection were used to characterize all cell lines. The cell line was regularly cultured in Dulbecco's modified Eagle media (DMEM; Gibco; Thermo Fisher Scientific; Waltham, MA, USA) supplemented with 10% fetal bovine serum (FBS; Gibco; Thermo Fisher Scientific) at 37\u0026deg;C in a humid environment with 5% CO\u003csub\u003e2\u003c/sub\u003e.\u003c/p\u003e \u003cp\u003eTransfection of miRNA inhibitors were performed using Hieff TransTM \u003cem\u003ein vitro\u003c/em\u003e siRNA/miRNA Transfection Reagent (Yeasen, China) according to the manufacturer\u0026rsquo;s instructions. U251 cells were plated in T-25 cells culture flasks and 96 wells plates, 24 wells plates, and transfected with miRNA inhibitors. At 48\u0026ndash;72 h after transfection, the cells were collected and used for experiments.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec18\" class=\"Section2\"\u003e \u003ch2\u003eWestern blotting\u003c/h2\u003e \u003cp\u003eThe membrane was incubated with the following primary antibodies (1:1000) at 4\u0026deg;C overnight after being blocked with 5% non-fat milk at room temperature for 1 hour: anti-PD-L1 (Cat#13684), anti-β-actin (Cat#4970), anti-CTLA4 (Cat#53560), and anti-FOXP3 (Cat#19713) antibodies. The membrane was washed three times, followed by treatment with horseradish peroxidase-conjugated anti-mouse or anti-rabbit secondary antibody. After that, a Tanon 5200 system (Tanon, Shanghai, PR China) was used to photograph and wash the membrane. Using Quantity One software from Tanon (Shanghai, PR China), the intensity of protein bands was assessed using densitometric analysis. Target gene protein expression was calibrated to β-actin protein expression. A triplicate of the experiment was run.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec19\" class=\"Section2\"\u003e \u003ch2\u003eCCK-8 assay\u003c/h2\u003e \u003cp\u003eCells were seeded, cultivated for an overnight period, and then given the CCK8 solution in 96-well plates. The cells were then treated with the CCK-8 reagent at 37\u0026deg;C in a humid environment that contained 5% CO\u003csub\u003e2\u003c/sub\u003e. The absorbance at 450 nm was measured on an enzyme analyzer (Tecan infinite M2009PR, Tecan, M\u0026auml;nnedorf, Switzerland). The experiment was carried out three times.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec20\" class=\"Section2\"\u003e \u003ch2\u003eTranswell assay\u003c/h2\u003e \u003cp\u003eCell migration and invasion were performed with Transwell chambers (BD Biosciences, USA). After 48 h of corresponding treatment, U251 cells in each group were inoculated into the upper chamber of the Transwell chamber at the density of 4\u0026times;10\u003csup\u003e4\u003c/sup\u003e cells/well (100 \u0026micro;L culture medium containing 5% fetal bovine serum). Besides, 500 \u0026micro;L culture medium containing 10% fetal bovine serum was added to the 24-well culture plate of the lower chamber. After 24 h of routine culture, the chamber was removed, and cells from the upper layer of microporous membranes were wiped with cotton swabs. After that, cells were immobilized in 4% paraformaldehyde solution for 10 minutes at room temperature and stained with 0.5% crystal violet solution (Sigma-Aldrich; Merck KGaA) for 15 minutes. In the last step, 5 visual fields were randomly selected and observed under an optical microscope (Nikon, Japan) to count the number of cells invading the sub-layer of the microporous membranes of the chamber.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec21\" class=\"Section2\"\u003e \u003ch2\u003eWound healing assay\u003c/h2\u003e \u003cp\u003eCell migration was observed using a wound-healing assay. When the transfected U251 cells were maintained in a 6-well plate, achieving 90\u0026ndash;95% confluence, scratches were generated using the micropipette tips. The wound state was observed at 0 h and 24 h after scratching with an X71 inverted microscope (Olympus, Tokyo, Japan).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec22\" class=\"Section2\"\u003e \u003ch2\u003eStatistical analysis\u003c/h2\u003e \u003cp\u003eUnless stated otherwise, all statistical tests were performed by R software (version:4.1.3) The significance of differences between two groups was determined by Wilcoxon rank sum test. Hypergeometric test was used to determine the significance of overlap between two groups of genes. For patient survival analyses, Kaplan-Meier plots were created, we used the Cox proportional hazard model and the log-rank test to determine the difference of survival between different groups of the patients. Unless stated otherwise, a p-value\u0026thinsp;\u0026lt;\u0026thinsp;0.05 was considered significant. If necessary, p-values were corrected for multiple tests with the Benjamini-Hochberg procedure to calculate false discovery rate (FDR).\u003c/p\u003e \u003c/div\u003e"},{"header":"Results","content":"\u003cdiv id=\"Sec24\" class=\"Section2\"\u003e \u003ch2\u003eMiRNA deregulation mediated by copy number variation and DNA methylation\u003c/h2\u003e \u003cp\u003eGenomic and epigenomic profiling data of copy number variation, DNA methylation, and miRNA expression were obtained for 500 LGG patients from TCGA. We first explored the global CNV characteristics of LGG on the whole genome. Based on the copy number amplitude calculated by GISTIC2, we found that there is a concentrated copy number amplification on q32.1 of chromosome 1, p11.2 of chromosome 7 and 12q.14 of chromosome 11 and copy number deletion on q37.3 of chromosome 2, p21.3 of chromosome 3, q23.31 of chromosome 10, and q13.2 of chromosome 9 (Fig. \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e). The copy number amplitude data was then mapped to miRNA loci to build copy number profiles for 2265 miRNAs. Next, by mapping the 8825 450K array probes to miRNA promoter regions, we obtained DNA methylation profiling data for 1977 miRNAs (Table \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e and Table \u003cspan refid=\"MOESM2\" class=\"InternalRef\"\u003eS2\u003c/span\u003e). The metaplot was used to depict the methylation level around miRNA transcription start sites. In general, the methylation level around transcription start sites was significantly lower than that from gene body and intergenic regions for miRNAs (Fig. \u003cspan refid=\"MOESM2\" class=\"InternalRef\"\u003eS2\u003c/span\u003e). To further investigate the coordination between CNV and DNAm abnormalities, we also compared the frequencies of the miRNA associated aberrant CNV and DNAm in the genome. We define the CNV value\u0026thinsp;\u0026gt;\u0026thinsp;0.3 as CNV amplification (CNVamp); CNV value \u0026lt; -0.3 as CNV deletion (CNVdel), for DNA methylation, the β-value\u0026thinsp;\u0026gt;\u0026thinsp;0.8 was defined as hypermethylation (HyperM); β-value\u0026thinsp;\u0026lt;\u0026thinsp;0.2 was defined as hypomethylation (HypoM). The frequencies of CNVamp, CNVdel, HyperM and HypoM for each sample were calculated. We found that the frequencies of CNVamp are mildly correlated with the frequencies of CNVdel (Fig. S3A, r\u0026thinsp;=\u0026thinsp;0.1, p\u0026thinsp;=\u0026thinsp;0.03), whereas the frequencies of HyperM are significantly negatively correlated with the frequencies of HypoM (Fig. S3B, r = -0.7, p\u0026thinsp;\u0026lt;\u0026thinsp;2.2e\u0026thinsp;\u0026minus;\u0026thinsp;16). In addition, the frequencies of all the aberrant copy number are significantly correlated with the that of all aberrant DNA methylation (Fig. S3C, r\u0026thinsp;=\u0026thinsp;0.11, p\u0026thinsp;=\u0026thinsp;0.015). Abnormalities in directional CNVamp, CNVdel, HyperM and HypoM are also tightly correlated (Fig. S3D-G). In summary, our observation suggests that the LGG patients with frequent aberration of miRNA copy number are more likely to have frequent aberration of promoter methylation. These correlated frequencies of aberrant CNV and DNAm highlights the widespread involvement of genetic and epigenetic regulations for miRNAs in the pathogenesis of LGG.\u003c/p\u003e \u003cp\u003eTo assess the global effects of genomic and/or epigenomic aberrations on miRNA expression, the correlation coefficients between CNV or DNAm with corresponding miRNA expression were calculated and then normalized by using Fisher\u0026rsquo;s Z-transformation. It was found that the overall correlation coefficients between CNV and the corresponding miRNA expression shifted to the right (skewness\u0026thinsp;=\u0026thinsp;2.43, p\u0026thinsp;\u0026lt;\u0026thinsp;2.2e\u0026thinsp;\u0026minus;\u0026thinsp;16, D\u0026rsquo;Agostino test), whereas that between DNAm and the miRNA expression significantly shifted to the left (skewness = -1.11, p\u0026thinsp;\u0026lt;\u0026thinsp;2.2e\u0026thinsp;\u0026minus;\u0026thinsp;16, D\u0026rsquo;Agostino test) (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eA). This indicated the overall positive and negative regulatory effect on miRNA expression by CNV and DNAm, respectively, which is similar to that of protein coding genes. We next focus on the positively correlated miRNAs for CNV (CNV-miR) and the negatively correlated miRNAs for DNAm (DNAm-miR). A total of 351 CNV-miRs (Table S3) and 541 DNAm-miRs (Table S4) were identified at the thresholds of FDR\u0026thinsp;\u0026lt;\u0026thinsp;0.05. Genome-wide landscape for the associations between CNV and DNAm for the two different groups of miRNAs are depicted as Fig. S4. Some top ranked CNV-miRs and DNAm-miRs with significant associations between the expression with CNV and DNAm are presented as Fig. S5 and Fig. S6, respectively. Many of these miRNAs have been identified that play a key role in different cancer types. One of the examples is the hsa-miR-26a, whose expression is highly correlated with CNV value (r\u0026thinsp;=\u0026thinsp;0.62, p\u0026thinsp;\u0026lt;\u0026thinsp;2.2e\u0026thinsp;\u0026minus;\u0026thinsp;16). This miRNA has been found as a key driver gene for glioblastoma, its copy number amplification associated overexpression was observed in glioblastoma to transform cells and to promotes tumor growth (\u003cspan citationid=\"CR63\" class=\"CitationRef\"\u003e63\u003c/span\u003e, \u003cspan citationid=\"CR64\" class=\"CitationRef\"\u003e64\u003c/span\u003e). Another example is miR-155-5p for its negative correlated pattern between expression and promoter methylation (r = -0.74, p\u0026thinsp;\u0026lt;\u0026thinsp;2.2e\u0026thinsp;\u0026minus;\u0026thinsp;16). This negative correlation has been observed in several other studies with regard to glioma, which was indicated that significantly correlated with patients\u0026rsquo; survival (\u003cspan citationid=\"CR65\" class=\"CitationRef\"\u003e65\u003c/span\u003e, \u003cspan citationid=\"CR66\" class=\"CitationRef\"\u003e66\u003c/span\u003e). Similar phenomenon has been observed in many other cancer types including breast cancer (\u003cspan citationid=\"CR67\" class=\"CitationRef\"\u003e67\u003c/span\u003e), gastric cancer (\u003cspan citationid=\"CR68\" class=\"CitationRef\"\u003e68\u003c/span\u003e) and T-cell lymphoma (\u003cspan citationid=\"CR69\" class=\"CitationRef\"\u003e69\u003c/span\u003e), which suggests miR-155-5p act as a pan-cancer regulator with functions regardless of tumor origin.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eWe next inspected the genome distribution of CNV-miRs and DNAm-miRs and found that both groups of miRNAs are mainly distributed on chromosome 14, followed by chromosome 1 and chromosome 19 (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eB). The CNV-miRs and DNAm-miRs present an overlap of 173 miRNAs, which suggests that a significant number of miRNAs are coordinated regulated by CNV and DNAm, whereas others are mutually exclusive regulated (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eC). Additionally, we also checked the positional distribution of the probes within miRNA genes and found that DNAm-miRs probes are more frequently located in promoter regions of 2000 bp upstream from the TSS to 200 bp upstream from the TSS (TSS2000) rather than within 200 bp upstream of the transcription start site (TSS200) and gene body (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eD). Moreover, DNAm-miRs probes are more frequently located in N-shores rather than in S-Shores, OpenSea and CpG island regions (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eE). We further investigated the functions of the miRNA by performing the Gene Ontology enrichment analysis based on the target genes of the two different groups of miRNAs. To do this, we first dissect the influence of DNAm and CNV on miRNA expression, respectively. We use DNAm and CNV as covariates in a regularized linear regression model and subsequently inferring the regulatory connectives. Specifically, we propose expression level of specific miRNA as Y (\u003cem\u003eY\u003c/em\u003e\u003csub\u003e\u003cem\u003emiR\u003c/em\u003e\u003c/sub\u003e), and its methylation level as M (\u003cem\u003eM\u003c/em\u003e\u003csub\u003e\u003cem\u003emiR\u003c/em\u003e\u003c/sub\u003e), as well as the copy number level as N (\u003cem\u003eN\u003c/em\u003e\u003csub\u003e\u003cem\u003emiR\u003c/em\u003e\u003c/sub\u003e), we infer its association coefficients with the DNAm (\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({\\beta }_{m}\\)\u003c/span\u003e\u003c/span\u003e) and CNV (\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\({\\beta }_{n}\\)\u003c/span\u003e\u003c/span\u003e) for each miRNA across all samples as:\u003cdiv id=\"Eque\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Eque\" name=\"EquationSource\"\u003e\n$${Y}_{miR}\\approx {\\beta }_{0}+{\\beta }_{m}*{M}_{miR}+{\\beta }_{n}*{N}_{miR}$$\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003eWe filtered those miRNAs with negative coefficients for DNAm and positive coefficients for CNV under FDR\u0026thinsp;\u0026lt;\u0026thinsp;0.05. In this way, we identified a subset of 425 DNAm-miRs and 286 CNV-miRs that uniquely influenced by these two different mechanisms. Enrichment analysis indicated these DNAm-miRs are mainly involved in nervous system development and brain morphogenesis as well as some tumor related pathways, such as Wnt signaling and Notch signaling (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eF). Whereas for CNV-miRs, functions are mainly about neuron development, cell proliferation and immune response pathways (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eG). These observations indicated that CNV and DNAm mediated miRNA dysregulation may fulfil diverse functions in LGG, and thus it could help to define patients as more appropriate subtypes and to find suitable treatment methods.\u003c/p\u003e \u003cdiv id=\"Sec25\" class=\"Section3\"\u003e \u003ch2\u003eMolecular subtypes defined by CNV-miRs and DNAm-miRs\u003c/h2\u003e \u003cp\u003eNext, we investigated whether the expression profiles of CNV-miRs and DNAm-miRs can be used for patients\u0026rsquo; stratification that associated with prognosis. We applied the non-negative matrix factorization (NMF) clustering analysis to miRNA expression profile with the cluster number k set as 2\u0026ndash;10, the k values were then determined for both CNV-miRs and DNAm-miRs, respectively. The optimal number of k was selected as 4 based on profiles of DNAm-miRs (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eA and Fig. S7), and as 3 based on profiles of CNV-miRs (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eB and Fig. S8). For both clustering results that identified by DNAm-miRs and CNV-miRs, Kaplan\u0026ndash;Meier plot analyses demonstrated that the molecular subtypes among each classification were significant associated with overall survival of the patients (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eC-D). In addition, we compared the molecular subtypes obtained and the histological types and found that these classifications are distinct from each other (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eE), whereas the comparison between the classifications generated from miRNA profiles indicated that high proportions of samples from CNV-miRs based subtypes were remarkably overlapped with the DNAm-miRs based ones (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eF). These results indicate the heterogeneous for LGG samples at molecular level.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eWe next sought to identify molecular subtypes that reflected the multi-layered pattern of the CNV-miRs and DNAm-miRs. We employed an integrated clustering method MOVICS for the combined analysis of the genomic profiling data regarding miRNA expression, DNAm and CNV. We identified four integrative clusters (ICs) from the 10 multi-omics ensemble clustering algorithms. The number of ICs was determined by referring to the cluster prediction index, gap statistical analysis and silhouette score, and was comprehensively evaluated by the consensus clustering matrix (Fig. S9A-C). Then, the clustering results were further combined by the consensus ensemble approach with distinctive molecular patterns of the multi-omics data (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eA). Our classification system was closely related to overall survival of the patients (Log\u0026thinsp;\u0026minus;\u0026thinsp;rank p\u0026thinsp;=\u0026thinsp;3.15e\u0026thinsp;\u0026minus;\u0026thinsp;34; Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eB). Notably, cancer subtype 2 (IC2) exhibited the worst prognosis among all ICs evaluated (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eC). For purpose of comparison, we also examined the overall survival among histological types, and only mild differences for their prognosis can be observed (Log\u0026thinsp;\u0026minus;\u0026thinsp;rank p\u0026thinsp;=\u0026thinsp;0.005; Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eD). There are higher proportions of overlapping samples between histological type of astrocytoma and IC1 cluster, and between oligodendroglioma and IC3 cluster (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003eE). These results indicated that the integrated analysis of multi-omics data associated miRNAs identified clinically-relevant molecular subtypes in which miRNA dysregulation caused by epigenomic and genomic aberrance influences prognostic outcomes.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eTo further identify the clinical characteristics of the four integrative clusters, we compared the tumor grade, 1p/19q co-deletion, MGMT promoter methylation and TERT expression status among the clusters identified. We found that patients in these groups exhibited asymmetric distributions for these characteristics. For example, the major proportion of IC2 samples are grade III tumor, which is significantly higher than that of other groups (p\u0026thinsp;\u0026lt;\u0026thinsp;0.0001; Fig. S10A). In addition, all the IC1 and IC2 cluster samples present 1p/19q co-deletion, in contrast, nearly none of the samples in IC3 present co-deletion, whereas 64% of the samples in IC4 present co-deletion (p\u0026thinsp;\u0026lt;\u0026thinsp;0.0001; Fig. S10B). There are also significant more patients in IC2 present MGMT promoter methylation (p\u0026thinsp;\u0026lt;\u0026thinsp;0.0001; Fig. S10C), and fewer patients in IC2 and IC3 present TERT expression (p\u0026thinsp;\u0026lt;\u0026thinsp;0.0001; Fig. S10D). We next evaluated the prognostic predictive value of the tumor subtypes in different grades. We performed survival analyses for the patients with grade II and III glioma respectively. Stratified survival analyses showed that patients in the IC2 group have a worse prognosis than patients in the other groups in both grade II and III samples (Fig. S10E-F), which indicated the prognostic independence of our glioma classifications method for patient survival.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec26\" class=\"Section3\"\u003e \u003ch2\u003eDistinct genomic features among different molecular subtypes\u003c/h2\u003e \u003cp\u003eFinally, we examined the single nucleotide variation (SNV) in different subtypes obtained to inspect the molecular heterogeneity of LGG samples. We screened a set of genes that altered significantly among subtypes, surprisingly, the somatic mutational spectrum of these genes agrees fantastically with the clusters identified. As shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eA, we found that \u003cem\u003eIDH1\u003c/em\u003e mutates predominantly in clusters of IC1, IC3 and IC4, and the mutation frequency of this gene was further showed statistically significant differences among clusters (chi-square test p\u0026thinsp;=\u0026thinsp;3.53e\u0026thinsp;\u0026minus;\u0026thinsp;64, Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eB, left panel), whereas the \u003cem\u003eTP53\u003c/em\u003e predominates in IC1 (chi-square test p\u0026thinsp;=\u0026thinsp;1.85e\u0026thinsp;\u0026minus;\u0026thinsp;66, Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eB, middle panel). In contract, IC2 cluster bears nearly all the mutations for \u003cem\u003eEGFR\u003c/em\u003e (chi-square test p\u0026thinsp;=\u0026thinsp;3.33e\u0026thinsp;\u0026minus;\u0026thinsp;29, Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eB, right panel). A mutually exclusive mutation pattern for the \u003cem\u003eIDH1\u003c/em\u003e and \u003cem\u003eEGFR\u003c/em\u003e can be observed. In addition, some distinctly different mutant genes including \u003cem\u003eATRX\u003c/em\u003e, \u003cem\u003eCIC\u003c/em\u003e and \u003cem\u003eFUBP1\u003c/em\u003e were also identified in our analysis. We next checked the associations between mutational spectrums of these genes and the patients\u0026rsquo; survival. In general, only mild or no survival differences can be observed between groups defined by the mutational status of these genes (Fig. S11). This analysis indicates that our suggested classification system is superior to that based on mutations of individual genes for prognosis of LGG patients. We further extended our analysis to all mutated genes and calculated tumor mutation burden (TMB) for all the samples, as a result, IC2 presents highest mutation load across all clusters (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eC), which indicated genetic heterogeneity are highly correlated with the subtypes defined by DNA methylation and copy number associated miRNAs.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eInspired by the correlations between the mutational spectrum and the multi-omics clustering results for LGG, we further checked the pathways affected in order to obtain the potential targets for precise treatment, particularly for IC2. By using maftools, we identified the most significant pathways including \u0026ldquo;RTK-RAS\u0026rdquo;, \u0026ldquo;WNT\u0026rdquo;, \u0026ldquo;NOTCH\u0026rdquo; etc., which are differentially affected between IC2 and other groups for the genes involved (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eD). The potential druggable gene categories for IC2 and other groups are shown in Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003eE-F. We could see that the top list genes for IC2 includes the \u003cem\u003eCOL6A3\u003c/em\u003e, \u003cem\u003eEGFR\u003c/em\u003e, \u003cem\u003eMUC16\u003c/em\u003e, \u003cem\u003eOBSCN\u003c/em\u003e and \u003cem\u003ePIK3CA\u003c/em\u003e, whereas \u003cem\u003eARID1A\u003c/em\u003e, \u003cem\u003eATRX\u003c/em\u003e, \u003cem\u003eBCOR\u003c/em\u003e, \u003cem\u003eCIC\u003c/em\u003e and \u003cem\u003eFUBP1\u003c/em\u003e are included for other groups. These intriguing results strongly suggest that the distinct mutational events might synergistically involve in miRNA based LGG subtyping and thus provides novel insights to development of potential therapeutic strategies for clinical utility in the future.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec27\" class=\"Section3\"\u003e \u003ch2\u003eImmune landscape of LGG molecular subtypes\u003c/h2\u003e \u003cp\u003eThe composition of the tumor immune microenvironment (TIME) has been identified as the critical factor of tumor\u0026ndash;immune interactions to influence clinical outcome of patients and also the response to immune treatment (\u003cspan citationid=\"CR70\" class=\"CitationRef\"\u003e70\u003c/span\u003e). Many studies have been conducted on the TIME to uncover the underlying mechanisms of cancer (\u003cspan citationid=\"CR71\" class=\"CitationRef\"\u003e71\u003c/span\u003e). To evaluate whether the difference in multi-omics pattern imply immune heterogeneity and microenvironment of LGG, we thus investigated immune cell infiltrations and immune functions discrepancy among different clusters. We first focus on the 17 immunologically relevant pathways derived from ImmPort (\u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e). We calculated the gene expression fold changes between IC2 and other clusters and performed the GSEA for the collected immune pathways. Analysis results showed that more than half of the pathways are enriched by the genes present dysregulation between glioma clusters, such as interleukins, chemokines, cytokines and cytokine receptors (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eA). Next, we calculated the immune score by ESTIMATE for each tumor sample. We observe the immune scores were significantly increased in IC2 cluster, which suggests patients in IC2 cluster present higher infiltration of immune cells in tumor sample (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eB). In addition, IC2 cluster also yields higher stromal scores than other groups, but lower tumor purity scores obtained by ESTIMATE (Wilcoxon p\u0026thinsp;\u0026lt;\u0026thinsp;0.05, Fig. S12), which indicate that IC2 cluster samples are more heterogeneous than that of other groups. Moreover, we also calculated the antitumor immunoactivity for glioma patients, including CYT, MHC and CTL scores. All the three immunoactivity scores were significantly increased in IC2 cluster (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eC, left panel), indicating patients in IC2 cluster displayed greater tumor cell cytotoxic activity and stronger antigen recognition capacities. In addition, we collected eight different immune signature gene sets and the ssGSEA was uses to estimate the enrichment score of tumor samples (Table S5), also, IC2 cluster also presents higher enrichment of different immune signatures (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eC, right panel). For instance, higher IFN-γ score were observed in IC2 cluster, which indicated that IC2 may have stronger anti-tumor immunity and inhibited angiogenesis activities in tumor, thus to influence the efficiency of cancer immunotherapies (\u003cspan citationid=\"CR72\" class=\"CitationRef\"\u003e72\u003c/span\u003e).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eWe next evaluated the relative proportion of immune cells of glioma patients using TIMER. Similarly, we found the infiltrate levels of CD4\u0026thinsp;+\u0026thinsp;T, neutrophil, macrophage and dendritic cells were highest in IC2 among all cluster (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eD), which presents the typical phenotypes for immune \u0026ldquo;hot\u0026rdquo; tumor. Finally, we collected 63 well-known marker genes of 21 immune cells from public resources (Table S6), the clustering heatmap shows that these genes exhibited diverse expression patterns among clusters, where samples in IC2 showed higher expression in general. (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003eE). Collectively, these results suggested the DNAm and CNV mediated the heterogeneity of miRNA expression is highly correlated with tumor immune status and could determine the immunophenotypes of glioma.\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv id=\"Sec28\" class=\"Section2\"\u003e \u003ch2\u003eImmune phenotype associated miRNA identification and prognosis model construction\u003c/h2\u003e \u003cp\u003eHaving demonstrated the relationship between immune phenotype heterogeneity and the DNAm as well as CNV associated miRNA expression dysregulation, we next sought to comprehensively identify the miRNAs that determine the immune cell infiltrations and immune functions. By calculating the Pearson correlation coefficient between cell fractions and expression of DNAm-miRs and CNV-miRs, we identified 517 and 325 miRNAs significantly correlated with infiltrating immune cells respectively. In general, the number of positively correlated miRNAs was higher than that of negatively correlated ones for both groups (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eA). Interestingly, we observed that the number of miRNAs with immune cell infiltration associations increases with that of cell types, and there are 195 miRNAs were identified that correlated with all the six immune cell infiltrations (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eB, Table S7). In addition, we observe that miRNA expression correlated with diverse immune cell type infiltration in majority of the clusters identified. Most miRNAs were positively associated with immune cells, whereas the negative correlations are mainly occurring in CD4\u0026thinsp;+\u0026thinsp;T cells (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eC). In order to further explore effects of the immune-related key miRNAs on LGG, we collected experimentally verified immune related miRNAs from the RNA2Immune database, which is a public repository to catalogue ncRNA-immune associations (\u003cspan citationid=\"CR56\" class=\"CitationRef\"\u003e56\u003c/span\u003e). Here we selected 464 and 574 miRNAs from the sections of \u0026ldquo;immune cell function\u0026rdquo; and \u0026ldquo;cancer immunity\u0026rdquo;, respectively (Table S8). By comparing with those miRNAs that correlated with immune cell infiltrations, a total of 64 miRNAs were identified that overlap with the external validations (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eD), in which 47 miRNAs are from DNAm-miRs and 17 are from CNV-miRs (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eE).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eStudies have extensively indicated that miRNAs present strong potential for predicting cancer prognosis and exhibited great potential for clinical applications as therapeutic targets (\u003cspan citationid=\"CR73\" class=\"CitationRef\"\u003e73\u003c/span\u003e). Thus, we integrated the survival information to evaluate the potential ability of the identified miRNA as candidate prognostic marker for LGG. We performed a Cox proportional hazard regression analyses to identify the survival related miRNAs. Under a stringent threshold of p\u0026thinsp;\u0026lt;\u0026thinsp;0.0005, a total of 34 miRNAs were identified that correlated with patients\u0026rsquo; survival. Interestingly, we found that almost all of these miRNAs are positively correlated with immune cell infiltration levels, and only miR-9-5p exhibited negative correlations (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eF). To construct the prognostic prediction model, a LASSO proportional hazards regression analysis was then performed to further determine the relationship between patients\u0026rsquo; survival and the expressions of the identified 34 miRNAs. We found a significant correlation with the overall survival of LGG patients with the following 8 miRNAs: (hsa-miR-10b-5p, hsa-miR-155-5p, hsa-miR-196a-5p, hsa-miR-196b-5p, hsa-miR-200a-3p, hsa-miR-204-5p, hsa-miR-503-5p, hsa-miR-15b-5p). These markers were considered to function as prognostically significant immune-related miRNAs. We then implemented these 8 miRNAs in development of a miRNA based prognostic risk score that determine the degree of tumor immune infiltration. These 8 miRNAs were weighted based on the LASSO regression coefficient as follows:\u003c/p\u003e \u003cp\u003eRisk Score = (0.1617*\u003cem\u003eExp\u003c/em\u003e\u003csub\u003emiR\u0026minus;10b\u0026minus;5p\u003c/sub\u003e) + (0.3419*\u003cem\u003eExp\u003c/em\u003e\u003csub\u003emiR\u0026minus;155\u0026minus;5p\u003c/sub\u003e) + (0.0695*\u003cem\u003eExp\u003c/em\u003e\u003csub\u003emiR\u0026minus;196a\u0026minus;5p\u003c/sub\u003e) + (0.1391*\u003cem\u003eExp\u003c/em\u003e\u003csub\u003emiR\u0026minus;196b\u0026minus;5p\u003c/sub\u003e) + (0.1474*\u003cem\u003eExp\u003c/em\u003e\u003csub\u003emiR\u0026minus;200a\u0026minus;3p\u003c/sub\u003e) + (0.0840*\u003cem\u003eExp\u003c/em\u003e\u003csub\u003emiR\u0026minus;204v5p\u003c/sub\u003e) + (0.0996*\u003cem\u003eExp\u003c/em\u003e\u003csub\u003emiR\u0026minus;503\u0026minus;5p\u003c/sub\u003e) + (0.2083*\u003cem\u003eExp\u003c/em\u003e\u003csub\u003emiR\u0026minus;15b\u0026minus;5p\u003c/sub\u003e)\u003c/p\u003e \u003cp\u003eWe observed that the IC2 cluster samples present height risk scores among all groups as expected (Wilcoxon p\u0026thinsp;\u0026lt;\u0026thinsp;2.2e\u0026thinsp;\u0026minus;\u0026thinsp;16, Fig. S13). Besides, all the samples in this training set were separated into high-risk (N\u0026thinsp;=\u0026thinsp;84) and low-risk groups (N\u0026thinsp;=\u0026thinsp;411) according to the level of the risk scores, the distribution of risk score among high-risk and low-risk groups and expression heatmap of the miRNAs are indicated in Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eG. We can observe a clear distinction between for the signature miRNA expression from the two groups of patients.\u003c/p\u003e \u003cp\u003eIn order to determine whether the prognostic model was robust, a validation set obtained from CGGA was utilized to verify the validity and accuracy of the prognostic model based on the 8 miRNAs. We collected miRNA microarray data for 198 samples from this cohort, 188 of which have the survival data available, each sample of the validation set also acquired a risk score according to the same formula, and was then stratified as being either high-risk (N\u0026thinsp;=\u0026thinsp;121) or low-risk groups (N\u0026thinsp;=\u0026thinsp;67). The risk score distribution and gene expression data for the validation set are shown in Fig. S14. Patients with higher risk scores were noted to have worse survival in contrast to those with lower risk scores, the overall survival analysis revealed significant difference between the two groups (Log\u0026thinsp;\u0026minus;\u0026thinsp;rank p\u0026thinsp;=\u0026thinsp;4.24e\u0026thinsp;\u0026minus;\u0026thinsp;08, Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eH). The predictive potential of the prognostic model using time-dependent ROC curves was shows as Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eI. The area under the ROC curve (AUC) of the prognostic model for OS was 0.65 at 1 year, 0.75 at 2 years, 0.78 at 3 years and 0.76 at 5 years, suggesting that our prediction model has favorable efficacy for predicting both short- and long-term prognosis.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec29\" class=\"Section2\"\u003e \u003ch2\u003eqRT-PCR validation of signature miRNAs and clinical utility of risk model\u003c/h2\u003e \u003cp\u003eWe next validated the miRNAs as novel prognosis markers for glioma in our in-house cohort. We detected the expression level of the eight miRNAs form tumor and normal samples of 30 LGG patients by qRT-PCR (Table S9). The normalized expression levels from qRT-PCR showed that all these miRNAs present elevated expression in LGG tumor sample, although with different significance obtained. Two miRNAs of miR-10b-5p and miR-204-5p present mild upregulation for about 1.5 times of fold change, whereas others are significantly up-regulated in about 20\u0026ndash;50 times, Notably, miR-155-5p showed highest level of fold change in 56.4 times (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eA). These cross platform and cross-racial data indicated the robustness of these signature miRNAs presented here for patient\u0026rsquo;s prognosis.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eFurthermore, we conducted univariate and multivariate Cox regression analyses to explore whether the prognostic value of our risk model was independent of other clinical features in the TCGA LGG cohort. For the univariate Cox regression model, our prediction classifier was the strongest variable correlated with prognosis among clinical features including histological type, IDH mutation, 1p/19q codeletion, MGMT promoter methylation and TERT expression status. After multivariate adjustment by the clinical features, our classifier remains the strongest and independent prognostic factor in survival analysis (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eB). This results indicated that our model performed better than conventional clinicopathological features for clinical outcome prediction. To provide clinically correlated quantitative approach to predicting the prognosis of LGL patients, a nomogram that integrated the risk model and independent clinical risk factors (tumor grade and 1p/19q codeletion) was constructed (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eC). Our risk model was found to contribute the most risk points (ranging from 0 to 100) compared with the other clinical features, which was consistent with the Cox multivariate regression results. The calibration plot indicated that the nomogram performed well compared with an ideal model (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eD). The AUC score at 1-, 2-, 3- and 5- years were 0.86, 0.84, 0.83 and 0.80 for the nomogram, respectively (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eE). In summary, these findings suggest that the nomogram was a better model for predicting short-term or long-term survival of LGG patients than individual prognostic factors.\u003c/p\u003e \u003c/div\u003e\n\u003ch3\u003emiRNAs associated with immune evasion of LGG\u003c/h3\u003e\n\u003cp\u003eSince the IC2 cluster exhibited higher levels of immune cell infiltration, whereas patients in this cluster had poorer prognosis than others, we conducted in-depth functional investigations to explore possible reasons for this contradiction. We first performed functional enrichment analysis of GSVA on Hallmark gene sets obtained from MSigDB. Specially, besides the higher immune activities, the IC2 cluster is also endowed with higher energy expenditure states, such as glycolysis, hypoxia and angiogenesis, whereas IC4 is more closely related to metabolic activities, such as the bile acid and fatty acid metabolism, as well as protein secretion. In contract, IC3 and IC1 clusters are more distinguished by proliferative-related pathways, such as mitotic spindle, E2F targets, and G2M checkpoints, and also by classical tumor pathways including the MYC and Wnt-signaling (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003eA). Next, we further investigated whether the clusters obtained are correlated with immunosuppressive effects. We used the methodology from Shuichi \u003cem\u003eet al\u003c/em\u003e. to calculate three different scores in order to measure tumor proliferation and glycolysis rates (T-score), immune cell and antitumor immunity level (I-score), as well as immunosuppressive level (S-score), respectively. Here, the IC2 cluster conveys the highest for all three different scores, particularly for the S-score, which is much higher in IC2 than other groups (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003eB), which indicated an immune evasion state may occurred in the patients of IC2 cluster. We next inspected \u003cem\u003ePD-L1\u003c/em\u003e expression and found its expression is much higher in IC2 than other groups (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003eC), which verify the hypothesis of an immune evasion status for the IC2 cluster.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eTo identify key miRNAs that involved in LGG immune evasion, we collected experimentally validated regulatory associations between \u003cem\u003ePD-L1\u003c/em\u003e and miRNAs from miRTarBase, and further screened those by performing a Pearson correlation analysis between the expression of \u003cem\u003ePD-L1\u003c/em\u003e and DNAm-miRs as well as CNV-miRs. Base on negative correlation of PCC\u0026thinsp;\u0026lt;\u0026thinsp;0 and the threshold of FDR\u0026thinsp;\u0026lt;\u0026thinsp;0.05, we discovered many miRNA modulators that play key roles in \u003cem\u003ePD-L1\u003c/em\u003e mediated immune evasion. One of the examples is the miR-338-5p, our analysis indicated this miRNA presents significant negative correlation with \u003cem\u003ePD-L1\u003c/em\u003e for their expressions (r\u0026thinsp;=\u0026thinsp;\u0026minus;\u0026thinsp;0.1, p\u0026thinsp;=\u0026thinsp;0.0309) (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003eD), this is possibly due to the DNA methylation dysregulation in its promoter region, a negative correlation between expression of this miRNA and the promoter methylation was observed (r\u0026thinsp;=\u0026thinsp;\u0026minus;\u0026thinsp;0.15, p\u0026thinsp;=\u0026thinsp;0.0011) (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003eE). This regulatory relationship between miR-338-5p and \u003cem\u003ePD-L1\u003c/em\u003e has been identified for regulatory T cells (Tregs) mediated immune evasion (\u003cspan citationid=\"CR74\" class=\"CitationRef\"\u003e74\u003c/span\u003e). Another example is the miR-17-5p, whose expression is significantly negatively correlated with that of \u003cem\u003ePD-L1\u003c/em\u003e (r\u0026thinsp;=\u0026thinsp;\u0026minus;\u0026thinsp;0.12, p\u0026thinsp;=\u0026thinsp;0.0056) (Fig. S15A), and its expression dysregulation was identified that caused by its copy number variation (r\u0026thinsp;=\u0026thinsp;0.24, p\u0026thinsp;=\u0026thinsp;7.14e\u0026thinsp;\u0026minus;\u0026thinsp;08) (Fig. S15B). Inverse correlation between \u003cem\u003ePD-L1\u003c/em\u003e and miR-17-5p has been confirmed in many cancer types. For instance, up-regulation of \u003cem\u003ePD-L1\u003c/em\u003e due to post-transcriptional control by miR-17-5p confers resistance to BRAF or MEK inhibitors (BRAFi or MEKi) for metastatic melanoma (\u003cspan citationid=\"CR75\" class=\"CitationRef\"\u003e75\u003c/span\u003e). In addition, ceRNA regulatory loops involving miR-17-5p/PD-L1 have been observed in lung cancer and breast cancer to promote tumor progression and immune checkpoint blockades (ICB) resistance (\u003cspan citationid=\"CR76\" class=\"CitationRef\"\u003e76\u003c/span\u003e, \u003cspan citationid=\"CR77\" class=\"CitationRef\"\u003e77\u003c/span\u003e). Besides these, similar regulatory interactions were also identified for eight miRNAs including let-7a-3p, miR-30b-3p, miR-9-3p, miR-548f-3p, miR-139-3p miR-9-5p, miR-128-3p and miR-4797-3p. These miRNAs present significant negative correlations for their expression with \u003cem\u003ePD-L1\u003c/em\u003e (Fig. S16), which indicated this miRNA regulatory module may play critical roles in \u003cem\u003ePD-L1\u003c/em\u003e associated immune evasion for LGG (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003eF).\u003c/p\u003e \u003cp\u003eTo provide a more comprehensive view of miRNA related tumor immune evasion mechanism, we used TIDE method to evaluate the efficacy of ICB based immunotherapy for LGG. We found patients in the IC2 and IC3 clusters presents significant higher TIDE scores than other groups, which results in a lower efficacy of immunotherapy (Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003eA). According to TIDE analysis, only about one quarter of the patients in IC2 and IC3 may response to \u003cem\u003ePD-1\u003c/em\u003e and \u003cem\u003eCTLA-4\u003c/em\u003e based treatment, in contract, about 64%-79% of the patients in IC1 and IC4 could benefit from the treatment (Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003eB). By Pearson correlation analysis between TIDE scores and miRNA expression, we identified 225 miRNAs with significant correlation under FDR\u0026thinsp;\u0026lt;\u0026thinsp;0.05, which may involve in immunotherapy response. One such example is miR-10b-5p, a significant positive correlation was observed between TIDE score and miR-10b-5p expression (Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003eC). This miRNA has been identified that impairs TET2-mediated PD-L1 transcription inhibition, by which to promote immune evasion and tumor progression in glioblastoma (\u003cspan citationid=\"CR78\" class=\"CitationRef\"\u003e78\u003c/span\u003e). By further integration of the TIDE score correlated miRNAs and the putative targets obtained from public resources, we constructed a miRNA immune evasive network to describe the mechanism explaining roles of miRNA involved in immune envision regulation (Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003eD). Within this network, 29 unique targets of 82 immune envision related miRNAs that under the CNV and DNAm burden in LGG were identified. We screened genes that simultaneously targeted by two or more miRNAs and found \u003cem\u003ePLAGL2\u003c/em\u003e, \u003cem\u003eSCD1\u003c/em\u003e, \u003cem\u003ePIK3CB\u003c/em\u003e, \u003cem\u003eELOVL6\u003c/em\u003e, \u003cem\u003eNR2E1\u003c/em\u003e and \u003cem\u003eMET\u003c/em\u003e demonstrated the highest connectivity, which suggests the importance of these genes in immune envision and the associated immunotherapy response for LGG. Indeed, many of these genes have been validated for their functions in immune regulation and cancer development. For instance, the \u003cem\u003ePLAGL2\u003c/em\u003e and the \u003cem\u003ePIK3CB\u003c/em\u003e were identified as key hub gene for the regulatory network that significantly associated with tumor immune infiltration and sensitivity of chemotherapeutic drugs for cancer (\u003cspan citationid=\"CR79\" class=\"CitationRef\"\u003e79\u003c/span\u003e, \u003cspan citationid=\"CR80\" class=\"CitationRef\"\u003e80\u003c/span\u003e). \u003cem\u003eSCD1\u003c/em\u003e was found that expressed in cancer cells and immune cells causes immune-resistant conditions, inhibition of \u003cem\u003eSCD1\u003c/em\u003e enhances the recruitment of dendritic cells into tumors and the subsequent induction and tumor accumulation of antitumor CD8\u0026thinsp;+\u0026thinsp;T cells, thus this gene was identified as a potential target to enhance the antitumor effects of an anti-PD-1 therapy (\u003cspan citationid=\"CR81\" class=\"CitationRef\"\u003e81\u003c/span\u003e). For miRNAs, miR-214-3p, miR-25-3p and miR-27b-3p were identified as key hub regulators with highest connectivity, and also have been validated for their roles in immune escape (\u003cspan citationid=\"CR82\" class=\"CitationRef\"\u003e82\u003c/span\u003e, \u003cspan citationid=\"CR83\" class=\"CitationRef\"\u003e83\u003c/span\u003e). In summary, our analysis revealed miRNAs act as key regulators for genes that involved in cancer immune evasion, and thus could act as an adjuvant to several major immunotherapies by synergistically increasing their efficacy potential.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003cb\u003eValidation of the association between prognosis related miRNAs and immune evasion\u003c/b\u003e \u003cb\u003ein vitro\u003c/b\u003e \u003cb\u003eexperiments\u003c/b\u003e\u003c/p\u003e \u003cp\u003eWe next evaluated the capacities of prognosis related miRNAs as novel immune evasion regulators for glioma \u003cem\u003ein vitro\u003c/em\u003e experiments. miRNA inhibitor technology was used to knock down the miRNAs in the U251 cell line and measure the expression of miRNAs. We found that miRNA inhibitor significantly reduced the expression of miRNAs compared with transfected NC inhibitor cells (Fig.\u0026nbsp;\u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e10\u003c/span\u003eA). We further performed CCK8 to investigate the effect of miRNAs on the viability of glioma cells. Silencing of miR-155-5p, miR-196a-5p, miR-196b-5p, miR-200a-3p, miR-503-5p, and miR-15b-5p delayed the proliferation of U251 cell lines, while silencing of miR-10b-5p or miR-204-5p has no effect (Fig.\u0026nbsp;\u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e10\u003c/span\u003eB). Wound healing and transwell assays were employed to explore whether cell migration and invasion were influenced after silencing of miRNAs. We found that the knockdown of miR-155-5p, miR-196a-5p, miR-196b-5p, miR-200a-3p, miR-503-5p, and miR-15b-5p significantly decreased cell migration and invasion compared with the control, while silencing of miR-10b-5p or miR-204-5p has no effect (Fig.\u0026nbsp;\u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e10\u003c/span\u003eC-F). Moreover, the knockdown of miR-155-5p, miR-196a-5p miR-196b-5p, miR-200a-3p, miR-503-5p, and miR-15b-5p significantly inhibited the well-known immune evasion markers of \u003cem\u003ePD-L1\u003c/em\u003e, \u003cem\u003eCTLA4\u003c/em\u003e, and \u003cem\u003eFOXP3\u003c/em\u003e at the protein level, while silencing of miR-10b-5p or miR-204-5p has no effect (Fig.\u0026nbsp;\u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e10\u003c/span\u003eG). Taken together, these results suggested that miR-155-5p, miR-196a-5p, miR-196b-5p, miR-200a-3p, miR-503-5p, and miR-15b-5p significantly promote \u003cem\u003ein vitro\u003c/em\u003e cell migration, invasion, proliferation, and are associated with immune evasion in glioma.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e"},{"header":"Discussion","content":"\u003cp\u003eMiRNAs expression deregulation has been widely observed in cancer development and progression. However, less attention has been paid towards understanding the underlying mechanisms of miRNA expression dysregulation. In recent years, more and more evidences have indicated the genetic and epigenetic machineries are extensively involved in the process of miRNA biogenesis and impair its expression in cancer cells (\u003cspan citationid=\"CR84\" class=\"CitationRef\"\u003e84\u003c/span\u003e, \u003cspan citationid=\"CR85\" class=\"CitationRef\"\u003e85\u003c/span\u003e). Here we focus on two widely-known genetic and epigenetic processes of DNA methylation and copy number variation, in order to get deeper insights into this intricate regulation as well as their roles in cancer development. We developed strategies to obtain promoter methylation and copy number variation for miRNA genes from high-throughput profiling, combined with miRNA expression, we obtained miRNAs whose expression affected by DNAm and CNV dysregulation in lower grade glioma, and also identified prognostic subtypes based on multidimensional omics analysis. Different with previous studies which only investigate individual mechanism (\u003cspan citationid=\"CR86\" class=\"CitationRef\"\u003e86\u003c/span\u003e, \u003cspan citationid=\"CR87\" class=\"CitationRef\"\u003e87\u003c/span\u003e), we sought to inspect and give a more comprehensive description for both effect of DNAm and CNV on miRNA expression simultaneously. Our analysis not only provides novel molecular subtyping system for LGG, but also reveals that the regulatory mechanism surrounding miRNA may be closely related to tumor microenvironment and immune evasion, as well as the resulting prognostic outcome differences. Furthermore, this tumor classification around miRNAs is also closely related to genome wide single nucleotide polymorphisms. This work indicated that integrative analysis of multi-omics profiling focus on miRNA regulation could identify meaningful molecular subtypes which are in connection with mechanisms and sample heterogeneity in tumors as well as potential biomarkers and therapy targets.\u003c/p\u003e \u003cp\u003eThe relationship between genetic and epigenetic regulatory mechanisms has been extensively investigated by using integrative multi-omics study for molecular subtyping (\u003cspan citationid=\"CR88\" class=\"CitationRef\"\u003e88\u003c/span\u003e, \u003cspan citationid=\"CR89\" class=\"CitationRef\"\u003e89\u003c/span\u003e). Nonetheless, illustrating the complex interactions for DNAm- or CNV- related miRNA dysregulation remains limited, mainly due to that obtaining high-throughput profiles around miRNA at genetic and epigenetic layers is still challenging. Toward this end, we developed novel strategy and successfully extracted miRNA promoter methylation information by using credible miRNA promoter annotation and repurposing the 450K microarray data, which were originally designed to detect the methylation of protein coding genes. Although many miRNA genes reside within intronic regions and are thought to co-transcribed with their host genes, increasing evidences have indicated that intronic miRNAs carry their own promoters, therefore have independent transcription (\u003cspan citationid=\"CR90\" class=\"CitationRef\"\u003e90\u003c/span\u003e). In this case, we decided to use experimentally verified promoter region for individual miRNAs but without that from host genes being used. In addition, we also obtained miRNA CNV profiling by examining miRNAs \u003cem\u003eloci\u003c/em\u003e which localize in somatic copy number variation regions. Our results confirmed such associations between DNAm or CNV aberration and miRNA expression. Akin to protein coding genes, the methylation at promoter regions and the expression for corresponding miRNAs are often observed to have negative correlation, whereas for CNVs, the positive correlations dominant for all associations, although few opposite trends also can be observed. We found that about one third of the miRNAs in the human genome tend to be affected by either copy number or methylation alterations, and many of them were driven by a combination of both mechanisms. This phenomenon was also observed for protein coding genes from glioblastoma and ovarian cancer (\u003cspan citationid=\"CR91\" class=\"CitationRef\"\u003e91\u003c/span\u003e), suggesting the crosstalk between genetic and epigenetic layers are widespread in cancer cell, this synergistic regulation between different regulatory mechanisms of DNAm and CNV may function in maintaining the precise expression of genes.\u003c/p\u003e \u003cp\u003eThe concomitant regulation for miRNA expression by alterations in genomic CNV and DNAm was also indicted by the genomic distribution of miRNA genes affected. Surprisingly, both the CNV-miRs and DNA-miRs exhibited genomic preference for three different chromosomes, especially the chromosome 14. In addition, the frequencies of the aberrant DNA copy number and DNA methylation are also correlated. The concentration distribution of miRNA genes may raise from expansion of miRNA gene families in the genome. Several studies have indicated that some miRNA families may evolved from whole genome or tandem duplication events, their family members underwent fast amplification and rapid divergence, thus may have conserved regulatory elements (\u003cspan citationid=\"CR92\" class=\"CitationRef\"\u003e92\u003c/span\u003e). This could be another mechanism for elaborating the regulation of miRNA family members. This co-regulation in the genome may allow them to target the genes within key pathways that involved in cancer, and also provides us another perspective to understand the mechanism of miRNA evolution.\u003c/p\u003e \u003cp\u003eOur classification analysis based on CNV and DNAm associated miRNAs revealed novel molecular key features that represent novel precision therapy targets and biomarkers for LGG management. These LGG samples were subdivided into four subtypes which shows significant prognostic difference. Survival analysis showed that patients\u0026rsquo; survival was much lower in the IC2 subtype than that in other three subtypes. Results indicated that our classification is superior to those of protein-coding gene based mutational profiling or multi-omics classification methods for prediction of patients\u0026rsquo; prognosis (\u003cspan citationid=\"CR93\" class=\"CitationRef\"\u003e93\u003c/span\u003e, \u003cspan citationid=\"CR94\" class=\"CitationRef\"\u003e94\u003c/span\u003e). As a kind of highly effective gene expression regulator, one miRNA may regulate dozens to hundreds of target genes. One could hypothesize that disturbance to only a few miRNAs may cause huge downstream runaway effect (\u003cspan citationid=\"CR95\" class=\"CitationRef\"\u003e95\u003c/span\u003e). In such case, miRNAs show a higher value for tumor classification, development of relevant tools around miRNA omics data may provide us a better alternative for early diagnosis and accurate prognosis predictions in cancers.\u003c/p\u003e \u003cp\u003eAnalysis of different biological characteristics among subtypes revealed intra-tumor heterogeneity and molecular interpretability of LGG. Significant differences in the tumor immune microenvironment of the molecular subtypes were observed. Studies have extensively shown that tumor infiltrating lymphocytes (TILs) are involved in tumor progression and invasiveness. TILs include various lymphocytes with different activities, particularly the CD8\u0026thinsp;+\u0026thinsp;and CD4\u0026thinsp;+\u0026thinsp;T cells, which are usually associate with favorable prognostic in patients (\u003cspan citationid=\"CR96\" class=\"CitationRef\"\u003e96\u003c/span\u003e). As found in our study, higher level of TILs but the poor survival of IC2 indicated this group as an immunosuppressive subtype. miRNAs are emerging as important regulators involved in the cancer immune process and play essential roles in cancer immunotherapy (\u003cspan citationid=\"CR97\" class=\"CitationRef\"\u003e97\u003c/span\u003e). Therefore, it is necessary to determine the relationship between miRNAs and immune regulation in LGG. By estimating the correlations with infiltrating immunocytes and immune pathway activities, we identified immune-related miRNAs that may have prognosis value for LGG patients with 8 essential immune miRNAs were prioritized, many of them have been reported to be associated with immune regulation and cancer therapy from previous studies. For instance, miR-15b-5p mediates the degradation of \u003cem\u003ePD-L1\u003c/em\u003e mRNA through interaction with its 3'-untranslated region, by which to induce CD8\u0026thinsp;+\u0026thinsp;T and NK cell activation and cytotoxicity against neuroblastoma (\u003cspan citationid=\"CR98\" class=\"CitationRef\"\u003e98\u003c/span\u003e). Similar observations also apply to another miRNA of miR-155-5p, by targeting \u003cem\u003ePD-L1\u003c/em\u003e, this miRNA was found to activated CD8\u0026thinsp;+\u0026thinsp;T cell function and suppress ovarian cancer progression (\u003cspan citationid=\"CR99\" class=\"CitationRef\"\u003e99\u003c/span\u003e). The prognosis model constructed based on the expression of these miRNAs was validated by an independent cohort. It is worth noting that despite the different platforms and technologies being used to assess expression status, these miRNAs were consistently found that associated with patient\u0026rsquo;s risk in the discovery and validation cohort, suggesting the robust of the proposed prediction model. These results indicated that the integration of DNAm and CNV is an effective method to study the immune regulation of miRNAs for cancer prognosis.\u003c/p\u003e \u003cp\u003eRemarkably, our proposed sample classification has a good concordance with the mutation profiles among the subtypes, particularly for \u003cem\u003eIDH1\u003c/em\u003e, \u003cem\u003eTP53\u003c/em\u003e and \u003cem\u003eEGFR\u003c/em\u003e. Somatic mutation is another mechanism that may impact the biogenesis of miRNAs. However, the detailed process remains unrevealed. One inspiring example is the \u003cem\u003eIDH1\u003c/em\u003e, somatic mutations of this gene are strongly associated with the DNA methylome and transcriptome reorganization in LGG (\u003cspan citationid=\"CR100\" class=\"CitationRef\"\u003e100\u003c/span\u003e). Whether DNA methylation mediates this indirect association between \u003cem\u003eIDH1\u003c/em\u003e mutation and miRNA expression remains further analysis. Besides, the association of miRNA expression and mutation frequencies of other genes within LGG remains unknown. Nonetheless, our present findings provide potential therapy targets for corresponding subtypes and brings novel enlightenment to the study of the intrinsic regulatory mechanisms for LGG.\u003c/p\u003e"},{"header":"Conclusions","content":"\u003cp\u003eIn summary, we have designed and implemented a novel strategy to facilitate the integration of copy number variation, DNA methylation and miRNA expression data. Analysis results demonstrate the added value of for molecular subtyping of LGG by applying our method to the multi-omics data. Systematic study also characterized the immune landscape for glioma patients based on miRNAs related multi-omics data and depicted the crosstalk among DNA methylation, copy number variation and miRNA expression for immune regulation and evasion. The repertoire of immune-related miRNAs will facilitate the development of immunotherapeutic targets for glioma in future.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eAuthors\u0026rsquo; contributions\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eZ.Y, K.L and H.Y conceived the study. X.L and\u0026nbsp;H.X\u0026nbsp;performed experimental validation analysis. A.E.T supervised data analysis.\u0026nbsp;L.X and J.L performed data analysis.\u0026nbsp;M.F, J.L, H.Z, Y.W, L.Z, G.Z and Y.H collected and processed samples.\u0026nbsp;Z.Y wrote the paper. A.E.T revised the paper. All authors read and approved the paper.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFunding\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe study was supported by the National Natural Science Foundation of China (91959106, 32071270,\u0026nbsp;82072370, 82103376); Outstanding Innovative Research Team for Molecular Enzymology and Detection in Anhui Provincial Universities (2022AH010012); Program for Excellent Sci-tech Innovation Teams of Universities in Anhui Province (2022AH010074); the Major Science and Technology Projects in Anhui Province (202003a06020009); Anhui Natural Science Foundation (2208085MC48); Key University Science Research Project of Anhui Province (2023AH051746, 2023AH051768); Climbing Peak Training Program for Innovative Technology team of Yijishan Hospital, Wannan Medical College (PF201904); Peak Training Program for Scientific Research of Yijishan Hospital, Wannan Medical College (GF2019T01,\u0026nbsp;GF2019G15); the Open Project of Key Laboratory of Non-coding RNA Transformation Research of Anhui Higher Education Institution, Wannan Medical College (RNA202205);\u0026nbsp;Science and Technology Application Basic Research Project of Wuhu (2022jc60); and Key Research and Development Foundation supported by Science and Technology Department of Sichuan Province (2023YFS0243). The sponsors have no role in study design, in the collection, analysis and interpretation of data, in the writing of the report, and in the decision to submit the article for publication.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAvailability of data and materials\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eData analyzed in this manuscript is already publicly available from The Cancer Genome Atlas (TCGA) data portal: https://portal.gdc.cancer.gov/ and the Chinese Glioma Genome Atlas (CGGA) data portal: http://www.cgga.org.cn/.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eEthics approval and consent to participate\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eExperimental methods using human glioma tissues were approved by the Medical Ethics Committee of Wannan Medical College and are in compliance with the Helsinki Declaration of the world medical association. All of the donors gave written informed consent for their sample donation could also be used for scientific investigations.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eConsent for publication\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNot applicable.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCompeting interests\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe authors declare that they have no competing interests.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\u003cli\u003e\u003cspan\u003eOstrom QT, Bauchet L, Davis FG, Deltour I, Fisher JL, Langer CE, et al. The epidemiology of glioma in adults: a state of the science review. Neuro Oncol. 2014;16(7):896\u0026ndash;913.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGoodenberger ML, Jenkins RB. Genetics of adult glioma. Cancer Genet. 2012;205(12):613\u0026ndash;21.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWeller M, Stupp R, Reifenberger G, Brandes AA, van den Bent MJ, Wick W, et al. MGMT promoter methylation in malignant gliomas: ready for personalized medicine? Nat Rev Neurol. 2010;6(1):39\u0026ndash;51.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003evan den Bent MJ. Chemotherapy for low-grade glioma: when, for whom, which regimen? Curr Opin Neurol. 2015;28(6):633\u0026ndash;938.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eYang K, Wu Z, Zhang H, Zhang N, Wu W, Wang Z, et al. Glioma targeted therapy: insight into future of molecular approaches. Mol Cancer. 2022;21(1):39.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLouis DN, Perry A, Reifenberger G, von Deimling A, Figarella-Branger D, Cavenee WK, et al. The 2016 World Health Organization Classification of Tumors of the Central Nervous System: a summary. Acta Neuropathol. 2016;131(6):803\u0026ndash;20.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004;116(2):281\u0026ndash;97.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eO'Brien J, Hayder H, Zayed Y, Peng C. Overview of MicroRNA Biogenesis, Mechanisms of Actions, and Circulation. Front Endocrinol (Lausanne). 2018;9:402.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eShivdasani RA. MicroRNAs: regulators of gene expression and cell differentiation. Blood. 2006;108(12):3646\u0026ndash;53.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHwang HW, Mendell JT. MicroRNAs in cell proliferation, cell death, and tumorigenesis. Br J Cancer. 2007;96 Suppl:R40-4.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTufekci KU, Oner MG, Meuwissen RL, Genc S. The role of microRNAs in human diseases. Methods Mol Biol. 2014;1107:33\u0026ndash;50.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePaul P, Chakraborty A, Sarkar D, Langthasa M, Rahman M, Bari M, et al. Interplay between miRNAs and human diseases. J Cell Physiol. 2018;233(3):2007\u0026ndash;18.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eOhno M, Matsuzaki J, Kawauchi J, Aoki Y, Miura J, Takizawa S, et al. Assessment of the Diagnostic Utility of Serum MicroRNA Classification in Patients With Diffuse Glioma. JAMA Netw Open. 2019;2(12):e1916953.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhang Y, Dutta A, Abounader R. The role of microRNAs in glioma initiation and progression. Front Biosci (Landmark Ed). 2012;17(2):700\u0026ndash;12.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi Y, Xu J, Chen H, Bai J, Li S, Zhao Z, et al. Comprehensive analysis of the functional microRNA-mRNA regulatory network identifies miRNA signatures associated with glioma malignant progression. Nucleic Acids Res. 2013;41(22):e203.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWang Q, Li P, Li A, Jiang W, Wang H, Wang J, et al. Plasma specific miRNAs as predictive biomarkers for diagnosis and prognosis of glioma. J Exp Clin Cancer Res. 2012;31(1):97.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhang Y, Chen J, Xue Q, Wang J, Zhao L, Han K, et al. Prognostic Significance of MicroRNAs in Glioma: A Systematic Review and Meta-Analysis. Biomed Res Int. 2019;2019:4015969.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTumilson CA, Lea RW, Alder JE, Shaw L. Circulating microRNA biomarkers for glioma and predicting response to therapy. Mol Neurobiol. 2014;50(2):545\u0026ndash;58.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMondal I, Kulshreshtha R. Potential of microRNA based diagnostics and therapeutics in glioma: a patent review. Expert Opin Ther Pat. 2021;31(1):91\u0026ndash;106.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eAhmadpour S, Taghavi T, Sheida A, Tamehri Zadeh SS, Hamblin MR, Mirzaei H. Effects of microRNAs and long non-coding RNAs on chemotherapy response in glioma. Epigenomics. 2022;14(9):549\u0026ndash;63.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eXiao Y, Bi M, Guo H, Li M. Multi-omics approaches for biomarker discovery in early ovarian cancer diagnosis. EBioMedicine. 2022;79:104001.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHeo YJ, Hwa C, Lee GH, Park JM, An JY. Integrative Multi-Omics Approaches in Cancer Research: From Biological Networks to Clinical Subtypes. Mol Cells. 2021;44(7):433\u0026ndash;43.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePos O, Radvanszky J, Buglyo G, Pos Z, Rusnakova D, Nagy B, et al. DNA copy number variation: Main characteristics, evolutionary significance, and pathological aspects. Biomed J. 2021;44(5):548\u0026ndash;59.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKrijgsman O, Carvalho B, Meijer GA, Steenbergen RD, Ylstra B. Focal chromosomal copy number aberrations in cancer-Needles in a genome haystack. Biochim Biophys Acta. 2014;1843(11):2698\u0026ndash;704.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDe Carvalho DD, You JS, Jones PA. DNA methylation and cellular reprogramming. Trends Cell Biol. 2010;20(10):609\u0026ndash;17.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBlattler A, Farnham PJ. Cross-talk between site-specific transcription factors and DNA methylation states. J Biol Chem. 2013;288(48):34287\u0026ndash;94.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBorgel J, Guibert S, Li Y, Chiba H, Schubeler D, Sasaki H, et al. Targets and dynamics of promoter DNA methylation during early mouse development. Nat Genet. 2010;42(12):1093\u0026ndash;100.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eFialkova V, Vidomanova E, Balharek T, Marcinek J, Kudela E, Hanysova S, et al. DNA methylation as mechanism of apoptotic resistance development in endometrial cancer patients. Gen Physiol Biophys. 2017;36(5):521\u0026ndash;9.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSun W, Bunn P, Jin C, Little P, Zhabotynsky V, Perou CM, et al. The association between copy number aberration, DNA methylation and gene expression in tumor samples. Nucleic Acids Res. 2018;46(6):3009\u0026ndash;18.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eAli Syeda Z, Langden SSS, Munkhzul C, Lee M, Song SJ. Regulatory Mechanism of MicroRNA Expression in Cancer. Int J Mol Sci. 2020;21(5).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eCancer Genome Atlas, Research N, Brat DJ, Verhaak RG, Aldape KD, Yung WK, Salama SR, et al. Comprehensive, Integrative Genomic Analysis of Diffuse Lower-Grade Gliomas. N Engl J Med. 2015;372(26):2481\u0026ndash;98.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDi Lena P, Sala C, Prodi A, Nardini C. Missing value estimation methods for DNA methylation data. Bioinformatics. 2019;35(19):3786\u0026ndash;93.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTeschendorff AE, Marabita F, Lechner M, Bartlett T, Tegner J, Gomez-Cabrero D, et al. A beta-mixture quantile normalization method for correcting probe design bias in Illumina Infinium 450 k DNA methylation data. Bioinformatics. 2013;29(2):189\u0026ndash;96.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhao Z, Zhang KN, Wang Q, Li G, Zeng F, Zhang Y, et al. Chinese Glioma Genome Atlas (CGGA): A Comprehensive Resource with Functional Genomic Data from Chinese Glioma Patients. Genomics Proteom Bioinf. 2021;19(1):1\u0026ndash;12.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eXu F, Wang Y, Ling Y, Zhou C, Wang H, Teschendorff AE, et al. dbDEMC 3.0: Functional Exploration of Differentially Expressed miRNAs in Cancers of Human and Model Organisms. Genomics Proteom Bioinf. 2022;20(3):446\u0026ndash;54.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKozomara A, Birgaoanu M, Griffiths-Jones S. miRBase: from microRNA sequences to function. Nucleic Acids Res. 2019;47(D1):D155\u0026ndash;D62.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMermel CH, Schumacher SE, Hill B, Meyerson ML, Beroukhim R, Getz G. GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers. Genome Biol. 2011;12(4):R41.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ede Rie D, Abugessaisa I, Alam T, Arner E, Arner P, Ashoor H, et al. An integrated expression atlas of miRNAs and their promoters in human and mouse. Nat Biotechnol. 2017;35(9):872\u0026ndash;8.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMirzal A. Nonparametric Tikhonov Regularized NMF and Its Application in Cancer Clustering. IEEE/ACM Trans Comput Biol Bioinform. 2014;11(6):1208\u0026ndash;17.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLu X, Meng J, Zhou Y, Jiang L, Yan F. MOVICS: an R package for multi-omics integration and visualization in cancer subtyping. Bioinformatics. 2021;36(22\u0026ndash;23):5539\u0026ndash;41.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePierre-Jean M, Deleuze JF, Le Floch E, Mauger F. Clustering and variable selection evaluation of 13 unsupervised methods for multi-omics data integration. Brief Bioinform. 2020;21(6):2011\u0026ndash;30.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, et al. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res. 2020;48(W1):W509\u0026ndash;W14.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eYoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBhattacharya S, Dunn P, Thomas CG, Smith B, Schaefer H, Chen J, et al. ImmPort, toward repurposing of open access immunological assay data for translational and clinical research. Sci Data. 2018;5:180015.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eRooney MS, Shukla SA, Wu CJ, Getz G, Hacohen N. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell. 2015;160(1\u0026ndash;2):48\u0026ndash;61.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eNarayanan S, Kawaguchi T, Yan L, Peng X, Qi Q, Takabe K. Cytolytic Activity Score to Assess Anticancer Immunity in Colorectal Cancer. Ann Surg Oncol. 2018;25(8):2323\u0026ndash;31.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLauss M, Donia M, Harbst K, Andersen R, Mitra S, Rosengren F, et al. Mutational and putative neoantigen load predict clinical benefit of adoptive T cell therapy in melanoma. Nat Commun. 2017;8(1):1738.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLiu Y, Liang G, Xu H, Dong W, Dong Z, Qiu Z, et al. Tumors exploit FTO-mediated regulation of glycolytic metabolism to evade immune surveillance. Cell Metab. 2021;33(6):1221\u0026ndash;33e11.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eJu M, Bi J, Wei Q, Jiang L, Guan Q, Zhang M et al. Pan-cancer analysis of NLRP3 inflammasome with potential implications in prognosis and immunotherapy in human cancer. Brief Bioinform. 2021;22(4).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eThorsson V, Gibbs DL, Brown SD, Wolf D, Bortone DS, Ou Yang TH, et al. The Immune Landscape of Cancer. Immunity. 2018;48(4):812\u0026ndash;30. e14.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eRu B, Wong CN, Tong Y, Zhong JY, Zhong SSW, Wu WC, et al. TISIDB: an integrated repository portal for tumor-immune system interactions. Bioinformatics. 2019;35(20):4200\u0026ndash;2.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eShinohara S, Takahashi Y, Komuro H, Matsui T, Sugita Y, Demachi-Okamura A et al. New evaluation of the tumor immune microenvironment of non-small cell lung cancer and its association with prognosis. J Immunother Cancer. 2022;10(4).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLiberzon A, Birger C, Thorvaldsdottir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 2015;1(6):417\u0026ndash;25.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eJiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med. 2018;24(10):1550\u0026ndash;8.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWang J, Li S, Wang T, Xu S, Wang X, Kong X et al. RNA2Immune: A Database of Experimentally Supported Data Linking Non-coding RNA Regulation to the Immune System. Genomics Proteom Bioinf. 2022.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSmyth GK. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. 2004;3:Article3.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eNam JW, Rissland OS, Koppstein D, Abreu-Goodger C, Jan CH, Agarwal V, et al. Global analyses of the effect of different cellular contexts on microRNA targeting. Mol Cell. 2014;53(6):1031\u0026ndash;43.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eChen Y, Wang X. miRDB: an online database for prediction of functional microRNA targets. Nucleic Acids Res. 2020;48(D1):D127\u0026ndash;D31.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHuang HY, Lin YC, Cui S, Huang Y, Tang Y, Xu J, et al. miRTarBase update 2022: an informative resource for experimentally validated miRNA-target interactions. Nucleic Acids Res. 2022;50(D1):D222\u0026ndash;D30.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innov (Camb). 2021;2(3):100141.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018;28(11):1747\u0026ndash;56.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKim H, Huang W, Jiang X, Pennicooke B, Park PJ, Johnson MD. Integrative genome analysis reveals an oncomir/oncogene cluster regulating glioblastoma survivorship. Proc Natl Acad Sci U S A. 2010;107(5):2183\u0026ndash;8.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eFavero F, McGranahan N, Salm M, Birkbak NJ, Sanborn JZ, Benz SC, et al. Glioblastoma adaptation traced through decline of an IDH1 clonal driver and macro-evolution of a double-minute chromosome. Ann Oncol. 2015;26(5):880\u0026ndash;7.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSchliesser MG, Claus R, Hielscher T, Grimm C, Weichenhan D, Blaes J, et al. Prognostic relevance of miRNA-155 methylation in anaplastic glioma. Oncotarget. 2016;7(50):82028\u0026ndash;45.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWu X, Wan Q, Wang J, Hou P, Zhang Q, Wang Q, et al. Epigenetic Activation of lncRNA MIR155HG Mediated by Promoter Hypomethylation and SP1 is Correlated with Immune Infiltration in Glioma. Onco Targets Ther. 2022;15:219\u0026ndash;35.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eVrba L, Munoz-Rodriguez JL, Stampfer MR, Futscher BW. miRNA gene promoters are frequent targets of aberrant DNA methylation in human breast cancer. PLoS ONE. 2013;8(1):e54398.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi CL, Nie H, Wang M, Su LP, Li JF, Yu YY, et al. microRNA-155 is downregulated in gastric cancer cells and involved in cell metastasis. Oncol Rep. 2012;27(6):1960\u0026ndash;6.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSandoval J, Diaz-Lagares A, Salgado R, Servitje O, Climent F, Ortiz-Romero PL, et al. MicroRNA expression profiling and DNA methylation signature for deregulated microRNA in cutaneous T-cell lymphoma. J Invest Dermatol. 2015;135(4):1128\u0026ndash;37.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBinnewies M, Roberts EW, Kersten K, Chan V, Fearon DF, Merad M, et al. Understanding the tumor immune microenvironment (TIME) for effective therapy. Nat Med. 2018;24(5):541\u0026ndash;50.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTang T, Huang X, Zhang G, Hong Z, Bai X, Liang T. Advantages of targeting the tumor immune microenvironment over blocking immune checkpoint in cancer immunotherapy. Signal Transduct Target Ther. 2021;6(1):72.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGocher AM, Workman CJ, Vignali DAA. Interferon-gamma: teammate or opponent in the tumour microenvironment? Nat Rev Immunol. 2022;22(3):158\u0026ndash;72.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKim T, Croce CM. MicroRNA: trends in clinical trials of cancer diagnosis and therapy strategies. Exp Mol Med. 2023;55(7):1314\u0026ndash;21.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHolla S, Stephen-Victor E, Prakhar P, Sharma M, Saha C, Udupa V, et al. Mycobacteria-responsive sonic hedgehog signaling mediates programmed death-ligand 1- and prostaglandin E2-induced regulatory T cell expansion. Sci Rep. 2016;6:24193.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eAudrito V, Serra S, Stingi A, Orso F, Gaudino F, Bologna C, et al. PD-L1 up-regulation in melanoma increases disease aggressiveness and is mediated through miR-17-5p. Oncotarget. 2017;8(9):15894\u0026ndash;911.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eCheng G, Li Y, Liu Z, Song X. lncRNA PSMA3-AS1 promotes the progression of non-small cell lung cancer through targeting miR-17-5p/PD-L1. Adv Clin Exp Med. 2021;30(10):1043\u0026ndash;50.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSelem NA, Nafae H, Manie T, Youness RA, Gad MZ. Let-7a/cMyc/CCAT1/miR-17-5p Circuit Re-sensitizes Atezolizumab Resistance in Triple Negative Breast Cancer through Modulating PD-L1. Pathol Res Pract. 2023;248:154579.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDu W, Chen D, Wei K, Yu D, Gan Z, Xu G, et al. MiR-10b-5p Impairs TET2-Mediated Inhibition of PD-L1 Transcription Thus Promoting Immune Evasion and Tumor Progression in Glioblastoma. Tohoku J Exp Med. 2023;260(3):205\u0026ndash;14.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhao P, Zhen H, Zhao H, Huang Y, Cao B. Identification of hub genes and potential molecular mechanisms related to radiotherapy sensitivity in rectal cancer based on multiple datasets. J Transl Med. 2023;21(1):176.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eYe J, Zeng T. Mining database and verification of PIK3CB as a marker predicting prognosis and immune infiltration in renal clear cell carcinoma. Med (Baltim). 2022;101(22):e29254.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKatoh Y, Yaguchi T, Kubo A, Iwata T, Morii K, Kato D et al. Inhibition of stearoyl-CoA desaturase 1 (SCD1) enhances the antitumor T cell response through regulating beta-catenin signaling in cancer cells and ER stress in T cells and synergizes with anti-PD-1 antibody. J Immunother Cancer. 2022;10(7).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi H, Yang Z, Yang X, Zhang F, Wang J, Wu Z, et al. LINC01123 promotes immune escape by sponging miR-214-3p to regulate B7-H3 in head and neck squamous-cell carcinoma. Cell Death Dis. 2022;13(2):109.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eYou J, Wu W, Lu M, Xie Y, Miao R, Gu M, et al. Hepatic exosomes with declined MiR-27b-3p trigger RIG-I/TBK1 signal pathway in macrophages. Liver Int. 2022;42(7):1676\u0026ndash;91.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMarcinkowska M, Szymanski M, Krzyzosiak WJ, Kozlowski P. Copy number variation of microRNA genes in the human genome. BMC Genomics. 2011;12:183.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGulyaeva LF, Kushlinskiy NE. Regulatory mechanisms of microRNA expression. J Transl Med. 2016;14(1):143.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eRoy R, Chatterjee A, Das D, Ray A, Singh R, Chattopadhyay E, et al. Genome-wide miRNA methylome analysis in oral cancer: possible biomarkers associated with patient survival. Epigenomics. 2019;11(5):473\u0026ndash;87.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eVeerappa AM, Murthy MN, Vishweswaraiah S, Lingaiah K, Suresh RV, Nachappa SA, et al. Copy number variations burden on miRNA genes reveals layers of complexities involved in the regulation of pathways and phenotypic expression. PLoS ONE. 2014;9(2):e90391.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eChu G, Ji X, Wang Y, Niu H. Integrated multiomics analysis and machine learning refine molecular subtypes and prognosis for muscle-invasive urothelial cancer. Mol Ther Nucleic Acids. 2023;33:110\u0026ndash;26.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWang B, Li M, Li R. Identification and verification of prognostic cancer subtype based on multi-omics analysis for kidney renal papillary cell carcinoma. Front Oncol. 2023;13:1169395.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLiu B, Shyr Y, Cai J, Liu Q. Interplay between miRNAs and host genes and their role in cancer. Brief Funct Genomics. 2018;18(4):255\u0026ndash;66.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLouhimo R, Hautaniemi S. CNAmet: an R package for integrating copy number, methylation and expression data. Bioinformatics. 2011;27(6):887\u0026ndash;8.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhang F, Zhang Y, Lv X, Xu B, Zhang H, Yan J, et al. Evolution of an X-Linked miRNA Family Predominantly Expressed in Mammalian Male Germ Cells. Mol Biol Evol. 2019;36(4):663\u0026ndash;78.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLin WW, Ou GY, Zhao WJ. Mutational profiling of low-grade gliomas identifies prognosis and immunotherapy-related biomarkers and tumour immune microenvironment characteristics. J Cell Mol Med. 2021;25(21):10111\u0026ndash;25.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eCheng Q, Duan W, He S, Li C, Cao H, Liu K, et al. Multi-Omics Data Integration Analysis of an Immune-Related Gene Signature in LGG Patients With Epilepsy. Front Cell Dev Biol. 2021;9:686909.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBartel DP. MicroRNAs: target recognition and regulatory functions. Cell. 2009;136(2):215\u0026ndash;33.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eQuigley DA, Kristensen V. Predicting prognosis and therapeutic response from interactions between lymphocytes and tumor cells. Mol Oncol. 2015;9(10):2054\u0026ndash;62.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePaladini L, Fabris L, Bottai G, Raschioni C, Calin GA, Santarpia L. Targeting microRNAs as key modulators of tumor immune response. J Exp Clin Cancer Res. 2016;35:103.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePathania AS, Prathipati P, Olwenyi OA, Chava S, Smith OV, Gupta SC, et al. miR-15a and miR-15b modulate natural killer and CD8(+)T-cell activation and anti-tumor immune response by targeting PD-L1 in neuroblastoma. Mol Ther Oncolytics. 2022;25:308\u0026ndash;29.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi X, Wang S, Mu W, Barry J, Han A, Carpenter RL, et al. Reactive oxygen species reprogram macrophages to suppress antitumor immune response through the exosomal miR-155-5p/PD-L1 pathway. J Exp Clin Cancer Res. 2022;41(1):41.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTurcan S, Rohle D, Goenka A, Walsh LA, Fang F, Yilmaz E, et al. IDH1 mutation is sufficient to establish the glioma hypermethylator phenotype. Nature. 2012;483(7390):479\u0026ndash;83.\u003c/span\u003e\u003c/li\u003e\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"Lower grade glioma, microRNA, DNA methylation, Copy number variation, Omics-data integration, Tumor immune microenvironment, Regulatory network","lastPublishedDoi":"10.21203/rs.3.rs-3935802/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-3935802/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003ch2\u003eBackground\u003c/h2\u003e \u003cp\u003eLower grade glioma (LGG) is the most frequent primary tumors of the central nervous system and has been a major healthcare burden, however, the specific molecular mechanism underlying its initiation and progression remains to be elucidated. Although it is known that microRNAs (miRNAs) are widely involved in the regulation of malignant phenotypes of glioma, the underling mechanism for miRNA dysregulation remains largely unanswered.\u003c/p\u003e\u003ch2\u003eMethods\u003c/h2\u003e \u003cp\u003eIn the present work, we developed a novel strategy to obtain the genome wide copy number variation (CNV) and promoter DNA methylation (DNAm) data of miRNAs and performed a systematic integrative study for the multi-omics data to identify mechanisms underlying miRNA dysregulation molecular subtyping in LGG. The relationship between LGG subtypes, prognosis, molecular features, tumor immune microenvironment and response to immune therapy was further analyzed. We also developed a prognostic model based on immune-related miRNAs that were differentially expressed between LGG samples. Then, the influence of the prognostic model on the immune microenvironment in LGG was comprehensively analyzed.\u003c/p\u003e\u003ch2\u003eResults\u003c/h2\u003e \u003cp\u003eWe identify 719 miRNAs whose expression was associated with alterations of copy number variation or promoter methylation. Integrative multi-omics analysis of the copy number and methylation related miRNAs revealed four subtypes with differing prognoses, which were validated with independent cohort data. These glioma subtypes exhibited distinct immune-related characteristics as well as clinical and genetic features. We further screened immune-related miRNAs through investigation of their correlation with immune cell infiltrations and immune microenvironment. By construction of a miRNA regulatory network, we identified candidate miRNAs associated with immune evasion and response to glioma immunotherapy. We finally evaluated the associations between prognosis related miRNAs and immune regulation. Among them, miR-155-5p, miR-196a-5p, miR-196b-5p, miR-200a-3p, miR-503-5p, and miR-15b-5p were validated as immunoevasive biomarkers and to promote cell migration, invasion and proliferation for glioma through \u003cem\u003ein vitro\u003c/em\u003e experiments.\u003c/p\u003e\u003ch2\u003eConclusions\u003c/h2\u003e \u003cp\u003eOur study systematically reveals the crosstalk among DNA methylation, copy number variation and miRNA expression for immune regulation in glioma, and could have important implications for patient stratification and development of novel biomarkers for immunotherapy approaches.\u003c/p\u003e","manuscriptTitle":"Integrative analysis of genomic and epigenomic regulation reveals microRNA regulatory network mediated tumor heterogeneity and immune evasion in lower grade glioma","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2024-02-09 21:19:58","doi":"10.21203/rs.3.rs-3935802/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"d2273d40-1073-411c-abbe-1d0558d6b2fe","owner":[],"postedDate":"February 9th, 2024","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[],"tags":[],"updatedAt":"2024-02-12T00:32:05+00:00","versionOfRecord":[],"versionCreatedAt":"2024-02-09 21:19:58","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-3935802","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-3935802","identity":"rs-3935802","version":["v1"]},"buildId":"qtupq5eGEP_6zYnWcrvyt","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}