On the GEO database (http://www.ncbi.nlm.nih.gov/geo/), the datasets meeting the following criteria were included: (1) datasets with data of lung function and (2) datasets with prognostic data (also including the data of decline of lung function). In addition, in order to compare the CD247 expression between IPF and control and clarify the change process of Cd247 expression of T cells in the process of fibrosis, the GSE110147 (15), GSE33566 (16), and GSE141259 (17) datasets (bleomycin-induced mouse model of pulmonary fibrosis) were chosen. In the GSE141259 dataset, single-cell RNA-sequencing (scRNA-seq) data of 0, 3, 7, 10, 14, 21, and 28 days were provided after bleomycin was used. Finally, we selected 11 datasets: four came from lung tissue samples (GSE32537 (18), GSE47460 (19), GSE110147 (15), GSE141259 (17)), six came from blood samples (GSE93606 (20), GSE27957 (21), GSE28042 (21), GSE38958 (22), GSE132607 (23), GSE33566 (16)), and one came from bronchoalveolar lavage fluid (BALF, GSE70866 (24)). Approval of the Ethics Committee was not necessary for the datasets from the GEO database.

The gene expression matrix of GSE32537 (18), GSE47460 (19), GSE110147 (15), GSE93606 (20), GSE28042 (21), GSE38958 (22), GSE132607 (23), GSE33566 (16), and GSE70866 (24) and the raw data of GSE27957 (21) were downloaded from GEO. “Affy” package (25) (v.1.68.0) was used to normalize the array data according to the robust multi-array average (RMA) method (GSE27957). In addition, lung function [percent predicted forced vital capacity (FVC% predicted) and percent predicted diffusion capacity of the lung for carbon monoxide (Dlco% predicted)] was extracted from the GSE93606, GSE38958, GSE132607, GSE32537, and GSE47460 datasets ( Tables 1 , 2 ). The data of lung function were not complete in the GSE33566 dataset. Therefore, the data of lung function were not used in this study. Transplant-free survival (TFS) was extracted from the GSE27957 and GSE28042 datasets ( Table 1 ). The composite end point (CEP, death, or decline in FVC >10% over a 6-month period) and survival were extracted from the GSE93606 and GSE70866 datasets, respectively ( Tables 1 , 2 and Supplementary Table 1 ). Furthermore, follow-up data were also extracted from the GSE132607 (transcriptome data and lung function data) and GSE93606 datasets (transcriptome data), respectively.

Pearson correlation coefficients were used to determine the association between CD247 and other genes or lung function according to R package “stats” (v.4.0.5) among these datasets. R packages “gplots” (v.3.1.1), “pheatmap” (v.1.0.12), and “RColorBrewer (v.1.1-2)” were used to construct the heatmap. R package “forestplot” (v.1.10.1) was used to construct the forest plot (https://CRAN.R-project.org/package=forestplot).

The computational analysis of the GSE141259 dataset was performed using R package “Seurat” (4.0.3) (26). Quality control had been finished by the authors of this GSE141259 dataset; therefore, 29,297 cells were analyzed. The Seurat SCTransform () function was used to normalize the scRNA-seq data. Principal component analysis (PCA) was calculated using the Seurat RunPCA () function. UMAP embedding and Louvain clusters were calculated using the first 50 principal components with the Seurat RunUMAP () and FindClusters () functions, respectively. Resolution was set as 0.9. The Seurat FindAllMarkers () function was used to find markers of 31 clusters, and cell types were identified based on markers of each cluster according to the CellMarker (27) and PanglaoDB databases (28). Expression and distribution of Cd247 were visualized according to Seurat DotPlot () and FeaturePlot () functions. Dynamic Cd247 expression was visualized by the web (https://theislab.github.io/LungInjuryRegeneration/).

We used STRING (http://www.string.embl.de/, version: 11.0 b) (29) to analyze proteins that purportedly interact with CD247 [“medium confidence (0.400)”, meaning of network edges (“evidence”), max number of interactors to show (“no more than 50 interactors” in 1st shell), and active interaction sources (“experiments”)]. The MEM database (30) (https://biit.cs.ut.ee/mem/index.cgi) was used to verify the correlation between CD247 and genes obtained from STRING.

Differentially expressed genes (DEGs) were defined as expression levels of genes that were significantly diverse in IPF patients with low-expression CD247 compared with those with high-expression CD247 (|log Fold Change|>0.5 and false discovery rates (FDR) < 0.05). “Limma” package (v.3.46.0) (31) was used for the analysis of DEGs. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were analyzed and visualized according to R package “clusterProfiler” (32). P values were adjusted with the Benjamini and Hochberg (BH) correction method. The single-sample gene set enrichment analysis (ssGSEA) score of 19 immune cells and the activity of 15 immune-related pathways (33) were calculated by the “GSVA” R package (v.1.38.2) (34). CIBERSORT is a large-scale analysis tool of RNA mixtures for cellular biomarkers and therapeutic targets according to the gene expression feature sets of 22 immune cell subtypes (http://cibersort.stanford.edu/) (35). Subsequently, the 20 immune cell subtypes were classified into four types: lymphocytes (B cells naive, B cells memory, plasma cells, T cells CD8, T cells CD4 naive, T cells CD4 memory resting, T cells CD4 memory activated, T cells follicular helper, T cells regulatory, T cells gamma delta, NK cells resting, NK cells activated), macrophages (monocytes, macrophages M0, macrophages M1, macrophages M2), dendritic cells (dendritic cells resting, dendritic cells activated), and mast cells (mast cells resting and mast cells activated).

SPSS Statistics 23 (IBM SPSS) and R software (Version 4.0.3) were used for statistical analysis. Categorical variables were described as number (%) and were compared by the chi-square test or the Fisher exact test. Continuous variables were compared by independent group t tests. The diagnostic accuracy of CD247 for IPF or Dlco% predicted decline ≥15% over 12 months (Dlco15) were estimated by receiver operating characteristic (ROC) analysis. Youden index was used to calculate the optimal cutoff values of CD247 for Dlco15. Logistic regression analysis was used to estimate the odds ratio (OR) for the statistically significant correlation between gene expression and Dlco15. Kaplan–Meier analysis was used to compare the TFS or CEP or survival between different groups. R package “survminer” (https://CRAN.R-project.org/package=survminer, v.0.4.8) was used to estimate the optimal cutoff expression value of CD247 for the survival analysis. The hazard ratio (HR) for the statistically significant correlation between CD247 expression and TFS or CEP was estimated by univariate Cox regression. The R package “survivalROC” (https://CRAN.R-project.org/package=survivalROC, v.1.0.3) was used to construct the time-dependent ROC curve to evaluate the predictive value of CD247. Some statistical analyses were visualized by GraphPad Prism 9. The bilateral test was used.

The association between CD247 expression and lung function was respectively determined in the GSE38958, GSE132607, and GSE93606 datasets (blood samples). The significant correlations (p < 0.05) between genes and lung function are shown in Supplemental Material 1 . After intersecting the genes with significant correlations in the three blood datasets, 33 genes associated with Dlco% predicted and 17 genes associated with FVC% predicted were selected for further analysis ( Supplementary Tables 2 , 3 and Figure 1 ). Furthermore, of the 33 genes, 14 genes with consistent positive correlations (ATP10A, CD244, CD247, DDX19A, DIS3L, OSBPL3, RFX5, SH3YL1, TDRKH, TMTC4, TTC39B, TXK, UBE3C, and UTP15) and 11 genes with consistent negative correlations (BNIPL, C9orf131, DUSP13, GJA5, GUCY2D, MYL4, RHAG, SLC28A1, SLC6A7, TNFRSF19, and TULP2) were selected as the candidate hub genes ( Figure 1B , pink rectangles). Of the 17 genes, one gene with consistent positive correlations (PMPCB) and five genes with consistent negative correlations (ALLC, ANLN, DRD3, GLT8D2, and NECAB1) were also selected as the candidate hub genes ( Figure 1D , pink rectangles).

The follow-up data (4, 8, 12 months) were extracted from the GSE132607 dataset (blood). In the fourth month, only MYL4 was significantly negatively associated with Dlco% predicted ( Table 3 ). In the eighth month, six genes (CD247, DIS3L, OSBPL3, TDRKH, TTC39B, and UTP15) were significantly positively associated with Dlco% predicted, and three genes (GUCY2D, MYL4, and RHAG) were significantly negatively associated with Dlco% predicted ( Table 3 ). In the 12th month, five genes (CD244, CD247, DDX19A, RFX5, and TXK) were significantly positively associated with Dlco% predicted, and MYL4 was significantly negatively associated with Dlco% predicted ( Table 3 ). However, no genes were significantly associated with FVC% predicted in the follow-up data. Therefore, CD247 and MYL4 with a greatly consistent significant correlation were chosen for further analysis.

Dlco% predicted decline ≥15% over 12 months (Dlco15) is a useful index to evaluate whether patients with IPF were progressive. Patients with Dlco15 were more likely to have the lower CD247 expression and higher MYL4 expression compared with those without Dlco15 in the GSE132607 dataset ( Figures 2A, C and Supplementary Table 4 ). In addition, Dlco% predicted in patients with low-expression CD247 classified by the median value of CD247 at baseline was significantly lower than that in patients with high-expression CD247 after visiting at 4, 8, and 12 months in the GSE132607 dataset ( Table 4 and Figure 2B ). Likewise, Dlco% predicted in patients with high-expression MYL4 classified by the optimal cutoff value (4.37) of MYL4 at baseline was significantly lower than that in patients with low-expression MYL4 after visiting at 4, 8, and 12 months ( Table 4 and Figure 2D ). Subsequently, two lung tissue datasets (GSE32537 and GSE47460) were used. As shown in Figure 3 , CD247 expression was significantly positively associated with Dlco% predicted in both blood and lung tissue. However, MYL4 was significantly negatively associated with Dlco% predicted in blood samples, not in lung tissue samples ( Supplementary Figure 1 ). Furthermore, CD247 was not significantly correlated with GAP (gender, age, and physiological index) in the GSE70866 dataset (BALF sample, data not shown).

The expression of CD247 in patients with IPF was lower than that in controls in the GSE93606 (blood), GSE33566 (blood), GSE47460 (lung tissue), and GSE110147 (lung tissue) datasets, whereas the expression of CD247 did not show a difference in the GSE38958 (blood) and GSE32537 (lung tissue) datasets ( Figure 4 ). The expression of CD247 in patients with IPF (n = 75) was lower than that in controls (n = 19) in the GSE28042 dataset (13.24 vs. 13.56, p = 0.0043). Immunosuppressive therapy was not used before the blood samples were collected in the GSE93606 and GSE33566 datasets, implying that the expression of CD247 in the blood samples was not affected by the immunosuppressive therapy ( Table 1 ). The diagnostic values of CD247 for IPF were variable in different datasets ( Supplementary Figure 2 ). Besides, the expression of MYL4 did not show the significant difference in these datasets (data not shown). Therefore, CD247 was considered as the key gene.

According to the ROC curve analysis and R package “survminer”, the optimal cutoff value of CD247 was chosen for the logistic regression analysis in the GSE132607 (blood) dataset and prognosis-related analysis in the GSE93606 (blood), GSE27957 (blood), and GSE28042 (blood) datasets, respectively. According to logistic regression analysis, a low expression of CD247 at visiting at 0, 8, and 12 months was the risk factor of Dlco15 in the GSE132607 dataset, whereas the low expression of CD247 was not the risk factor of Dlco15 in the fourth month ( Figure 5A ). According to Cox regression analysis and Kaplan–Meier analysis, low-expression CD247 at visiting at 0, 1, 3, 6, and 12 months was the risk factor of CEP in the GSE93606 dataset and significantly associated with shorter TFS time in the GSE27957 and GSE28042 datasets ( Figure 5B , Supplementary Figures 3A–G ). Furthermore, the ROC curve showed that the areas under the curve (AUC) were 0.736 at 1 year and 0.741 at 2 years for CEP in the GSE93606 dataset ( Supplementary Figure 3H ). Moreover, AUCs for non-TFS were respectively 0.889, 0.787, and 0.702 at 1, 2, and 3 years in the GSE27957 dataset, whereas the AUC was relatively low in the GSE28042 dataset ( Supplementary Figures 3I, J ). In the GSE70866 dataset (BALF samples), CD247 was not significantly associated with mortality (p > 0.05, Supplementary Table 5 ).

The protein–protein interaction (PPI) network was constructed based on CD247 according to the STRING database (average node degree: 3.71, PPI enrichment p-value: <1.62e^-3; Figure 6A ) (29). According to the PPI network, 13 protein-coding genes (CD3E, SHC2, DOCK2, JAK3, PTPN6, LCK, FYN, SHC1, CSK, PTPN3, ZAP70, SHC4, and SHC3) were found to interact with CD247. The MEM database (30) gathers several hundreds of publicly available gene expression data sets from the ArrayExpress database, which is a useful tool for the co-expression analysis across hundreds of datasets. Results from the MEM database (30) revealed a significant co-expression relationship between CD247 and the 9 of 13 interacted genes (CD3E, ZAP70, LCK, FYN, JAK3, DOCK2, PTPN6, CSK, and SHC1; Figure 6B ). In addition, in patients with IPF, results from blood, lung tissue, and BALF samples revealed a significant co-expression association between CD247 and the 6 of 13 interacted genes (CD3E, ZAP70, LCK, FYN, JAK3, and PTPN3; Figure 6C ). GO and KEGG analysis implied that CD247 and its 13 interacted genes were enriched in the immune response especially the T cell-related biological processes (BP), cellular component (CC), and pathways ( Figures 6D, E ). Furthermore, CD247, ZAP70, LCK, FYN, PTPN6, SHC1, SHC4, SHC2, and SHC3 were associated with NK cell-mediated cytotoxicity ( Supplementary Table 6 ). Interestingly, five genes (CD247, CD3E, ZAP70, LCK, PTPN6) were associated with the hsa05235 pathway (PD-L1 expression and PD-1 checkpoint pathway in cancer, Supplementary Table 6 ).

In order to reveal the biological significance correlated with CD247, the DEGs between the patients with high-expression and low-expression CD247 based on the median value were used to conduct the GO enrichment and KEGG pathway analysis in the GSE38958, GSE132607, and GSE93606 datasets (blood samples). In the three datasets, DEGs were mainly enriched in inflammation- and immune-related response and pathways such as neutrophil activation involved in immune response, T cell activation, T cell differentiation, leukocyte chemotaxis, IL-17 signaling pathway, T cell receptor signaling pathway, Th17 cell differentiation, PD-L1 expression, and PD-1 checkpoint pathway in cancer and so on ( Figure 7 ). Furthermore, related functions or pathways with ssGSEA showed that patients with low-expression CD247 were more likely to have lower immune activity [lower T cell general, Th1 cells, TIL (tumor-infiltrating lymphocytes), checkpoint, cytolytic activity, T cell co-inhibition, and T cell exhaustion scores)] and higher degree of inflammation response [higher DCs (dendritic cell), M2 macrophages, and neutrophils] compared with patients with high-expression CD247 in the three datasets ( Figure 8 ), which was consistent with the results of CIBERSORT analysis ( Supplementary Figures 4A, C, E ). In addition, low-expression CD247 was significantly associated with low lymphocyte score ( Supplementary Figures 4B, D, F ).

According to the LungMAP database (36), CD247 is mainly expressed by the T cells and NK cells in the human lung based on scRNA-seq data (37) ( Supplementary Figure 5A ). In the mouse lung, Cd247 is also mainly expressed by the T cells and NK cells, and the expression of Cd247 of T cells is increased in the acute inflammatory stage, then decreased in the fibrotic stage after bleomycin injection ( Supplementary Figures 5B–D ).