{'fraction_non_alphanumeric': 0.037174927761724826, 'fraction_numerical': 0.028395198933096243, 'mean_word_length': 5.534858387799564, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Significance of macrophage infiltration in the prognosis of lung adenocarcinoma patients evaluated by scRNA and bulkRNA analysis

Purpose To investigate the significance of macrophage infiltration to the prognosis of lung adenocarcinoma. Methods R language bioinformatics analysis technology, was used to obtain macrophage infiltration-related module genes through WGCNA (Weighted Gene Co-Expression Network Analysis). Marker genes of macrophage subtypes were identified using single-cell sequencing of lung adenocarcinoma tissue. Risk score models were constructed and validated using external data cohorts and clinical samples. Results Analysis of cohorts TCGA-LUAD, GSE11969, GSE31210, GSE50081, GSE72094 and GSE8894, revealed a negative correlation between macrophage infiltration and survival. Immunohistochemical analyses of clinical samples were consistent with these data. Based on cell-cluster-markers and TAMs-related-genes, TOP8 genes were obtained (C1QTNF6, CCNB1, FSCN1, HMMR, KPNA2, PRC1, RRM2, and TK1) with a significant association to prognosis. Risk score models including 9 factors (C1QTNF6, FSCN1, KPNA2, GLI2, TYMS, BIRC3, RBBP7, KRT8, GPR65) for prognosis were constructed. The efficacy, stability and generalizability of the risk score models were validated using multiple data cohorts (GSE19188, GSE26939, GSE31210, GSE50081, GSE42127, and GSE72094). Conclusions Macrophage infiltration negatively correlates with prognosis in patients with lung adenocarcinoma. Based on cell-cluster-markers and TAMs-related-genes, both TOP8 genes (C1QTNF6, CCNB1, FSCN1, HMMR, KPNA2, PRC1, RRM2, TK1) and risk score models using C1QTNF6, FSCN1, KPNA2, GLI2, TYMS, BIRC3, RBBP7, KRT8, GPR65 could predict disease prognosis.


Introduction
Lung cancer remains the most common malignancy worldwide and a leading cause of cancer-related death, despite advances in screening and treatment (1,2). Whether it was for the non-small cell lung cancer (NSCLC) or small cell lung cancer (SCLC) patients, immunotherapy was the most shining one among many treatment methods, which had changed the landscape of anti-tumor therapy and brought anti-tumor therapy into a new era (3)(4)(5)(6)(7). However, there were still many details in the screening of immunotherapy benefit populations and related predictors needed to be further elucidated (8- 16). Specific macrophage phenotypes can act as indicators of lung cancer prognosis and the efficacy of immunotherapy (17)(18)(19)(20)(21)(22)(23)(24). Sequencing technologies and R language based bioinformatics, formerly reported (25-27), can be used for studies in this area (28-30). Based on our previous studies (28-30), we performed bioinformatics analysis and clinical sample validation to identify specific macrophage signatures that can act as indicators of therapeutic efficacy.

Immune infiltration analysis
Immune infiltration for each sample was calculated using IOBR of the R package for the training set TCGA expression matrix and GEO data, respectively (method = 'cibersort').

Survival analysis
For survival assessments, R packages "survminer" and "survival" were analyzed and survival curves were constructed based on survival time and status. Differences in prognosis among the groups were assessed.

Screening of modules corresponding to macrophages using WGCNA
Hierarchical clustering analysis was performed on the TCGA expression matrix using the R package "hclust", "method=average". Phenotypic information was obtained using the infiltration ratio of macrophages. A correlation between different modules and macrophages was obtained.

Clustering analysis of samples
The R package "ConsesusClusterPlus" was used to perform consensus clustering analysis. After clustering on the TCGA and GEO data, the optimal number of categories were determined according to the change of area under the CDF curve. The k value of the cluster category ranged from 2 to 6.

Analysis of single-cell data
Single-cell data were filtered using the R package "seurat" to remove cells with ≥ 20% mitochondrial expression. Data were analyzed using the "seurat" normalization pipeline. To identify tumor-associated macrophage (TAM) populations, marker genes from published studies were used to identify corresponding clusters. TAM populations were selected for standardization analysis using "Seurat".

Trajectory analysis of single-cell data
Trajectory analysis was performed on TAM subclasses using the R package "monocle" with default parameters. This resulted in differentiation trajectories and key genes determining these trajectories.

Gene set variation analysis
To investigate differences in the expression patterns of specific TAM isoforms in biological processes, GSVA enrichment analysis was performed using the R package "GSVA". GSVA is a nonparametric, unsupervised method primarily used to assess alterations in signaling pathways and biological processes in samples.

Construction of risk scoring model
Univariate cox regression analysis was performed on "cellcluster-markers" and "TAMs-related-genes", and genes significantly associated with OS survival were screened at the p<0.05 level. According to the identified prognosis-related genes, the R package 'glmnet' was used to construct a prognosis model (or classifier model) with a 10-fold cross-validation fold using the cox method. Characteristic factors were then screened. Kaplan-Meier survival analysis and ROC curves were used to evaluate the predictive power of the prognostic model.

Clinical sample validation (sample collection and immunohistochemistry)
Lung Cancer samples were collected from the First Affiliated Hospital of Shandong First Medical University & Shandong Provincial Qianfoshan Hospital from June 2012 to February 2020. Written informed consent was provided by all participants. Tumor tissues were surgically resected, formalin fixed and paraffin embedded (FFPE) for histological evaluation. HE-stained and immunohistochemical (IHC) slides were examined by two independent and experienced pathologists according to the WHO criteria.
Samples were IHC stained with mouse anti-human CD68 monoclonal antibodies (MAB-0863, clone MX075) and mouse anti-human CD163 monoclonal antibodies (MAB-0869, clone MX081). CD68 was used as a general surface marker for macrophages, whilst CD163 was used as a marker for M2 macrophages (31). Double-labeled immunohistochemical staining was performed using alkaline phosphatase and horseradish peroxidase conjugated secondary antibodies. Substrates were fast red (AP-Red) and diaminobenzidine (DAB) (Roche Ltd) stained. Slides were processed using an automated Roche BenchMark XT staining system according to the manufacturer's protocol.

Proportion of immune infiltrating cells and the prognostic efficacy of macrophages
CIBERSORT was used to evaluate the levels of immuneinfiltration from different lung adenocarcinoma datasets (TCGA-LUAD, GSE11969, GSE31210, GSE50081, GSE72094, and GSE8894). According to the median macrophage ratio, samples were divided into high-and low-levels of macrophage infiltration. Survival differences between high-and low-groups showed a significant correlation with macrophage infiltration (Figure 1; Supplementary

Screening of modules corresponding to macrophages
To identify macrophage-related genes related to infiltration, WGCNA module analysis was performed on the training dataset Genes corresponding to red modules were named "TAMs-relatedgenes" for subsequent analysis.

TAMs-related-gene-based clustering analysis, molecular typing and prognostic assessments
We analyzed the expression profiles of TAMs-related-genes in samples from different lung adenocarcinoma datasets (TCGA-LUAD, GSE13213, GSE31210, GSE72094, and GSE8894) to construct consistent clustering profiles. Based on cumulative distribution functions and incremental area maps, we selected stable clusters of TAMs-related-genes to obtain multiple subtypes ( Supplementary

Preprocessing of single-cell data
To further investigate the role of macrophages in lung adenocarcinoma, published single-cell sequencing data of lung adenocarcinoma patients was analyzed (PMC7210975) (32). Gene distribution and mitochondrial gene expression were screened (Supplementary Figures 3A-C). Cells with mitochondrial expression ≥20% were identified as dead and removed.

Identification of TAMs in total cells
Markers were used to detect the presence of TAMs in the lung adenocarcinoma single-cell datasets (Supplementary Figures 3D-I). TAMs were then extracted and subtype analysis performed to obtain a TAMs subtype map (Supplementary Figure 3I).

Screening of differential expression genes among TAMs subsets
To identify marker genes amongst the different TAM subgroups, samples were screened in "Seurat". Dot and violin plots revealed the top5 marker genes for each TAM subtype (Supplementary Figure 4; Supplementary Table 5_TAM_marker_genes.txt).
Survival curves of high and low macrophage infiltration in lung adenocarcinoma datasets. Horizontal axis: survival time. Vertical axis: survival probability. Color: level of macrophage infiltration. Survival analysis using

Simulation of dynamic changes in macrophages
"Monocle" was used to identify dynamic changes of macrophages in the tumors and cell polarization (Supplementary Figures 5A-C). Cluster 0 could be divided into Cluster 1 and Cluster 2 amongst TAM subtypes. The identified genes were found to regulate differentiation (Supplementary Figure 5D). Gene enrichment analysis on the subtypes of TAM showed that Cluster 2 positively correlated with E2F TARGETS and G2M CHECKPOINT, whilst Cluster 4 negatively correlated with these pathways (Supplementary Figure 5E).

Screening of prognostic factors based on cell-cluster-markers and TAMS-related-genes using univariate cox regression analysis
Markers of each TAM subtype and TAMs-related-genes were used to identify genes related to the prognosis. Samples were divided into high-and low-expression groups according to the median of gene expression. Univariate Cox analysis was performed and survival curves of the top8 prognostic genes were displayed (Figures 2A-H 3.9 Construction of risk score models and evaluation of the prognostic efficacy Based on the "GLMNET" of the R package, LASSO (Least Absolute Shrinkage and Selection Operator) regression analysis was used to construct a regression model for the expression matrix of prognosis related genes corresponding to "Cell-Cluster-Markers" and "Tams-Related-Genes". By analysis, when the value of the freedom degree was 9, the model was accurate (Figures 2I-O; Supplementary Table 7_ forest.univariate_cox.txt, Supplementary Table 7_Signature_ Coef.txt). The calculation formula of the risk score model are listed as follows: Risk Score = 0.0354754835*C1QTNF6 (Expression Value) + 0.0023344103* FSCN1 (Expression Value) + 0.0022298189*GLI2 (Expression Value) + 0.0001616254 * KPNA2 (Expression value) + 0.0005176419*TYMS (Expression Value) + 0.0037498174 *BIRC3 (Expression Value) + 0.0033257017*RBBP7 (Expression Value) + 0.0002465129 *KRT8 (Expression Value) -0.0263442444 *GPR65 (Expression Value). Kaplan-Meier survival curves indicated a significant difference in survival between high and low risk groups. The ROC curve indicated high performance of the risk score model.

Validation of risk score prognostic models in external datasets
To further verify the stability of the risk score model, external and independent data GSE19188, GSE26939, GSE31210, GSE50081, GSE42127 and GSE72094 were used to verify predictive efficacy. Through Kaplan-Meier survival analysis, the constructed risk score model performed well for all external data predictions (Figure 3).

Robust principal component analysis of risk scoring models in clinical factors
To confirm the stability of the risk score model according to clinical characteristics, differences in survival status between high-and low-risk groups in terms of age, gender, radiotherapy, clinical characteristics and Pathlogic M were explored. Significant differences in survival between high-and low-risk groups were observed in those aged ≥ 60 and ≤ 60 years (Figures 4A, B; Supplementary Table 8_clinical_inf.txt). Similar differences were observed between gender subgroups ( Figures 4C, D). In the radiotherapy group, differences between high-and low-risk groups were more pronounced (Figures 4E, F). In Pathologic M ( Figures 4G, I), significant differences between high-and low-risk groups were observed for M0, indicative of higher stability.

Differences in risk score models among cancer clinical factors
To investigate the relationship between the risk score model and clinical characteristics, specific features were selected for analysis. The risk score was found to be related to radiation therapy, pathologic T and tumor stage. No significant relationship to age or gender were observed (Supplementary Figure 6).

Evaluation of risk score models through univariate and multivariate cox regression analysis
To determine whether the risk score model could act as an independent prognostic factor for cancer, the "coxph()" function in the R package "survival" was adopted for univariate and multivariate regression analysis on training and test sets, respectively. We found that in all validation and test sets, the p value of the risk score was ≤ 0.05 ( Figure 5

Construction of a nomogram model of risk scores and clinical factors to predict cancer progression
We next sought to apply the risk scoring model to the prediction of cancer progression in the clinic. The R package "rms" was adopted to construct a nomogram using a variety of clinical features. Calibration curves were used to calculate 1, 2, 3,

Prediction of immunotherapy efficacy amongst subtypes
We next investigated whether the risk score model could predict the prognosis of immunotherapy. Data were calculated using the risk score model and the K-M survival status between high-and low-risk groups evaluated (Supplementary Figure 7A,  Supplementary Table 11_Immune_treatment.xlsx). Upon statistical analysis of the distribution of CR/PR and PD/SD, the proportion of treatment response rates significantly differed between high-and low-risk groups (Supplementary Figure 7B, chi-square test p= 0.004133). No significant differences in the risk scores between the different treatment response groups were observed (Supplementary Figure 7C).

Overall survival analyses of M1 and M2 macrophage subtypes in patients with lung cancer
A total of 32 patients with lung cancer were evaluated for M1 and M2 macrophage subtypes. Samples were stained using double-labeled IHC. The majority of patients were in p at ho lo g ic al St a g e I I ( 6 2 . 5 % ) a n d t h e d o m in an t histopathological type was adenocarcinoma (68.8%). The clinicopathological characteristics of the lung cancer patients are shown in (Table 1).
To identify M1 and M2 macrophage subtypes, CD68 and CD163 antibodies were used for double-labeled IHC staining. CD68 (brown/yellow) as a surface marker for all macrophages primarily localized to the cytoplasm, whilst CD163 (red) localized to the plasma membrane. M2 macrophages were identified through double staining for CD68 and CD163. M1 macrophages were identified through staining with CD68 alone. Representative IHC images are shown in (Figures 6A, B).
The prognostic value of macrophage infiltration was next evaluated. Total macrophages, M2 to M1, and M2 macrophage infiltration were identified as detrimental to patient survival ( Figures 6C, E, F), whilst M1 macrophage infiltration was beneficial to prognosis ( Figure 6D). The infiltration of M1 macrophages in adenocarcinoma was significantly higher than that in squamous cell carcinoma of lung cancer. No significant differences in M2 nor total macrophage infiltration were observed between these two histological subtypes ( Figure 6G).

Discussion
Macrophages with different phenotypes are frequently cited as indicators of the prognosis of lung cancer patients and the efficacy of immunotherapy (17)(18)(19)(20)(21)(22)(23)(24). In our preliminary analyses, macrophage infiltration, rarely reported in lung cancer, had a significant detrimental effect on the prognosis of lung cancer patients ( Figure 1). These data were consistent across cohorts ( Figures 1A-F ) and further verifi ed in follow-up immunohistochemical analysis of clinical samples ( Figure 6C). Collectively, these data highlight how macrophages not only act as innate immune cells to regulate immunological responses (23), but play an important role in the prognosis of lung cancer. This lays the foundation for subsequent module analysis based on macrophage infiltration (Supplementary Figure 1).
Based on the expression profiles of TAMs-related-genes, a consistent clustering profile was constructed ( Supplementary  Figures 2A-E). Significant differences in both survival analysis and PCA (Supplementary Figures 2F-O) were observed. These apparent differences were further identified in single-cell data (Supplementary Figure 3) confirming the importance of macrophages to the prognosis of lung cancer patients (17)(18)(19)(20)(21)(22)(23)(24)32). These data also highlight the need for further refinement of relevant factors to more favorably evaluate patient prognosis.
We used WGCNA to identify macrophage infiltrationrelated module genes and single-cell sequencing of lung adenocarcinoma tissue to identify marker genes of macrophage subtypes. This permitted the construction of a risk assessment model with high prognostic efficacy. The model performed well on external and independent datasets. Immunohistochemistry analysis of clinical samples were consistent with our data. We therefore infer that the risk score has both high clinical practicability and application.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.