We sourced clinical and RNA-Seq data specific to Acute Myeloid Leukemia (AML) from the TCGA dataset, accessible through the UCSC Xena portal (https://xenabrowser.net). Each tumor’s gene expression data was aligned with the human hg38 genome annotation. For subsequent analyses, we converted the gene expression measurements to Transcripts Per Kilobase Million (TPM) followed by log2 transformation ([TPM] + 1). We excluded samples that were missing either gene expression details or clinical annotations, resulting in a primary cohort of 142 AML patients. For our validation efforts, we turned to two public datasets: GSE37642 (17–20) and GSE12417 (21, 22). We incorporated only the primary AML samples with available clinical details and normalized gene expression metrics. We omitted probes without gene annotation, and for genes associated with multiple probes, the median expression was determined. This process led to the creation of validation sets with 417 and 163 AML samples, respectively. Additionally, we integrated 2293 Ligand-Receptor (LR) pairs from the connectome DB2020 database (23) detailed in ( Supplementary Table S1 ). A schematic outline of our research methodology is presented in ( Figure 1A ).

For each set of patient data, individuals were grouped into “high” or “low” categories based on the aggregate expression levels of specific ligand-receptor (LR) gene pairs. A patient was labeled as “high” if their total LR pair gene expression met or exceeded the median value for the entire patient group; otherwise, they were designated as “low.” To assess the influence of these classifications on patient outcomes, we utilized R’s “Survival” package (version 4.3.2). The log-rank test was employed to establish statistical significance, and hazard ratios (HRs) were derived using Cox regression models. Each patient group’s survival outcomes were scrutinized independently. To amalgamate insights from various cohorts, we turned to a meta-analysis approach, deploying the “Edgington” method through the “sump” function in R’s “metap” package (version 1.4). We selected the 94 significant LR pairs based on two primary criteria: 1) a Storey-adjusted q-value less than 0.1, ensuring statistical significance across multiple hypothesis tests, and 2) a consistent hazard ratio, either surpassing or falling below 1, across all evaluated cohorts, indicating a robust association with patient survival outcomes. To adjust for multiple hypothesis testing, we adopted Storey’s method (24) via the “qvalue” package in R (version 2.18.0).

Using consensus clustering, we categorized the samples based on their gene expression patterns. We generated a consistency matrix with the “ConsensusClusterPlus” package in R (25). After narrowing down to significant ligand-receptor (LR) pairs, we used these to ascertain the molecular subtypes of the samples.

For the clustering process, the “pam” algorithm was employed, and the “Canberra” method was chosen as the distance metric. We executed 500 bootstrap replications, with each cycle including 80% of the patients from the training set. We considered a cluster range from 2 to 10. The most stable clustering solution was pinpointed by examining both the consistency matrix and the consensus cumulative distribution function (CDF), as detailed by Senbabaoglu et al. (26).

To delve into the distinct gene expression landscapes across various molecular subtypes, we implemented Gene Set Enrichment Analysis (GSEA v4.0) (27, 28). We utilized comprehensive gene sets from the Hallmark database (29) and applied specific criteria for statistical significance, including a normalized p-value less than 0.01 and a False Discovery Rate (FDR) below 0.05. For functional characterization of genes that are either upregulated or downregulated, we employed the “clusterProfiler” package (30) to conduct Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses. This helps to identify the biological pathways that are most impacted by the gene expression differences in each subtype.

Regarding the immune landscape within the samples, we applied the deconvolution technique known as Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) (31). This enabled us to approximate the proportions of 22 different immune cell types within each sample. Furthermore, we used the ESTIMATE algorithm to assess the relative abundance of immune and stromal cells in the tumor microenvironment.

To evaluate the functionality and presence of T cells within tumors, we utilized the default settings of the TIDE program, which provides scores for T-cell dysfunction and exclusion (32). Additionally, we calculated the single-sample Gene Set Enrichment Analysis (ssGSEA) scores using the GSVA package in R (33), which allowed us to represent the relative enrichment levels of each KEGG pathway. By integrating these multifaceted analyses, we aim to provide a comprehensive overview of the functional roles of genes, the immune cell composition, and their potential implications for patient prognosis and treatment strategies.

For our analysis of genomic variations, we sourced the Simple Nucleotide Variation (SNV) dataset from the TCGA, specifically the Level 4 data processed via the “MuTect2” algorithm. This data was retrieved from the Genomic Data Commons (GDC) portal (https://portal.gdc.cancer.gov/).

To examine and illustrate the Single Nucleotide Polymorphisms (SNPs), we utilized the “oncplot” function available in the R package “maftools.” (34), allowing for comprehensive visualization and analysis of these genomic alterations, aiding in the interpretation of their potential biological significance.

To develop a personalized prognostic model, we focused on ligand-receptor (LR) pairs that demonstrated significant relevance for patient outcomes, incorporating them into a penalized Cox regression model using the L1-penalized LASSO (Least Absolute Shrinkage and Selection Operator) methodology, implemented with the R package “glmnet.” The optimal λ value, which controls the strength of the penalty applied to the model, was determined using ten-fold cross-validation, where the data was divided into ten subsets, with each subset used as a validation set once. The optimal λ was chosen based on minimizing the partial likelihood deviance, thereby retaining the most important predictor variables by shrinking the coefficients of less relevant ones to zero, ensuring the best predictive performance while avoiding overfitting. After identifying the predictor variables, we further refined the model using stepwise multivariate regression analysis with the Akaike Information Criterion (AIC) via the “stepAIC” function in the “MASS” package in R. This approach iteratively removed variables to minimize the AIC value, achieving an optimal balance between model complexity and goodness of fit. The final set of 10 LR pairs included in the LR.score were those that remained stable across multiple models and consistently contributed to predicting overall survival. The final risk score for each patient, termed the LR.score, was computed as follows:

To investigate the relationship between drug responsiveness and the LR.score prognostic model, we obtained drug sensitivity data for nearly a thousand cancer cell lines from the Genomics of Drug Sensitivity in Cancer (GDSC) database. We used the Area Under the Curve (AUC) values for various anti-cancer drugs as indicators of drug efficacy. Spearman’s rank correlation analysis was employed to assess the correlation between each drug’s sensitivity and the LR.score. We set a significance threshold at an absolute Spearman’s correlation coefficient (|ρ|) greater than 0.2 and adjusted for False Discovery Rate (FDR) using the Benjamini-Hochberg method with a significance level set at less than 0.05. To further refine drug response predictions, we utilized the “pRRophetic” R package (35). Additionally, for a more targeted approach, we incorporated transcriptomic and clinical data from the IMvigor210 cohort of patients with metastatic bladder cancer who were treated with the anti-PD-L1 drug Atezolizumab. This data was accessed using the “IMvigor210CoreBiologies” R package (36). The immune checkpoints list waw derived from the HisgAtlas database (37).

Bone marrow specimens were collected from 24 individuals with a primary diagnosis of AML at the Zhuji Affiliated Hospital of Wenzhou Medical University. Additionally, 12 healthy bone marrow samples were donated from individuals undergoing total hip arthroplasty, serving as a control group. All participants provided informed consent, and the study was conducted in line with the Declaration of Helsinki principles. The research was approved from the Ethics Committee of Zhuji People’s Hospital Affiliated to Wenzhou Medical University. Diagnosis for AML was made using the French-American-British (FAB) system, alongside tests such as immunophenotyping, cytogenetic analysis, and molecular genetic profiling. A complete response (CR) to treatment was identified by several criteria, including bone marrow blasts below 5%, no Auer rod-positive blasts, no extramedullary leukemia, an absolute neutrophil count above 1.0 × 10^{^9}/L, and a platelet count over 100 × 10^{^9}/L. RNA was extracted from various tissues and cells using Trizol reagent as per the manufacturer’s instructions. The reverse transcription step was carried out using the PrimeScript RT reagent Kit. Quantitative real-time PCR (qRT-PCR) was then performed with SYBR Prime Script RT-PCR Kits, following the prescribed procedure. Expression levels of CLEC11A, ICAM4, ITGA4, and AVP were quantified using the 2^-ΔΔCt method and normalized against GAPDH mRNA. These expression levels were reported relative to a control level set to a baseline of 1.0. The primer sequences used for amplifying CLEC11A, ICAM4, ITGA4, and AVP genes were as follows: for CLEC11A, forward (5′-GGG CCT CTA CCT CTT CGA AA-3′) and reverse (5′-CAG TTC TCG AGC GTG CCA CC-3′); for ICAM4, forward (5’-GCCTACAGTGAGGGACAGG-3’) and reverse (5’-ATCACGGGCTGCCAGAAG-3’); for ITGA4, forward (5′-TTCCAGAGCCAAATCCAAGAGTAA-3′) and reverse (5′-AAGCCAGCCTTCCACATAACAT-3′); and for AVP, forward (5′-GGGCAGGTAGTTCTCCTCCT-3′) and reverse (5′-CACCTCTGCCTGCTACTTCC-3′).

In addition to the mRNA analysis, protein levels of CLEC11A, ICAM4, ITGA4, and AVP were measured using enzyme-linked immunosorbent assays (ELISAs). Bone marrow plasma was separated by centrifugation at 1500 × g for 10 minutes at 4°C from the same 24 AML patients and 12 healthy controls as described above and stored at -80°C until analysis. ELISA kits specific for human CLEC11A (#ELK6393), ICAM4 (#ELK3048), ITGA4 (#ELK9691), and AVP (#ELK5414) were used according to the manufacturer’s protocols. All reagents and samples were brought to room temperature (18-25°C) before use. The 25× Wash Buffer was diluted to 1× with double-distilled water. Standard dilutions were prepared from a stock concentration of 60 ng/mL to create a standard curve with points at 60, 30, 15, 7.5, 3.75, 1.88, and 0.94 ng/mL. 100 μL of each standard, control, or sample was added in duplicate to the appropriate wells of a pre-coated microplate, which was then incubated at 37°C for 80 minutes. After incubation, the wells were washed three times with 200 μL of 1× Wash Buffer, followed by the addition of 100 μL of Biotinylated Antibody Working Solution to each well. The plate was incubated for another 50 minutes at 37°C, and the washing step was repeated. Next, 100 μL of Streptavidin-HRP Working Solution was added to each well, followed by incubation at 37°C for 50 minutes and another washing step repeated five times. The colorimetric reaction was developed by adding 90 μL of TMB Substrate Solution to each well, incubating the plate in the dark at 37°C for 20 minutes. The reaction was stopped by adding 50 μL of Stop Solution, changing the color from blue to yellow. Absorbance was immediately measured at 450 nm using a microplate reader (R&D, MN, USA). The concentration of each protein in the samples was then determined by comparing the corrected absorbance values to the standard curve generated from the known concentrations of the standards.

For the evaluation of statistical differences between two continuous and normally distributed variables, we utilized the unpaired Student’s t-test. In cases where the variables were not normally distributed, the Wilcoxon rank-sum test was employed instead. For comparisons involving three or more groups with non-parametric distributions, Kruskal-Wallis tests were conducted. For examining relationships between categorical variables, Fisher’s exact test was implemented. Correlation analyses were conducted using Spearman’s rank correlation test to assess the strength and direction of associations between variables. All visualizations and statistical computations were carried out using R software, version 4.3.2, provided by the R Foundation for Statistical Computing in Vienna, Austria.

To pinpoint Ligand-Receptor (LR) pairs significantly correlated with the survival outcomes in AML, we incorporated data from three distinct AML cohorts: TCGA-LAML, GSE37642, and GSE12417. Initially, survival analysis for these LR pairs was executed separately for each cohort. This was followed by a meta-analysis where we amalgamated the P-values pertaining to the prognostic importance of these LR pairs from all three cohorts (adjust P-value< 0.01). Upon completing multiple hypothesis testing corrections, we isolated 94 LR pairs that demonstrated substantial prognostic relevance. Of these, 56 LR pairs were indicative of a poor prognosis, while 38 suggested a favorable outcome ( Figure 1B ; Supplementary Table S2 ). The interactive network involving these LR pairs is depicted in ( Figure 1C ).

Additionally, we undertook KEGG pathway enrichment analysis targeting the specific ligands and receptors from these pairs. Remarkably, these 94 LR pairs were predominantly enriched in key biological pathways. These include Cytokine−cytokine receptor interaction, Viral protein interaction with cytokine and cytokine receptor, PI3K−Akt signaling pathway, Proteoglycans in cancer, Cell adhesion molecules (CAMs), and ECM−receptor interaction ( Figure 1D ).

Further, we summed the expression levels of the receptor and ligand genes to represent the expression intensity of each LR pair. We then used the gene expression levels of these LR pairs for molecular subtyping. In this step, we included the 94 LR pairs that were identified in the previous analysis as being significantly correlated and prognostically relevant. Using Consensus Clustering, we clustered 142 AML samples from the TCGA cohort. Based on the Cumulative Distribution Function (CDF), we determined the optimal number of clusters. The CDF Delta area curve suggested that choosing three clusters provided the most stable clustering results (as shown in Figures 2A, B ). Ultimately, we selected k=3 to obtain three molecular subtypes ( Figure 2C ). Upon further analysis of the prognostic features of these three subtypes, we observed significant prognostic differences among them. Specifically, subtype C3 showed better prognosis, while subtype C1 had poorer outcomes (P = 0.00052; Figure 2D ). Additionally, using the same methodology, we conducted molecular subtyping for the AML patient cohort in GSE37642, and observed similar significant prognostic differences among these three subtypes (P< 0.0001; Figure 2E ), consistent with the training set. The same phenomenon was also observed in the GSE12417 cohort (P< 0.0001; Figure 2F ).

Additionally, in the TCGA dataset, we compared the distribution of various clinical pathological features across the three molecular subtypes to examine if these features varied among the subtypes. We found that there were differences in “age” groups, “CALGB Cytogenetics Risk Category,” and “FAB Category” among the three molecular subtypes. Specifically, the C1 subtype had a significantly higher proportion of older patients compared to the C3 subtype (Fisher’s exact test, -log10 P-value = 4.74). Most patients in the C3 subtype fell into the “Favorable” (Fisher’s exact test, -log10 P-value = 7.53) ( Figure 3A ). Similarly, we examined differences in clinical information among different molecular subtypes in the GSE37642 and GSE12417 cohorts ( Figures 3B, C ). We observed that the majority of patients did not have RUNX1 mutations. These analyses further support the notion that these molecular subtypes not only have prognostic value but also correlate with various clinical and pathological features, thereby potentially aiding in more targeted treatment approaches.

We further investigated the differences in genomic alterations among the three molecular subtypes within the TCGA cohort. Our analysis revealed that there were no significant differences among the molecular subtypes in terms of Aneuploidy Score, Homologous Recombination Defects, Fraction Altered, Number of Segments, and Tumor Mutation Burden ( Figure 4A ).

Moreover, we scrutinized the variations in gene mutations across the distinct molecular subtypes. We highlighted the top 20 genes that showed substantial differences in their mutation rates ( Figure 4B ). Of particular interest were genes like DNMT3A, NPM1, and RUNX1, which exhibited noticeable disparities in mutation frequency among the three molecular subgroups. Intriguingly, a higher prevalence of RUNX1 and DNMT3A mutations was observed in subtypes C1 and C2, and these mutations were correlated with unfavorable prognostic outcomes. This observation aligns with existing scientific literature that associates these specific mutations with poorer survival rates in AML (38, 39).

These findings imply that while some genomic features like Aneuploidy Score and Tumor Mutation Burden may not differ significantly among the subtypes, specific genes do show variations in mutation rates. This could have implications for understanding the biological distinctions among the subtypes and potentially for targeted therapeutic strategies.

Next, we investigated whether there were differentially activated pathways within the various molecular subgroups. To identify these pathways, we performed Gene Set Enrichment Analysis (GSEA) using all candidate gene sets from the Hallmark database, with FDR threshold less than 0.05 for significant enrichment. In the TCGA cohort, when compared to the C3 subtype, the C1 subtype exhibited activation of 13 pathways and suppression of 3 pathways, as shown in ( Figure 5A ). These pathways were mainly related to immune responses, such as INTERFERON_GAMMA_RESPONSE, INFLAMMATORY_RESPONSE, INTERFERON_ALPHA_RESPONSE, ALLOGRAFT_REJECTION, and COMPLEMENT. In addition, we analyzed the significant enriched gene sets in the two validation cohorts when comparing the C1 subtype with the C3 subtype. Generally, activated pathways primarily included immune markers, such as INTERFERON_GAMMA_RESPONSE, IL6_JAK_STAT3_SIGNALING, IL2_STAT5_SIGNALING, and INFLAMMATORY_RESPONSE ( Figure 5B ). Furthermore, we compared differential pathways between C1 and C2, C1 and C3, as well as C2 and C3 in TCGA cohort. As shown in ( Figure 5C ),

we identified a notable concentration of immune-related pathways among the diverse subtypes. GSEA analysis across these subtypes revealed that individuals in the C1 subtype typically displayed heightened activity in immune regulatory pathways. This leads us to hypothesize that the ligand-receptor pairs utilized for molecular classification might have significant influence in shaping the immune microenvironment.

To delve deeper into the variations in the immune microenvironment across patients from different molecular groups, we assessed the relative proportions of 22 immune cell types across the three AML cohorts using the CIBERSORT algorithm. This revealed pronounced differences in the presence of specific immune cells across the three molecular subtypes (as illustrated in ( Figures 6A, C, E ). Notably, the C3 subtype exhibited elevated levels of CD8 T cells and plasma cells, suggesting heightened adaptive immune responses. Additionally, we employed the ESTIMATE algorithm to assess immune cell infiltration, as shown in ( Figures 6B, D, F ). We found that the “ImmuneScore” was consistently highest in the C2 subtype across the TCGA, GSE37642, and GSE12417 cohorts, indicating a relatively higher level of immune cell infiltration in the C2 subtype, importantly, C2 also has highest myeloid lineage cells compare to C1 and C2.

We found that molecular subtypes based on LR pairs exhibit different mutation landscapes, distinct pathway features, and varying levels of immune infiltration. To further refine our risk model, we employed lasso regression on the 94 significant LR pairs identified from the meta-analysis within the TCGA-LAML cohort. We showed that as the value of lambda increases, the number of coefficients tending towards zero also increases ( Supplementary Figure S1A ). Using 10-fold cross-validation, we identified an optimal lambda value of 0.0868 and selected 13 LR pairs for further analysis ( Supplementary Figure S1B ). Subsequently, we used stepwise multivariate regression analysis and employed the Akaike Information Criterion (AIC) to further refine our model, ultimately identifying 10 key LR pairs. These pairs include “AVP->AVPR1B,” “CALCA->CALCRL,” “CCL7->ACKR4,” “CLEC11A->ITGA11,” “CXCL12->ITGA4,” “HGF->MET,” “ICAM4->ITGA4,” “IL2->IL2RA,” “NCAM1->ROBO3,” and “NLGN1->NRXN2” ( Supplmentary Figure S1C ).

Building on the 10 distinguished ligand-receptor (LR) pairs, we developed an LR scoring system, termed LR.score, to quantifiably gauge the activity patterns of these LR pairs in AML patients. We observed a notably elevated LR.score in patients of the “C1” subtype in comparison to their “C3” counterparts ( Figure 7A ). Probing the clinical relevance of the LR.score, we categorized patients into high and low LR.score groups, using “0” as the demarcation point. Interestingly, those with a diminished LR.score exhibited a pronounced survival advantage ( Figure 7B ; log-rank test, P< 0.0001). The receiver operating characteristic (ROC) curve’s area under the curve (AUC) for 1-, 3-, and 5-year overall survival stood at 0.85, 0.84, and 0.87, respectively ( Figure 7C ). In our validation datasets, GSE37642 and GSE12417, a similar pattern emerged where the LR.score for “C1” was appreciably higher than for “C2” and “C3” (Wilcoxon rank-sum test, P< 0.001) ( Figures 7D, G ). In alignment with our initial observations, patients with a reduced LR.score in these validation sets also showcased a marked survival benefit ( Figures 7E, H ; log-rank test, P< 0.0001). The AUC values from the ROC assessments were 0.67, 0.71, and 0.69 for 1-, 3-, and 5-year overall survival in GSE37642, and 0.73, 0.74, and 0.72 in GSE12417, respectively. The consistency in AUC values across these intervals underscores the reliability of the model’s prognostic potential ( Figures 7F, I ).

To examine whether the LR.score could serve as an independent prognostic factor, we performed both univariate and multivariate Cox regression analyses using patient clinical features such as age, gender, cytogenetics risk category, and FAB category. We found that the LR.score is a reliable and independent prognostic biomarker ( Figures 7J, K ; HR=3.26, 95% CI 2.29-4.65, P = 5.68E-11).These results suggest that the LR.score can reflect the LR-pairs patterns in AML patients and predict prognosis.

To investigate the relationship between the LR.score and clinical characteristics of AML, we analyzed the differences in the LR.score based on various clinical-pathological features in the TCG dataset. Our findings indicated that as age increases, the LR.score also rises. Moreover, higher “cytogenetics risk” levels were associated with elevated LR.score values ( Figure 8A ). We further compared the relationship between patients’ LR.score and their clinical-pathological characteristics in the GSE37642 and GSE12417 cohorts ( Figures 8B, C ). It was observed that older patients tended to have higher LR.score as well as RUNX1 mutation patients. In summary, the higher the clinical staging level of the patient, the greater the LR.score.

We further delved into the distribution differences of scores for 22 immune cell types across LR.score groupings in the TCGA cohort, as depicted in ( Figure 9A ). Generally, scores for most immune cells do not display significant disparities between LR.score groups. However, certain cells like T_cells_gamma_delta, NK_cells_activated, and Mast_cells_resting exhibited notable differences. Moreover, when comparing immune infiltration, the ImmuneScore was consistently lower in the low LR.score group than in the high LR.score group ( Figure 9B ). Additionally, Employing Pearson’s correlation coefficient, we evaluated the association between the LR.score and immune cell infiltration, the results of which are shown in ( Figure 9C ). Interestingly, a strong positive correlation was observed between LR.score and T_cells_regulatory (Tregs), while a significant negative correlation was noted with Mast_cells_resting.

To understand the relationship between LR.score and biological functionality, we employed single-sample GSEA analysis (ssGSEA) on the gene expression profiles of AML samples from the TCGA cohort. This enabled us to compute scores for each sample across various functionalities, leading to ssGSEA scores for each function across individual samples. Further correlation analysis between these functional scores and LR.score revealed functions with correlations greater than 0.35, as shown in ( Figure 9D ) and ( Supplementary Table S3 ). As a result, 17 pathways positively correlated with the LR.score of samples, while 2 pathways demonstrated a negative correlation. Notably, cancer-related pathways such as KEGG_GLYCEROPHOSPHOLIPID_METABOLISM, KEGG_GLYCEROLIPID_METABOLISM, KEGG_REGULATION_OF_ACTIN_CYTOSKELETON, and KEGG_MAPK_SIGNALING_PATHWAY were positively correlated with LR.score.

Furthermore, we investigated if there were disparities between LR.score groupings concerning their responses to treatment. Initially, we compared the expression of immune checkpoints between the LR.score groups. As shown in ( Figure 10A ), certain immune checkpoint genes exhibited differential expression across the LR.score groups. Specifically, high gene expressions of ARHGEF5, CD274, CD80, CTLA4, LAG3, PDCD1, and VISTA are associated with high LR.score (P< 0.05, Wilcoxon rank-sum test).

We subsequently examined the disparities in immunotherapy outcomes between the LR.score categories. Utilizing the TIDE algorithm, we gauged the probable clinical ramifications of immunotherapy among our designated high and low LR.score groups. Notably, a heightened TIDE prediction score implies an augmented propensity for immune evasion, hinting at potentially diminished benefits from immunotherapy for the patient. As illustrated in ( Figure 10B ) and (Supplementary Table S4 ), discernible distinctions between high- and low- LR.scores in MDSC, CAF, and TIDE scores were absent. Simultaneously, we contrasted the anticipated T-cell dysfunction and exclusion scores across various metabolic molecular subtypes within the TCGA cohort. The group with elevated LR.score exhibited the most pronounced T-cell dysfunction score (P< 0.01, Wilcoxon rank-sum test).

Delving deeper into the ramifications of LR.score on drug responsiveness, we explored its association with reactions in tumor cell lines to drugs. Through Spearman correlation analysis, we discerned four notable associations between LR.score and drug susceptibility in the Genomics of Drug Sensitivity in Cancer (GDSC) database ( Figure 10C ). Among these, three correlations indicated drug resistance in tandem with the LR.score, encompassing Forentinib, BPD-00008900, and Vinblastine. We also gauged the variances in responses to prevalent chemotherapy agents such as ‘Bexarotene’, ‘Bortezomib’, ‘Erlotinib’, and ‘Rapamycin’ between the LR.score categories ( Figure 10D ). Observations revealed that patients with a diminished LR.score manifested heightened sensitivity to drugs, notably Bortezomib (P = 0.0039, Wilcoxon rank-sum test) and Erlotinib (P = 0.028, Wilcoxon rank-sum test), as compared to their high LR.score counterparts. Collectively, these insights underscore a linkage between LR pairs and drug sensitivity, positing the LR.score as a potential biomarker to inform tailored therapeutic approaches.

Considering the LR.score model’s proven capability in predicting the risk and treatment outcomes for AML, we aim to validate the individual gene expression profiles to determine whether the model possesses the requisite sensitivity and specificity to be applicable for validation within our patient cohort. First of all, we investigated the cancer specificity of certain gene-protein pairs by checking the cell line database in The Human Protein Atlas (HPA), focusing on the following associations: “AVP->AVPR1B,” “CALCA->CALCRL,” “CCL7->ACKR4,” “CLEC11A->ITGA11,” “CXCL12->ITGA4,” “HGF->MET,” “ICAM4->ITGA4,” “IL2->IL2RA,” “NCAM1->ROBO3,” and “NLGN1->NRXN2”. Our analysis revealed that, among the genes studied, CLEC11A, ICAM4, ITGA4, and AVP exhibited the highest specificity to AML when compared to other cancer types. This suggests a strong association between the expression levels of these genes and AML, pointing to their potential role as specific biomarkers for this disease ( Figure 11A ). Next, we examined the expression of CLEC11A, ICAM4, ITGA4, and AVP in 12 normal and 24 AML bone marrow samples by quantitative reverse-transcription polymerase chain reaction (qRT-PCR). In comparison to the healthy bone marrow ( Figure 11B ), the levels of CLEC11A, ITGA4, and AVP were notably reduced in samples from AML patients, with a highly significant statistical difference (P< 0.0001). The expression of ICAM4 was also found to be slightly diminished (P< 0.05). These findings indicate that CLEC11A, ICAM4, ITGA4, and AVP may act as inverse biomarkers in the context of AML. Additionally, we used enzyme-linked immunosorbent assay (ELISA) to further validation these genes in the protein level. As shown in Figure 11C , the protein concentrations of CLEC11A and ITGA4 were significantly higher in normal samples compared to AML (P< 0.0001), and the same trend was observed for ICAM4 and AVP (P< 0.01). Furthermore, survival analysis reveals that patients with higher levels of CLEC11A has a better prognosis than those with lower expression ( Figure 11D , P< 0.0001), reinforcing the potential role of CLEC11A as a particularly significant negative biomarker.