Supplementary files: 1. Supplementary materials: Figure S1: Protein-Protein Interaction Networks (PPI) of 310 genes from 20 AAMR-related gene sets; Figure S2: Heatmap and correlation analysis between TMEM251 and enrichment scores; Figure S3: Connections between hub genes and clinicopathological factors; Figure S4: Significantly enriched AAMR-related gene sets in GSEA; Figure S5: Immune infiltration heatmap of different risk groups; Figure S6: Gene mutation patterns of samples. (A) Gene mutation patterns of all samples (n = 330). (B) Gene mutation patterns of the low-risk group. (C) Gene mutation patterns of the high-risk group. Datasheet S1: Pearson correlation analysis between TMEM251 and other genes; Datasheet S2: GO analysis results (BP); Datasheet S3: GO analysis results (CC); Datasheet S4: GO analysis results (MF); Datasheet S5: KEGG analysis results; Datasheet S6: 310 research genes from the 20 enriched gene sets; Datasheet S7: Unsupervised cluster analysis based on 1234 TMEM251-coexpression genes; Datasheet S8: Unsupervised cluster analysis based on 310 genes from the 20 enriched gene sets; Datasheet S9: Unsupervised cluster analysis based on GSVA enrichment scores; Datasheet S10: Differential gene expression analysis between cancer and normal samples; Datasheet S11: The intersection of 310 genes and DEGs; Datasheet S12: LASSO regression analysis of 55 research genes; Datasheet S13: Stepwise multivariate Cox regression analysis for the predictive model; Datasheet S14: Multivariate Cox proportional hazard regression model (without N stage); Datasheet S15: Gene Set Enrichment Analysis (GSEA) results between different risk groups; Datasheet S16: Immune infiltration analysis (CIBERSORT) results; Datasheet S17: Tumor mutation burden (TMB) scores of all samples. 2. Supplementary other files: Raw data: All raw data included in this article. Code and datasheet: Software codes and associated datasheet.

Bubble plots exhibited the top-ranked GSEA results for the top-ranked KEGG pathways in (A) recurrent POP compared with primary POP uterosacral ligaments, and (B) recurrent POP compared with non-POP uterosacral ligaments. The classic GSEA plots showed that both adipose- and inflammation-related pathways were activated in the two contrast matrices for (C) recurrent POP vs.vs primary POP uterosacral ligaments, and (D) recurrent POP vs.vs non-POP uterosacral ligaments. Heatmaps showed hierarchical clustering of ssGSEA scores of KEGG pathways differentially enriched in (E) recurrent POP vs.vs primary POP uterosacral ligaments and (F) recurrent POP vs.vs non-POP uterosacral ligaments. (G) Spearman correlation analysis of pathway ssGSEA scores revealed that PPAR signaling pathways was strongly associated with adipose- and inflammation-related pathways in non-POP, primary POP and recurrent POP uterosacral ligaments. GSEA, g: Gene set enrichment analysis; ssGSEA, the s: Single-sample gene set enrichment analysis; POP, p: Pelvic organ prolapses; NES, n: Normalized enrichment scores.

Object: To identify genes of prognostic value which associated with tumor microenvironment (TME) in acute myeloid leukemia (AML).Methods and Materials: Level 3 AML patients gene transcriptome profiles were downloaded from The Cancer Genome Atlas (TCGA) database. Clinical characteristics and survival data were extracted from the Genomic Data Commons (GDC) tool. Then, limma package was utilized for normalization processing. ESTIMATE algorithm was used for calculating immune, stromal and ESTIMATE scores. We examined the distribution of these scores in Cancer and Acute Leukemia Group B (CALGB) cytogenetics risk category. Kaplan-Meier (K-M) curves were used to evaluate the relationship between immune scores, stromal scores, ESTIMATE scores and overall survival. We performed clustering analysis and screened differential expressed genes (DEGs) by using heatmaps, volcano plots and Venn plots. After pathway enrichment analysis and gene set enrichment analysis (GESA), protein-protein interaction (PPI) network was constructed and hub genes were screened. We explore the prognostic value of hub genes by calculating risk scores (RS) and processing survival analysis. Finally, we verified the expression level, association of overall survival and gene interactions of hub genes in the Vizome database.Results: We enrolled 173 AML samples from TCGA database in our study. Higher immune score was associated with higher risk rating in CALGB cytogenetics risk category (P = 0.0396) and worse overall survival outcomes (P = 0.0224). In Venn plots, 827 intersect genes were screened with differential analysis. Functional enrichment clustering analysis revealed a significant association between intersect genes and the immune response. After PPI network, 18 TME-related hub genes were identified. RS was calculated and the survival analysis results revealed that high RS was related with poor overall survival (P < 0.0001). Besides, the survival receiver operating characteristic curve (ROC) showed superior predictive accuracy (area under the curve = 0.725). Finally, the heatmap from Vizome database demonstrated that 18 hub genes showed high expression in patient samples.Conclusion: We identified 18 TME-related genes which significantly associated with overall survival in AML patients from TCGA database.

The mechanisms underlying ETS-driven prostate cancer initiation and progression remain poorly understood due to a lack of model systems that recapitulate this phenotype. We generated a genetically engineered mouse with prostate-specific expression of the ETS factor, ETV4, at lower and higher protein dosages through mutation of its degron. Lower-level expression of ETV4 caused mild luminal cell expansion without histologic abnormalities and higher-level expression of stabilized ETV4 caused prostatic intraepithelial neoplasia (mPIN) with 100% penetrance within 1 week. Tumor progression was limited by p53-mediated senescence and Trp53 deletion cooperated with stabilized ETV4. The neoplastic cells expressed differentiation markers such as Nkx3.1 recapitulating luminal gene expression features of untreated human prostate cancer. Single-cell and bulk RNA-sequencing showed stabilized ETV4 induced a novel luminal-derived expression cluster with signatures of the cell cycle, senescence, and epithelial to mesenchymal transition. These data suggest that ETS overexpression alone, at sufficient dosage, can initiate prostate neoplasia. Methods Mouse prostate digestion: Intraperitoneal injection of tamoxifen was administered in 8-week-old mice. 2 weeks after tamoxifen treatment, the mouse prostate was digested 1 hour with Collagenase/Hyaluronidase (STEMCELL, #07912), and then 30 minutes with TrypLETM Express Enzyme (Thermo Fischer, # 12605028) at 37°C to isolate single prostate cells. The prostate cells were stained with PE/Cy7 conjugated anti-mouse CD326 (EpCAM) antibody (BioLegend, 118216) and then, CD326 and EYFP double positive cells were sorted out by flow cytometry, which are luminal cells mainly from the anterior prostate and dorsal prostate. The mRNA or genomic DNA were extracted from these double-positive cells and then were used for ATAC-sequencing and RNA-sequencing analysis. ATAC-seq and primary data processing: ATAC-seq was performed as previously described. Primary data processing and peak calling were performed using ENCODE ATAC-seq pipeline (https://github.com/kundajelab/atac_dnase_pipelines). Briefly, paired-end reads were trimmed, filtered, and aligned against mm9 using Bowtie2. PCR duplicates and reads mapped to mitochondrial chromosome or repeated regions were removed. Mapped reads were shifted +4/-5 to correct for the Tn5 transposase insertion. Peak calling was performed using MACS2, with p-value < 0.01 as the cutoff. Reproducible peaks from two biological replicates were defined as peaks that overlapped by more than 50%. On average 25 million uniquely mapped pairs of reads were remained after filtering. The distribution of inserted fragment length shows a typical nucleosome banding pattern, and the TSS enrichment score (reads that are enriched around TSS against background) ranges between 28 and 33, suggesting the libraries have high quality and were able to capture the majority of regions of interest. Differential peak accessibility: Reads aligned to peak regions were counted using R package GenomicAlignments_v1.12.2. Read count normalization and differential accessible peaks were called with DESeq2_v1.16.1 in R 3.4.1. Differential peaks were defined as peaks with adjusted p-value < 0.01 and |log2(FC)| > 2. For visualization, coverage bigwig files were generated using bamCoverage command from deepTools2, normalizing using the size factor generated by DESeq2. The differential ATAC-seq peak density plot was generated with deepTools2, using regions that were significantly more or less accessible in ETV4AAA samples relative to EYFP samples. Motif analysis: Enriched motif was performed using MEME-ChIP 5.0.0 with differentially accessible regions in ETV4AAA relative to EYFP. ATAC-seq footprinting was performed using TOBIAS. First, ACACCorrect was run to correct Tn5 bias, followed by ScoreBigwig to calculate footprint score, and finally BindDetect to generate differential footprint across regions. RNA-seq analysis: The extracted RNA was processed for RNA-sequencing by the Integrated Genomics Core Facility at MSKCC. The libraries were sequenced on an Illumina HiSeq-2500 platform with 51 bp paired-end reads to obtain a minimum yield of 40 million reads per sample. The sequenced data were aligned using STAR v2.3 with GRCm38.p6 as annotation. DESeq2_v1.16.1 was subsequently applied on read counts for normalization and the identification of differentially expressed genes between ETV4AAA and EYFP groups, with an adjusted p-value < 0.05 as the threshold. Genes were ranked by sign(log2(FC)) * (-log(p-value)) as input for GSEA analysis using ‘Run GSEA Pre-ranked’ with 1000 permutations (48). The custom gene sets used in GSEA analysis are shown in Table S2. Unsupervised hierarchical clustering: To get an overall sample clustering as part of QC, hierarchical clustering was performed using pheatmap_v1.0.10 package in R on normalized ATAC-seq or RNA-seq data. It was done using all peaks or all genes, with Spearman or Pearson correlation as the distance metric. To have an overview of the differential gene expression from the RNA-seq data, unsupervised clustering was also performed on a matrix with all samples as columns and scaled normalized read counts of differentially expressed genes between ETV4AAA and EYFP as rows. Integrative analysis of ATAC-seq, RNA-seq, and ChIP-seq data: ERG ChIP-seq peaks were called using MACS 2.1, with an FDR cutoff of q < 10-3 and the removal of peaks mapped to blacklist regions. Reproducible peaks between two biological replicates were identified as ETV4AAA ATAC-seq peaks. ERG ChIP-seq peaks and ETV4AAA ATAC-seq peaks were considered as overlap if peak summits were within 250bp. To determine whether the overlap was significant, enrichment analysis was done using regioneR_v1.8.1 in R, which counted the number of overlapped peaks between a set of randomly selected regions in the genome (excluding blacklist regions) and the ERG-ChIP seq peaks or ETV4AAA ATAC-seq peaks. A null distribution was formed using 1000 permutation tests to compute the p-value and z-score of the original evaluation. To assign ATAC-seq peaks to genes, ChIPseeker_v1.12.1 in R was used. Each peak was unambiguously assigned to one gene with a TSS or 3’ end closest to that peak. Differential gene expression between ETV4AAA and EYFP was evaluated using log2(FC) calculated by DESeq2. p-values were estimated with Wilcoxon rank t-test and Student t-test. scRNA-sequencing: Tmprss2-CreERT2, EYFP; Tmprss2-CreERT2, ETV4WT; Tmprss2-CreERT2, ETV4AAA; and Tmprss2-CreERT2, ETV4AAA; Trp53L/L mice were euthanized 2 weeks or 4 months after tamoxifen treatment (n=3 mice for each genotype and time point). After euthanasia, the prostates were dissected out and minced with scalpel, and then processed for 1h digestion with collagenase/hyaluronidase (#07912, STEMCELL Technologies) and 30min digestion with TrypLE (#12605010, Gibco). Live single prostate cells were sorted out by flow cytometry as DAPI-. For each mouse, 5,000 cells were directly processed with 10X genomics Chromium Single Cell 3’ GEM, Library & Gel Bead Kit v3 according to manufacturer’s specifications. For each sample, 200 million reads were acquired on NovaSeq platform S4 flow cell. Reads obtained from the 10x Genomics scRNAseq platform were mapped to mouse genome (mm9) using the Cell Ranger package (10X Genomics). True cells are distinguished from empty droplets using scCB2 package. The levels of mitochondrial reads and numbers of unique molecular identifiers (UMIs) were similar among the samples, which indicates that there were no systematic biases in the libraries from mice with different genotypes. Cells were removed if they expressed fewer than 600 unique genes, less than 1,500 total counts, more than 50,000 total counts, or greater than 20% mitochondrial reads. Genes detected in less than 10 cells and all mitochondrial genes were removed for subsequent analyses. Putative doublets were removed using the Doublet Detection package. The average gene detection in each cell type was similar among the samples. Combining samples in the entire cohort yielded a filtered count matrix of 48,926 cells by 19,854 genes, with a median of 6,944 counts and a median of 1,973 genes per cell, and a median of 2,039 cells per sample. The count matrix was then normalized to CPM (counts per million), and log2(X+1) transformed for analysis of the combined dataset. The top 1000 highly variable genes were found using SCANPY (version 1.6.1) (77). Principal Component Analysis (PCA) was performed on the 1,000 most variable genes with the top 50 principal components (PCs) retained with 29% variance explained. To visualize single cells of the global atlas, we used UMAP projections (https://arxiv.org/abs/1802.03426). We then performed Leiden clustering. Marker genes for each cluster were found with scanpy.tl.rank_genes_groups. Cell types were determined using the SCSA package, an automatic tool, based on a score annotation model combining differentially expressed genes (DEGs) and confidence levels of cell markers from both known and user-defined information. Heat-map were performed for single cells based on log-normalized and scaled expression values of marker genes curated from literature or identified as highly differentially expressed. Differentially expressed genes between different clusters were found using MAST package, which were shown in heat-map. The logFC of MAST output was used for the ranked gene list in GSEA analysis (48). The custom gene sets used in GSEA analysis are shown in Table S2. Gene imputation was performed using MAGIC (Markov affinity-based graph imputation of cells) package, and imputated gene expression were used in the heatmap. Analysis of public human gene expression datasets: To analyze TP53 RNA expression in human prostate cancer samples, we obtained normalized RNA-seq data from prostate cancer TCGA (www.firebrowse.org) (3). To assess the role of TP53 loss on

Additional file 2: Fig. S1 a, Expression levels of different class I HDAC genes in different pediatric small-round-blue-cell tumors, carcinomas and normal tissues by box plot presentation using a comparative study of the amc onco-genomics software tool ( https://hgserver1.amc.nl/cgi-bin/r2/main.cgi ). The number of samples in each cohort is given in brackets. b, Differential expression levels of class I HDAC genes in primary EwS at different tumor sites by box plot presentation using the GSE63157 study set and the amc onco-genomics software tool. The number of samples in each cohort is given in brackets, ND: not determined. p-value < 0.05. c, Retroviral gene transfer of EWS-FLI1 cDNA into MSC lines L87 and V54.2 [4] results in HDAC3 and HDAC8 induction as measured via qRT-PCR, while no change of HDAC1 and HDAC2 expression was observed. Induction of EWS-FLI1-dependent EZH2 expression served as control. Fig. S2 a, Tube formation assay with the EwS cell lines CHLA-10 and SK-N-MC after incubation with 3µM MS275 or 4nM FK228 over-night compared to WT control. Both HDACi clearly enhanced endothelial differentiation potential (scale bar 0.5mm). b, Analysis of neurogenic differentiation potential of the EwS cell lines CHLA-10, EW7 and SK-N-MC treated for six days with 0.5µM MS-275 or 0.2nM FK228. The neurogenic differentiation marker GFAP (glial fibrillary acidic protein) and GAP43 (growth associated protein 43) were significantly upregulated after incubation with both HDACi as demonstrated by qRT-PCR. Fig. S3 a, Cell cycle analysis of CRISPR/Cas9 HDAC1 or HDAC2 knock outs compared to their controls (Cas9) in three different EwS cell lines are shown. Distributional analysis of cell cycle phases of HDAC1 or HDAC2 knock outs compared to their control were performed by propidium iodine staining and flow cytometry measurement, respectively. b, To analyze apoptosis in HDAC1 and HDAC2 CRISPR/Cas9 knock outs, DNA double strand breaks were measured with anti-phospho-histone H2AX-FITC conjugated mAbs and counterstained with DAPI. Left, the frequency of g-H2AX positive foci per cell was summarized in bar graphs. Right, fluorescence images show a representative experiment with HDAC1 (top) and HDAC2 (bottom) in two different EwS cell lines each, compared to their controls. Fig. S4 a, Western blot analysis of class I HDAC protein levels and their compensation in CRISPR/Cas9 HDAC1 and HDAC2 knock outs compared to their controls (Cas9). Protein levels were detected by antibodies against HDAC1, HDAC2, HDAC3 and HDAC8. b-actin or GAPDH antibodies were used as loading control. b, Heat map of 229 genes differentially expressed in three different EwS lines CHLA-10, EW7 and SK-N-MC after CRISPR/Cas9 HDAC1 knock out, are shown. Each column represents one individual array. Microarray data with their normalized fluorescent signal intensities were used (robust multichip average (RMA); GSE162786). c, Circos plots of downregulated genes (left column) and heatmaps of pathways and ontology terms the downregulated genes are enriched for (right column). The plots are based on gene lists for three EwS cell lines (CHLA-10, EW7, SK-N-MC), containing the 300 strongest downregulated genes after HDAC1 or HDAC2 knock out, respectively. The lists of downregulated genes for HDAC knock out effects in the top row is based on averaged expression data from HDAC1 and HDAC2 knock outs. The circos plots show overlaps in the gene sets, where each gene is a spot on the inner arc. Purple lines indicate genes shared by the gene lists, and blue lines indicate functional overlaps in the lists. A blue line connects two different genes belonging to the same enriched ontology term. The strongest enriched ontology terms are depicted in the heatmaps. The cells are colored by p-value. Grey cells indicate that a term is not significantly enriched in a gene list. Hence, the heatmap shows common and unique enrichments for the three cell lines. Fig. S5 a, HDAC3 or HDAC8 expression after transient shRNA transfection measured by qRT-PCR in EwS cell lines CHLA-10, SK-N-MC or EW7, respectively. Induction of three different shRNAs was done with Doxycycline for 72 hours. b, Proliferation of EwS cells after transfection with HDAC3 (top) or HDAC8 (bottom, left) specific shRNA. Further proliferation of SK-N-MC HDAC1 knock out cells with transient HDAC3 or HDAC8 knock down (bottom, right). Control cells were transfected with irrelevant shRNA. Proliferation and cell impendence was measured by the xCELLigence assay every 4 hours. Data are shown as mean ± SEM (hexaplicates/group; p-value < 0.001, respectively < 0.0001). c, Analysis of the invasive potential of EwS cell line CHLA-10 after transient shRNA transfection with HDAC3 (top) or SK-N-MC HDAC1 knock out with HDAC8 (bottom) specific shRNA 48 hours after seeding. d, Evaluation of tumorigenicity of CRISPR/Cas9 knock outs of HDAC1 and their controls (Cas9) in EwS cell line CHLA-10. Immune deficient Rag2-/-γC-/-mice were injected s.c. with 4x106 EwS cells. Mice with an average tumor size >10 mm in diameter were considered positive and sacrificed. Kaplan-Meier plots of individual experiments with six mice per group are shown. Log-rank test was used to test for differences in survival. Fig. S6 a, Proliferation of SK-N-MC and CHLA-10 after treatment with Doxorubicin and/or HDACi (MS-275 or FK228) was analyzed with the xCELLigence system. Cell impedance was measured every 4 hours. Data are shown as mean ± SEM (hexaplicates/group; p-value < 0.0001). Fig. S7 a, Proliferation of EwS CRISPR/Cas9 HDAC 1 knock outs and their controls (Cas9) in CHLA-10 or SK-N-MC cells after treatment with Vincristine (top 2 panels) or combined treatment of SK-N-MC with MS-275 and Vincristine (bottom panel). Proliferation and cell impendence were measured by the xCELLigence assay every 4 hours. Data are shown as mean ± SEM (hexaplicates/group; p-value > 0.0001). b, Heatmaps of pathways and ontology terms that are enriched among up- and downregulated genes. The plots are based on gene lists for two EwS cell lines (EW7, SK-N-MC), containing the 300 strongest differentially expressed genes after FK228, Vincristine or combined treatment, compared to solvent controls, respectively. The strongest enriched ontology terms are depicted in the heatmaps. The cells are colored by p-value. Grey cells indicate that a term is not significantly enriched in a gene list. Hence, the heatmap shows common and unique enrichments for the two cell lines. c, Spheroid growth was monitored in Greiner bio-one CELLSTAR® Cell-Repellent Surface 96-well round bottom plates. Left, CHLA-10 or EW7 cells were plated in Matrigel-containing medium and cells were treated for 48 hours with the inhibitors as indicated. Results were compared to solvent controls. Right, primary EwS tumor cells derived from PDX mice. Cell viability was measured with CellTiter Glo® Luminescent assay (quadruplets/group). Fig. S8 a, Western blot analysis of apoptosis susceptibility after FK228 or MS-275 and/or A-395 treatment, respectively. Protein levels measured by antibodies against, PARP, CASP3, and GAPDH as loading control. CHLA-10, EW7 or SK-N-MC cells were treated for 48 hours with inhibitors. b, Left, heat map of 824 genes, 3-fold differentially expressed in different EwS tumor samples (CHLA-10 and SK-N-MC) at the end of treatment, are shown. Right, zoomed in heat map with 132 genes contains only those genes with a p-value < 0.05. Each column represents one individual array. Microarray data with their normalized fluorescent signal intensities were used (robust multichip average (RMA); GSE162788). Cells were treated for 27 hours with solvent control or EEDi (A-395), HDACi (FK228) or with both inhibitors. c, Volcano plot of differentially expressed genes of EwS cells at the end of treatment (CHLA-10, SK-N-MC). The plot shows fold changes of log2 expression values (log FC) and p-values obtained from differential expression analysis comparing tumors treated with A-395 + FK228 to solvent controls. Depicted in red are genes obtaining p-value < 0.05 and absolute log FC > 1; in blue, genes with p-value < 0.05 and absolute log FC ≤ 1; in green, genes with p-value ≥ 0.05 and absolute log FC > 1; and in black, genes with p-value ≥ 0.05 and absolute log FC ≤ 1. Positive log FCs indicate higher expression of the gene in the treated cell lines. D, GSEA enrichment plots of up- and downregulated gene sets after combined A-395 and FK228 treatment. NES: normalized enrichment score. GSEA: http://www.broadinstitute.org/gsea/index.jsp .

Additional file 1: Figure S1. Immune cell expression heterogeneity and cell subsets distribution across patients, related to Fig. 1. (A) UMAP projection of 2068 single T cells (left), 515 B cells (middle) and 126 macrophages (right) from 19 patients. Each dot corresponds to one single cell, colored according to cell cluster. (B) Heatmap of T cell clusters (left), B cell clusters (middle) and macrophage clusters (right) with unique signature genes. Top 20 specifically expressed genes are marked alongside, if available. (C-E) Bar plots showing the number (left panel) and fraction (right panel) of cells originating from the 19 patients for each subcluster of T cells (C), B cells (D) and macrophages (E). (F) The fractions of the 15 subclusters, NK cells, CAFs and endothelial cells in each patient. Figure S2. Cell subcluster characterization of functional status. (A) Top 100 ranked (based on fold change) differentially expressed genes indicative of the functional status in each T-cell cluster (top) and z-score normalized mean expression of known functional marker sets across single T cells (bottom). The numbers in parentheses correspond to the ranks and the key markers (Table S1) are highlighted by red color. (B) Heatmap showing the log2-transformed expression of selected T cell function-associated genes in single cells. (C) Violin plots showing the expression profile of selected genes involved in T-cell cytotoxicity (top) and exhaustion (bottom), stratified by T-cell clusters. (D) Top 100 ranked (based on fold change) differentially expressed genes indicative of the functional status in each cluster (C1, C2 and C3 for B cells; C0 and C1 for macrophages). The numbers in parentheses correspond to the ranks and the key markers (Table S1) are highlighted by red color. (E) Z-score normalized mean expression of known functional marker sets across single B cells (top) and the log2-transformed expression of selected B cell function-associated genes in single cells (bottom). (F-G) Heatmaps showing the z-score normalized mean expression of known functional marker sets across single macrophages and their log2-transformed expression in single cells. Blue boxes highlight the key markers and the numbers in brackets represent the total times appeared in literature. Figure S3. MM17 reference profile and performance assessment. (A) Heatmap of MM17 reference profile depicting z-score normalized expression of each gene across 17 tumor microenvironment (TME) cell subsets. (B-C) Correlation between predicted proportions and true proportions for each individual cell state (B) and for each individual patient (C). (D) Confusion matrix of all TME cell states. Figure S4. Functional associations of tumor microenvironment (TME) cell states. (A-D) Enriched GO biological processes of T_CD8_Cytotoxic (A), B_Non-regulatory (B), T_CD8_Mixed (C) and CAF (D) based on gene set enrichment analysis (GSEA). Figure S5. Associations between cell states and clinico-pathological variables. (A-C) Associations of molecular and clinical features with cell states. (A) Boxplots showing the cell fraction distribution of each cell state stratified by tumor type (left), gender (middle) and tumor status (right). (B) Boxplots showing the cell fraction distribution of each cell state stratified by integrative age (left), tumor stage (middle), and race (right). (C) The fraction distribution of cell states stratified by TCGA subtypes. Median value difference of cell fraction among subtypes was evaluated using Mood’s test. Wilcoxon rank sum tests were used to examine the significance of the differences between two groups. For tumor stage, patients with Stage 0, Stage I, IA, IB, Stage II, IIA, IIB and IIC are grouped as “LOW” (n=154), Stage III, IIIA, IIIB, IIIC and Stage IV are grouped as “HIGH” (n=162). * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001. Figure S6. Associations between cell states and immune phenotypes, related to Fig. 4. (A) Scatterplots showing relationships between T_CD8_Cytotoxic and M_M2 (top), B_Regulatory and T_CD4_Exhausted (middle), CAF and T_CD8_Mixed (Cytotoxic and Exhausted) (bottom). Pearson correlations and p values are indicated. For significant correlations, linear models are shown as blue lines. (B) Contributions of the cell states to CA-1 (top) and CA-2 (bottom). (C) Scatter chart of the Pearson correlations of CA-1 and CA-2 with cell states. Different colors indicate whether or not significant associations between CA scores and cell states were observed (p < 0.05). (D) Boxplots showing the cell fraction distribution of each cell state stratified by the median values of CA-1 (top) and CA-2 (bottom), respectively. Wilcoxon rank sum tests were used to examine the significance of the differences between two groups. (E) The distribution of cell states across the three immunophenotype groups classified by median values of CA-1 and CA-2. Median value difference of cell fraction among groups was evaluated using Mood’s test. Then the statistical significance between any two groups was evaluated by Wilcoxon rank sum test and p values are shown at the top of each panel. * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001. Figure S7. Assessment on association between tumor microenvironment immune phenotypes (TMIP) and response to immune checkpoint blockade (ICB) in melanoma. (A) Box plots showing differences of CA-1 (upper panel) and CA-2 (middle panel) scores between responders and non-responders in patients under immunotherapy in TCGA data. Bar charts showing numbers of responders and non-responders with different TMIPs in those patients (lower panel). (B) Projection of each patient of Riaz et al. dataset onto the first and second component of the correspondence analysis. Left panel showed pre-treatment samples and right panel denoted on-treatment patients. Non-responders were colored blue, and responders were colored orange. Points denoted Ipi-naive patients, and triangles denoted Ipi-progressed patients. (C) Box plots showing differences of CA-2 scores between responders and non-responders in anti-PD1 pre-treatment patients (upper panel) and on-treatment patients (lower panel) who progressed after a first-line anti-CTLA4 treatment (Ipi-progressed) in Riaz et al. data. (D-E) Comparison of each cell state proportion between responders and non-responders in Ipi-progressed patients based on pre-treatment (D) and on-treatment (E) transcriptomic profiles. ns: not significant; *: p < 0.05. Table S1. Gene lists used for functional analyses. Table S3. Demographics and characteristics of patients with melanoma. Table S4. Uni- and multivariate analysis for progress-free survival (316 sample). Table S5. Uni- and multivariate analysis for overall survival (316 sample).

Supplementary Material 1. Supplementary Fig. 1. Quality control and baseline data of each enrolled sample. (A). Principal component analysis before and after processing in Harmony. (B). t-SNE projections among different groups. (C). Quality control of the scRNA-seq data. (D). t-SNE projections from all enrolled samples. t-SNE in the control group (left), STS Pre group (middle), and STS Post group (right). (E). Boxplot comparing the proportion of plasmacytoid dendritic cells(pDCs) across the groups. The STS Pre vs. CT and STS Post vs. CT sample comparisons show exact P values determined by the Wilcoxon rank-sum test. Pre- vs. post-STS scores were calculated via the paired two-sample Wilcoxon signed-sum test. (F). The baseline information of patients and healthy controls enrolled in the scRNA-seq cohort. Supplementary Fig. 2. Focused analysis of T cells and pDCs. (A). UMAP embedding of T lymphocytes from all profiled samples in different groups. UMAP in the control group (left), STS Pre group (middle), and STS Post group (right). (B). Boxplot comparing the proportions of CRIP + CD4 + T cells, NK T cells, and TRGC2 + CD8 + T cells across the groups. The exact P values determined by the Wilcoxon rank-sum test are shown for the STS Pre vs. CT and STS Post vs. CT comparisons. Differences between STS Pre and STS Post were evaluated by the paired two-sample Wilcoxon signed-sum test. (C). Enriched pathways from Gene Ontology Biological Process Enrichment Analysis for TRGC + CD8 + T cells. (D). Enriched pathways identified by Gene Ontology Biological Process enrichment analysis in NEAT + T cells. (E). Heatmap representing the enrichment of MSigDB Hallmark gene sets for each T lymphocyte subtype across groups. (F). Pseudotime trajectory analysis of CD4 + T lymphocyte subtypes. (G). Heatmap represents DEGs within pDCs across groups. (H). Heatmap representing the enrichment of MSigDB Hallmark gene sets in the MSigDB of each group within pDCs. Supplementary Fig. 3. Focused analysis of B cells and myeloid cells. (A). Heatmap representing the enrichment of Hallmark gene sets in the MSigDB for each cell type within B lymphocytes across groups. (B). Boxplot comparing the proportions of myeloid cells across the groups. The exact P values determined by the Wilcoxon rank-sum test are shown for the STS Pre vs. CT and STS Post vs. CT comparisons. Differences between STS Pre and STS Post were evaluated by the paired two-sample Wilcoxon signed-sum test. (C). Enriched pathways from Gene Ontology Biological Process Enrichment Analysis for CD16 + monocytes. (D). Heatmap representing the enrichment of MSigDB Hallmark gene sets in each monocyte cell type across groups. Supplementary Fig. 4. The characteristics of IFN-related genes involved in pathogenesis. (A). Heatmap showing the differentially expressed genes (DEGs) in classical dendritic cells(cDCs) across groups. (B). Heatmap showing the genes differentially expressed in mast cells across groups. (C). Heatmap representing the enrichment of MSigDB Hallmark gene sets in the mast cells across groups. (D). Heatmap representing the enrichment of MSigDB Hallmark gene sets in the cDC across groups. (E). Heatmap representing the enrichment of Hallmark gene sets in the MSigDB for each cell type within Natural killer (NK) cells across groups. (F). Venn plot of the overlapping genes downregulated in the STS Pre group among B lymphocytes, T lymphocytes, monocytes, NK cells, cDCs and pDCs. (G). Receiver operating characteristic (ROC) curve and area under the curve (AUC) of overlapping genes expressed at lower levels before treatment in all cell types. (H). T-SNE analysis of CXCR4 expression in the three groups. (I). The relative expression of CXCR4 across groups was determined through qPCR. Statistical significance is denoted as P  0.05. Supplementary Fig. 7. Supplementary analysis from an extra INS cohort (GEO233277) also validates the activation of IFN. (A). UMAP dimensionality reduction embedding from GEO datasets. (B). Heatmap showing the expression levels of the markers across each cell type using scRNAseq from the GEO datasets. The color intensity indicates the marker of interest. (C). Violin plot of the ISGs across groups using scRNAseq from the GEO datasets. Significance was evaluated with the Wilcoxon rank-sum test. (D). ISG scores among cell subtypes across groups using scRNAseq from the GEO datasets. Significance was evaluated with the Wilcoxon rank-sum test. Statistical significance is denoted as P 

(A) Heat map representation of the Arabidopsis gene set whose regulation by ARR1ΔDDK:GR or CKs depends on DELLA proteins. The colour scale represents Z-scores. (B) Enrichment of GO categories of all ARR1 target genes in the presence of DELLAs, visualized with ReviGO. (C) Gene expression analysis by RT-qPCR in response to short-term ARR1ΔDDK:GR induction with or without PAC. (D) Gene expression analysis by RT-qPCR in response to short-term induction of gai-1 with or without BA. (E) ChIP analysis of RGA::GFP-RGA at the promoters of six representative common targets for ARR1 and DELLAs. (F) Fold enrichment of GFP-RGA at selected target promoters in the presence (+DEX) vs the absence (-DEX) of ARR1ΔDDK:GR, in F1 seedlings of a cross between RGA::GFP-RGA and 35S::ARR1ΔDDK:GR plants. In this experiment, qPCR values of ChIP samples were normalized per input in each condition (-DEX, and +DEX), and here we show the ratio between those two conditions. (G) ChIP analysis of endogenous RGA at the promoters of six representative common targets for ARR1 and DELLAs, in the wild type and in arr1 arr12 mutants. ChIP was performed with anti-RGA antibodies. For (C-G), data correspond to single biological samples analyzed in triplicates. A second biological sample showed equivalent results.

Supplementary Material 1: Fig. S1. Identification of Metformin-Sensitizing Genes via CRISPR-Cas9 Screening. A, The schematic diagram showed how CDTSL identified metformin-sensitizing genes and targeted inhibitors in TNBC. B, The schematic diagram showed the composition of the CDTSL sgRNA sequences. C, The schematic diagram showed the process of CDTSL library screening. D, The Venn diagram showed how 67 candidate genes were identified through MAGeCK analysis to meet the "metformin sensitization" model. E, Functional enrichment analysis revealed a significant enrichment of histone modification-related genes. Fig. S2. The sequencing results of 1462 breast cancer patients (5 cohorts) were displayed. The scatter plot in the upper-left corner compared the expression levels of HDAC10 in tumor tissues versus normal tissues. The remaining subplots analyzed survival differences between patients with high/low HDAC10 expression groups across different cohorts (GSE9893, GSE61304, GSE42568, GSE22219, and TCGA-BRCA) using Kaplan-Meier curves, covering endpoints such as overall survival (OS), disease-free survival (DFS), relapse-free survival (RFS), and progression-free survival (PFS). The p-values from the log-rank test were also annotated. Fig. S3. Combination Efficacy of SAHA and Metformin in TNBC Cell Lines. A, The IC50 curves for SAHA (purple curve) and metformin (orange curve) were shown. The left panel displayed the percentage inhibition (%), while the right panel presented the combination index (CI) at each drug concentration. MDA-MB-231 cells were treated with SAHA, metformin, or both at the indicated concentrations. B, The IC50 curves for SAHA (purple curve) and metformin (orange curve) were shown. The left panel displayed the percentage inhibition (%), while the right panel presented the combination index (CI) at each drug concentration. Hs578T cells were treated with SAHA, metformin, or both at the indicated concentrations. Fig. S4. Colony formation and quantification of MDA-MB-231 cells treated with SAHA, metformin, and combinations. The formed colonies and quantification of MDA-MB-231 cells treated with different concentrations of SAHA (0, 0.125, 0.25, 0.5, or 1 μM), metformin (0, 2 or 4 mM), or their combinations were displayed (in three replicates). Fig. S5. We analyzed differential gene expression, enrichment, and membrane receptor intersections for SAHA and metformin treatments. A, GO-KEGG enrichment analysis of differentially expressed genes (p-value < 0.05 and |Log2FoldChange| > 1) between the SAHA-treated group and the control group. B, Venn diagram showing the intersection of differentially expressed genes between the SAHA-treated group and the control group with membrane receptor genes. C, Heatmap displaying the intersection of differentially expressed genes between the SAHA-treated group and the control group with membrane receptor genes. D, Volcano plot of the differential expression analysis between the SAHA-treated group and the control group (p-value < 0.05 and |Log2FoldChange| > 1). E, GO-KEGG enrichment analysis of differentially expressed genes (p-value < 0.05 and |Log2FoldChange| > 0.8) between the Metformin-treated group and the control group. F, Venn diagram showing the intersection of differentially expressed genes between the Metformin-treated group and the control group with membrane receptor genes. G, Heatmap displaying the intersection of differentially expressed genes between the Metformin-treated group and the control group with membrane receptor genes. H, Volcano plot of the differential expression analysis between the Metformin-treated group and the control group (p-value < 0.05 and |Log2FoldChange| > 0.8). Fig. S6. Protein quantification data from Fig. 3. A, Protein quantification from Fig. 3C showing changes in FGFR4 and STAT3 phosphorylation. MDA-MB-231 cells pretreated with JQ1 for 24 hours were exposed to SAHA for an additional 12 hours. FGFR4 and STAT3 phosphorylation changes were detected by immunoblotting. All bands were quantified from experiments repeated three times. B, Protein quantification from Fig. 3F showing changes in FGFR4 and STAT3 phosphorylation. (Left) MDA-MB-231 cells were treated with the indicated concentrations of SAHA (0, 5, 10 μM) for 12 hours. (Middle) MDA-MB-231 cells were treated with the indicated concentrations of metformin (0, 10, 20, 40 mM) for 48 hours. (Right) MDA-MB-231 cells pretreated with metformin (20 mM) for 36 hours were exposed to SAHA (5 μM) for an additional 12 hours. FGFR4 and STAT3 phosphorylation changes were detected by immunoblotting. All bands were quantified from experiments repeated three times. Fig. S7. We introduced FGFR4 siRNA with SAHA and overexpressed FGFR4 with metformin to observe apoptosis, proliferation, and protein changes. A, Cell apoptosis and growth assays. MDA-MB-231 cells were transfected with non-targeting control (NC) or FGFR4 siRNAs for 48 hours, followed by SAHA treatment for an additional 24 hours for apoptosis analysis, and 72 hours for cell growth analysis. B, Immunoblotting and protein quantification of FGFR4 and MCL-1 change. MDA-MB-231 cells pretreated with non-targeting control (NC) or FGFR4 siRNAs for 48 hours, followed by SAHA (5μM) treatment for an additional 24 hours for Immunoblotting. FGFR4 and MCL-1 change was detected by immunoblotting. All the bands were quantified from experiments repeated three times. C, Cell apoptosis and growth assays. MDA-MB-231 cells were transfected with control or FGFR4 plasmid for 48 hours, followed by metformin treatment for an additional 48 hours for apoptosis analysis, and 72 hours for cell growth analysis. D, Immunoblotting and protein quantification of FGFR4 and MCL-1 change. MDA-MB-231 cells pretreated with control or FGFR4 plasmid for 48h, followed by metformin (20mM) treatment for an additional 48 hours for Immunoblotting. FGFR4 and MCL-1 change was detected by immunoblotting. All the bands were quantified from experiments repeated three times. Fig. S8. Cell apoptosis analysis shown in Fig. S7, performed in triplicates. Fig. S9. Multidimensional analysis of FGFR4 gene necessity, expression profiles, and prognostic impact in pan-cancer and breast cancer. A, DepMap database CRISPR-Cas9 whole-genome screening results: Displaying the top 200 pan-cancer cell lines ranked by FGFR4 CERES scores, reflecting the importance of FGFR4 for cell survival. B, FGFR4 expression levels across various cancer types: Comparison between tumor and normal tissues. C, Association analysis between FGFR4 expression and various clinical features: Including PR status, ER status, HER2 status, PAM50, Pathologic T stage, Pathologic N stage, Pathologic M stage, Pathologic stage, Age, Race, Menopause status, and Histological type. D, Correlation between FGFR4 expression levels and disease-specific survival (DSS) in breast cancer patients (TCGA-BRCA dataset). E, Relationship between FGFR4 expression levels and disease-free survival (DFS) in breast cancer patients (GSE21653 dataset). F, Association between FGFR4 expression levels and relapse-free survival (RFS) in breast cancer patients (GSE9893 dataset). G, Correlation between FGFR4 expression levels and relapse-free survival (RFS) in breast cancer patients (GSE22219 dataset). H-J, Comparison of FGFR4, SLC2A1, and LDHA gene expression between cancer and normal tissues. K, Differential expression of PFKL gene between high and low FGFR4 expression groups. L, Differential expression of SLC2A1 gene between high and low FGFR4 expression groups. Fig. S10. We quantified intratumoral molecular changes using randomly selected visual fields and immunohistochemical metrics. Three visual fields were randomly selected from each group after immunohistochemistry, and the molecular changes shown in Fig. 5 were quantified using two immunohistochemical metrics.