Integrative analyses of single-cell and bulk RNA sequencing to construct the tumor-associated macrophage-related prognostic signature in lung adenocarcinoma

Data collection and processing

The RNA-seq data and miRNA-seq data of LUAD, along with the corresponding clinical data, were downloaded from The Cancer Genome Atlas (TCGA, https://www.cancer.gov/ccg/research/genome-sequencing/tcga). After excluding the sample with incomplete follow-up information, a total of 486 LUAD and 59 normal samples from the RNA-seq data, as well as a total of 521 LUAD and 46 normal samples the miRNA-seq data, were included for subsequent analysis in this study. The mRNAs and lncRNAs of samples were annotated based on the gene transfer format (GTF) files containing the gene symbol. Additionally, transcriptomic profiles and corresponding clinical information in the GSE31210 dataset, which includes 226 tumor samples, were obtained from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) and generated by GPL570 platforms. Moreover, scRNA-seq files from nine LUAD tumor specimens were obtained from the GEO database (accession number GSE171145), generated by the GPL24676 platform (Illumina NovaSeq 6000). Finally, a total of 428 tumor-associated macrophage (TAM)-related genes with a score greater than five were obtained from the GeneCards database using keywords searches (Table S1).

Screening differentially expressed RNAs in LUAD

Limma R package was performed to identify the differentially expressed lncRNAs (DE-lncRNAs), mRNAs (DEGs), and miRNAs (DE-miRNAs) according to the criteria of |log2 (fold change, FC)|> 1 and adjusted P value < 0.05. Volcano plots were drawn using the ggplot2 R package to visualize the differentially expressed RNAs, and a heatmap was drawn using the pheatmap R package to illustrate the top50 DE-lncRNAs, top50 DEGs, and all DE-miRNAs.

Construction and validation of a TAM-related gene signature

The differentially expressed TAMGs (DE-TAMGs) were identified by overlapping the DEGs and 428 TAMGs. The 486 patients from the TCGA-LUAD cohort were randomly split into training and testing sets at a 7:3 ratio using the caret R package. Afterward, the DE-TAMGs that were associated with survival were identified using univariate Cox analysis, these DE-TAMGs were incorporated into a least absolute shrinkage and selection operator (LASSO) regression model to select the prognostic gene signature by glmnet R package. A risk score was calculated as follows, risk score = ${\sum }_1^i\left(CoefGenei\ast ExpGenei\right)$. Coef represents the regression coefficient, and ExpGene represents the expression values of the gene. All patients were stratified into high- and low-risk groups based on the median risk score. The correlation between risk score and clinical characteristics (age, gender, pathologic stage, and T/N/M stages) were analyzed using Pearson correlation analysis. The Kaplan–Meier method and log-rank test were applied to compare the overall survival (OS) of patients between high- and low-risk groups using the survival R package. The time-dependent receiver operator characteristic (ROC) curve was conducted to evaluate the predictive accuracy of the risk model using the survivalROC R package.

Gene Set Enrichment Analysis (GSEA)

The pathway enrichment between high- and low-risk groups was analyzed using GSEA. The c2.cp.kegg.v7.4.symbols.gmt and h.all.v7.4.symbols.gmt were selected from the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb) for the enrichment analysis, which was performed using the GSEA R package.

Tumor immune microenvironment analysis

The ESTIMATE algorithm was used to evaluate the immune score, stromal score, and ESTIMATE score using the estimate R package. Differences in immune score, stromal score, and ESTIMATE score between high- and low-risk groups were detected using the Wilcoxon-test. Besides, the CIBERSORT algorithm was performed to evaluate the abundance of 22 types of immune cells based on the gene signatures consisting of 547 genes. The significant differences in infiltrated immune cells between high- and low-risk groups were detected using the Wilcoxon test.

Immunotherapeutic and chemotherapeutic effect analysis

Human leukocyte antigen (HLA) expression plays a key role in tumor immunogenicity and responses to immunotherapies (Gong & Karchin, 2022). Thus, the differences in HLA genes and immune checkpoint levels between high- and low-risk groups were detected using Wilcoxon test. The Tumor Immune Dysfunction and Exclusion (TIDE) (Lu et al., 2019), Immunophenoscore (IPS) (Charoentong et al., 2017), and the submap algorithm were used to predict the response to immune checkpoint blockade. A P value < 0.05 was considered statistically significant. The pRRophetic algorithm was conducted to estimate the therapeutic response of the samples with different risk scores based on the Genomics of Drug Sensitivity in Cancer (GDSC).

Construction of a lncRNA-miRNA-mRNA network

The miRNA-mRNA pairs were identified based on risk signatures and DE-miRNAs in the online databases, including miRcode (https://bio.tools/miRcode), miRDB (https://mirdb.org/), miRTarBase (https://mirtarbase.cuhk.edu.cn/), Starbase (https://starbase.sysu.edu.cn/), and TargetScan (https://www.targetscan.org/). Additionally, miRNA-lncRNA pairs were screened out based on DE-lncRNA in Starbase (https://starbase.sysu.edu.cn/). Finally, an endogenous competitive network, lncRNA-miRNA-mRNA, was constructed based on the DE-lncRNAs, shared-miRNA that lncRNA-miRNA and miRNA-mRNA pairs, and risk gene signature. The network was visualized using Cytoscape Version 3.8.0.

Clinical specimen

Tumor tissues and paracancerous non-tumor tissue samples were collected from lung cancer patients (n = 28) at Yueqing People’s Hospital. Histological examination was used to confirm the diagnosis of LUAD. All procedures were approved by the Ethical Committee of Yueqing People’s Hospital (YQYY202300221). Written informed consent was obtained from all participants prior to their inclusion in the study. The tissue specimens were immediately frozen in liquid nitrogen after collection and stored at −80 °C for further RNA extraction.

Experimental validation of the TAM-related genes involved in ceRNAs

Total RNA was extracted and purified using TRIzol (TAKARA, Dalian, China) according to the manufacturer’s instructions. cDNA synthesis was performed using a Bestar qPCR RT Kit (DBI Bioscience, Shanghai, China) according to the manufacturer’s protocol. Data were collected as previously described in He et al. (2016), Specifically the experimental validation of the TAM-related genes involved in ceRNAs. Quantitative real-time PCR (RT-qPCR) analyses were performed with a Bestar^®SybGreen qPCR Mastermix (DBI Bioscience, Shanghai, China) following: 95 °C for 20 s, followed by 40 cycles of 95 °C for 1 s and 60 °C for 20 s. GAPDH was used as the reference gene for the normalization of all gene expression results. The average of three independent analyses for each gene was calculated. The fold changes were calculated through relative quantification (2^−ΔΔCt). All reactions were run in triplicate and repeated in three independent experiments. The primers used were as follows: GAPDH forward, 5′-TGTTCGTCATGGGTGTGAAC-3′ and reverse 5′-ATGGCATGGACTGTGGTCAT-3′. PDGFB forward, 5′-TTATCATGGGCCTCGGGGA-3′ and reverse 5′-CAGACGGACGAGGGAAACAA-3′. CD74 forward, 5′-GGCTACTGCTGGTGTGTCTT and reverse,5′-TCCAAGGGTGACGAAAGAGC-3′. KLF4 forward, 5′-GTCCCGGGGATTTGTAGCTC-3′ and reverse 5′-TGTAGTGCTTTCTGGCTGGG-3′. CCL20 forward, 5′- TGTCAGTGCTGCTACTCCAC-3′ and reverse 5′-ACAAGTCCAGTGAGGCACAA-3′.

Data processing and single cell analysis

The Serut R package was used for standard downstream processing for scRNA-seq data. Gene expression was required in at least three cells, with a gene count more than 300 and less than 7,000, the mitochondria proportion less than 10%, and the erythrocyte proportion less than 3%. Cells that did not meet these criteria were excluded from the analysis. Afterward, “NormalizeData” was used to perform the data normalization. PCA (Principal Component Analysis) was used to reduce the dimensionality of the data, after which gene expression value were converted to z-score using “ScaleData” function. Then, Uniform Manifold Approximation and Projection for Dimension Reduction (UMAP) was utilized for unsupervised clustering and unbiased visualization of cell populations on a two-dimensional map. To account for cell cycle differences, “CellCycleScoring” was used to detect cell cycle-related genes. “FindAllMarkers” was performed to identify marker genes of each cluster with the threshold value of |log (fold change, FC)| > 0.025 and a minimum cell population fraction of 0.25 in either of the two populations. The SingleR package was utilized for cell-type annotation. The classical gene markers of LUAD were visualized as bubble plots by ‘DotPlot’. The hub gene expression across different cell types, along with significant cell types, was selected for biological function (GO and KEGG pathway enrichment) analysis using clusterProfiler R package. Monocle2 R package was used for cell trajectory and pseudo-time analysis with the DDRTree method.

Identification of the differential expression of RNAs in LUAD

The flowchart of this study are shown in Fig. S1. In the present study, the limma R package was performed to identify the differential expression of RNAs in LUAD, resulting in a total of 316 DE-lncRNAs (116 upregulated and 200 downregulated DE-lncRNAs, Figs. 1A, 1D), 162 DE-miRNAs (125 upregulated and 37 downregulated DE-miRNAs, Figs. 1B, 1E), and 2,601 DEGs (924 upregulated and 1,677 downregulated DEGs, Figs. 1C, 1F) were identified in LUAD (Tables S2–S4).

Construction of a TAM-related gene signature in LUAD

A total of 147 DE-TAMGs were identified from the 2,601 DEGs (Fig. 2A) and subsequently incorporated into a univariate Cox regression model, which revealed 12 DE-TAMGs with significant associations (p-value < 0.05) in the training set, including CAT, KLF4, CLEC12A, GAPDH, PDGFB, HSPD1, TIMP1, ITGAL, CA9, TNFSF11, CD74, CX3CL1, FOLR1, CLEC10A, and CCL20 (Fig. 2B). LASSO regression analysis further narrowed down the candidate genes to KLF4, GAPDH PDGFB, TIMP1, CD74, and CCL20, which were ultimately selected to form the prognostic signature for LUAD (Figs. 2C–2D, Table S5).

The risk score for each patient was calculated based on the expression levels and corresponding regression coefficients of these six genes. Patients in the training set were subsequently divided into high-risk (n = 179) and low-risk (n = 163) groups (Table S6). We found that risk scores were associated with the ages and gender of patients (Figs. 2E–2G). Additionally, a higher risk score correlated with tumorigenesis, lymph node metastasis, and distant metastasis, and more aggressive progression (Figs. 2E, 2H–2K).

Validation of the TAM-related gene signature in LUAD

To validate the prognostic potential of the TAMG-related risk score, we applied it to several datasets, including the training TCGA-LUAD set, the test TCGA-LUAD set, the entire TCGA-LUAD set, and the external GSE31210 set. Based on the median risk score, patients were categorized into high- and low-risk groups. Patients in the high-risk group showed worse survival status compared to those in the low-risk group (Figs. 3A–3D). Kaplan–Meier OS survival curves indicated significantly poorer prognoses for patients in the high-risk group compared to those in the low-risk group (Figs. 3E–3H). Time-dependent ROC curves for those sets at 1-, 2-, 3-, 4-, and 5 years showed the risk model and TAM-related gene signature could predict the OS of LUAD patients with moderate sensitivity and specificity (Figs. 3I–3L).

TAM-related functional enrichment analysis

Based on previous analyses, we further conducted the enrichment analysis using GSEA to explore the functional pathways associated with the TAM-related risk score. The results indicated that several canonical cancer pathways were significantly enriched between the high- and low-risk groups (Figs. 4A–4B). Specifically, we found that the cell cycle, NOTCH signaling pathway, p53 signaling pathway, pathways in cancer, proteasome, hypoxia, PI3K/AKT/mTOR signaling pathway, reactive oxygen species pathway, and Wnt/beta-catenin signaling pathway were significantly enriched in the high-risk group compared to the low-risk group (Figs. 4A–4B).

Correlation analysis between TAM risk score and tumor immune microenvironment

We further investigated the immune characteristics within the tumor microenvironment and observed that the low-risk group had higher immune and ESTIMATE scores compared to the high-risk group (Figs. 5A–5C). In addition, we found significant differences s in the proportions of various immune cell types (Fig. 5D). In the high-risk group, the proportions of macrophages M0, T cells CD4 memory activated, neutrophils, NK cells resting, Mast cells activated was upregulated, but the fraction of T cells CD4 memory resting, mast cells resting, dendritic cells resting, T cells regulatory (Tregs), monocytes, and B cells memory were downregulated compared to the low-risk group (Fig. 5E). Macrophages M0 are unactivated macrophages, typically considered to be in a “resting” or undifferentiated state. They perform various immune surveillance and regulatory functions but usually require specific immune signals to differentiate into more immunologically active M1 or M2 macrophages. In the analysis of TAM-related high- and low-risk groups, M0 macrophages were significantly increased in the high-risk group, indicating a close association between TAMs and high-risk status.

Correlation analysis between TAM risk score and immunotherapeutic effect

We next compared the expression of HLA-related genes and immune checkpoint markers between high- and low-risk groups. Most HLA-related genes, including HLA class I (HLA-A/B/C/E/F) and HLA class II (HLA-DMA/B, HLA-DOA/B, HLA-DPA1/B1, HLA-DPQ1/B1, HLA-DRA1/B1, and HLA-DRB5), were downregulated in high-risk group than low-risk group (Fig. 6A). Similarly, the expression of immune checkpoint genes, such as CD80, CTLA4, and TIGIT was downregulated in a high-risk group than the low-risk group (Fig. 6B).

As expected, the TIDE score was significantly higher in a high-risk group than in a low-risk group (Fig. 6C), with a strongly positive correlation between risk score and TIDE score (Fig. 6D). In contrast, the IPS score was lower in the high-risk group compared to the low-risk group (Fig. 6E). These findings suggest that the low-risk group exhibits more active immunological functions and greater sensitivity to immunotherapy. However, the patients in the low-risk group did not show significant responses to anti-PD1 and anti-CTLA4 therapies (Fig. 6F).

Correlation analysis between TAM risk score and chemotherapeutic effect

We also identified the potential chemotherapeutic agents targeting the TAM-related prognostic signature. A total of 13 compounds showed significantly different half-inhibitory concentrations (IC50) between the high- and low-risk groups (Table S7). Most of these compounds exhibited higher IC50 values in the high-risk group than low-risk group (Fig. 7).

Development of a TAM-associated lncRNA-miRNA-mRNA network in LUAD

Using several online databases, including miRcode, miRDB, miRTarBase, Starbase, and TargetScan, a total of 66,130 miRNA-mRNA pairs were identified (Table S8), From these, 24 miRNA-TAMG (PDGFB, CD74, KLF4, and CCL20) pairs were selected for further analysis (Table S9). In addition, a total of 7,975 miRNA-lncRNA pairs were screened (Table S10).

Subsequently, a lncRNA-miRNA-mRNA network was constructed, incorporating 36 DE-lncRNAs, 23 shared-miRNA, and four TAMGs (PDGFB, CD74, KLF4, and CCL20, Fig. 8A). QPCR results confirmed that CCL20 expression was significantly upregulated in tumor tissues (Fig. 8E), while the expression of PDGFB, CD74, and KLF4 was downregulated in tumor tissues compared to adjacent normal tissues (Figs. 8B–8D). Additionally, a similar trend was observed in the TCGA-LUAD cohort, where CCL20 expression was significantly upregulated in tumor tissues, whereas PDGFB, CD74, and KLF4 were significantly downregulated compared to adjacent normal tissues (Figs. 8F–8I).

Validation of the key TAMG expression in monocyte/macrophage at the single-cell level

After the quality control of scRNA-seq data, a total of 30,874 cells were selected for subsequent analysis (Figs. S2A–S2B). After normalization and dimensionality reduction, 27 subpopulations were obtained (Figs. S2C–S2I). Then, the relative expression of marker genes in each cluster was presented in the heatmap (Fig. S2J)). Afterward, 17 cell types were annotated using a singleR package (Fig. 9A), including CD4 T cells, CD8 T cells, monocytes, B cells, secretory club cells (Club), cancer cells, macrophages, alveolar cells, NK cells, plasma, mast cells, smooth muscle cell, cycling cells, neutrophils, fibroblasts, endothelial cells, and ciliated airway epithelial cells (Ciliated).

The marker genes for each cell type, along with the top three markers, are shown in Figs. 9B–9C. Furthermore, we examined the expression of PDGFB, CD74, KLF4, and CCL20 in different cell types, with a focus on monocytes/macrophages. Our analysis indicated that CD74, KLF4, and CCL20 were significantly expressed in these cells, whereas PDGFB was downregulated (Fig. 9D, Figs. S3A–S3C).

GO enrichment analysis of macrophage marker genes indicated that significant enrichment in protein-macromolecule adaptor activity, nucleoside-triphosphatase regulator activity, GTPase regulator activity, and chaperone binding (Fig. 10A). For monocytes, marker genes were enriched in kinase regulator activity, DNA-binding transcription activator activity, and RNA polymerase II-specific functions (Fig. 10A).

KEGG pathway enrichment analysis indicated that the macrophage marker genes were mainly associated with calcium reabsorption regulation, galactose metabolism, and the IL-17 signaling pathway (Fig. 10B). In contrast, monocyte marker genes were significantly associated with serotonergic synapses, the Apelin signaling pathway, and non-small cell lung cancer (Fig. 10B).

Finally, cell trajectory and pseudo-time analysis of monocyte and macrophages using the monocle2 R package revealed a clear transition from monocytes to macrophages (Figs. 10C–10F). These results support the notion that tissue-resident macrophages originate from hematopoietic stem cells and monocyte-derived macrophages (Lazarov et al., 2023).