We downloaded two keloid RNA chip data sets (GSE44270 and GSE145725) and one keloid single-cell RNA sequencing data set (GSE163973) from the GEO database. After standardization, samples without grouping information were excluded, and finally, 32 samples were obtained in GSE44270, 19 samples in GSE145725, and 6 paired tissue samples in GSE163973.

1. We selected cells with more than 300 expressed genes, less than 5500 expressed genes in total, less than 10 percent mitochondrial genes, and less than 3 percent red blood cell genes. A total of 20,629 cells were retained. Then 3000 hypervariable genes were selected for analysis, and then the number of principal components (PCs) was set to 15 to obtain cell cluster clusters, and then these clusters were displayed in the form of a “tSNE” diagram. We then annotated the resulting cluster. First, we found the top 10 significant markers (Supplementary File 1) in each cluster using the “FindAllmarkers” function. Then all the cells were annotated by the “SingleR” package. SingleR is an R package commonly used for automatic cell type annotation of single-cell RNA-seq (scRNA-seq) data. Spearman correlation was calculated between the expression profile of each cell and that of the reference sample. Second, we define the score for each label as a fixed quantile of the correlation distribution (0.8 by default). Finally, we repeat this for all the labels, and the label with the highest score serves as an annotation for the cell. In this study, we annotated cells using the SingleR package, and finally obtained the 12 clusters (Supplementary File 2).

2. After all the cells were annotated by the “SingleR” package, by using the “subset” function of Monocle R package, we first extracted all the objects of fibroblasts, and randomly selected 1500 cells through the “sample” function for subsequent pseudotime analysis. Then we used the method of “DDRTree” to reduce the dimension of cells, and then calculated the type of cell differentiation state through the function of “reduceDimension”. We used “FindMarkers” function to find genes that had significant differences among different states. The top 5 most significant markers of each State are shown in Supplementary File 3. Finally, we use the “plot_cell_trajectory” function to display the graph of cell differentiation trajectory.

We used the “GOplot” package to perform gene function and pathway enrichment analysis on the differential genes obtained from the single-cell analysis and preserved the gene set with | logFC | >1& P <0.05. First, the top 20 most relevant gene sets were presented in bar chart. Then we extracted the top 5 most related gene sets and presented them in the form of a circle graph.

We conducted WGCNA analysis on GSE44270, selected the top 50 genes with variance, and eliminated the outlier samples. The samples are grouped into different modules by selecting the best soft threshold. At the same time, we also identified conservative modules and eliminated modules with Zsummary.qual<2. Finally, the sample groups were connected, and the correlation analysis between different modules and sample types was carried out by the “Spearson” method.

We conducted a differential analysis of GSE145275. The limma package was used to analyze the differences between groups of keloids and normal tissues after standardization and presented as heat maps. Then we set the filter (| logFC | > 1 & p < 0.05) to get differentially expressed genes, and display them in volcanic graph. Finally, we took the intersection of the differential genes obtained by single-cell analysis, the module genes obtained by WGCNA, and the differential genes obtained by GSE145725.

RNA-sequencing was performed on 3 keloid specimens as well as the patients’ normal skin (at least 1 cm from the keloid lesion). These 3 patients were enrolled between March 2021 and April 2021. All of them signed informed consent forms. And this study was approved by the Ethics Committee of the First Affiliated Hospital of Nanjing Medical University (No. 2021-SR-418). Keloid and the surrounding 1cm of normal tissue were resected by the same surgeon. TruSeq chain mRNA Library Prep Kit (Illumina) was used to generate the Library. Next-generation sequencing was performed on Illumina NovaSeq6000 (Illumina Inc., 100 cycles, single-read sequencing). Sample quality was assessed by FastQC. Illumina NovaSeq6000 analyzed RNA sequencing data for more global analysis of genomic abnormalities. Data is preprocessed using standard pipelines that contain quality control indicators, such as FastQC and MultiQC. Sequence alignment based on STAR-RNA sequencer and sequencing reads allocated to genome features via Featurecots and Vomo-Transformed.

First, expression analysis of potential biomarkers was performed in two datasets, GSE163973 and GSE44270. Subsequently, RNA-sequencing data from three patients at our center were used for validation. Student’s t test was used for statistical analysis.

We assessed the level of immune cell infiltration in all samples of GSE44270 and GSE145275 by using the analysis method of CIBERSORT and set the parameter “PERM” as 1000. Then we divided keloid patients into high-expression group and low-expression group with the median expression value of the potential biomarker as the node and used a violin diagram to show the distribution of immune cells in the high- and low-expression group.

First, we downloaded the gene set of “msigdb.v7.0.entrez.gmt” from GSEA’s official website. Then, after the differential genes were sorted down according to logFC, GSEA analysis was performed, and the result of P <0.05 was retained. Finally, the enrichment results were demonstrated using the function “Gseaplot2”.

We next conducted qRT-PCR experiment on 16 keloid patients, from whom the keloid tissue and peri-keloid normal tissue were taken for mRNA quantification. These 16 patients were enrolled between June 2021 and October 2021. All of them signed informed consent forms. And this study was approved by the Ethics Committee of the First Affiliated Hospital of Nanjing Medical University (No. 2021-SR-418). Keloid and the surrounding 1cm of normal tissue were resected by the same surgeon. Total cellular RNAs were isolated from cells using Trizol Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. The reverse transcription was conducted using the reverse transcription kit provided by Takara (Otsu, Shiga, Japan). Real-time polymerase chain reaction (RT-PCR) was performed using a QuantiTect SYBR Green PCR Kit (Takara), and on a Applied Biosystems QuantStudio 1 (Thermo, Waltham, MA, USA). Relative quantification was determined using the -2ΔΔCt method. The relative expression of messenger RNA (mRNA) for each gene was normalized to the level of glyceraldehyde-3-phosphate dehydrogenase (GAPDH) mRNA. The specific primer sequences adopted in this experiment are: TNC: Forward: 5’- AGTAACGGTGGTGGATTCTGG-3’; Reverse:5ʹ-CTTCCGGTTCGGCTTCTGTA-3ʹ.GAPDH:Forward:5ʹ-GAACGGGAAGCTCACTGG-3ʹ;Reverse: 5ʹ- GCCTGCTTCACCACCTTCT - 3ʹ.

All the statistical analyses were carried out based on R software. Comparisons between keloid and normal groups were performed using the student’s t-test. The Chi-squared test and Fisher’s exact test were performed to compare the clinical variables. ROC analysis was generated to determine the discrimination value of TNC. GSEA was conducted with the clusterProfiler package. Unless otherwise specified, P<0.05 is considered significant.

We performed quality control on the single-cell data set. As shown in Figure 2A, we excluded some cells and limited the percentage of mitochondrial genes and red blood cell genes to ensure the reliability of cell samples. Then we picked out 3,000 highly variable genes and tagged the top 10. In Figure 2B, all the hypervariable genes are marked in red. In Figure 2C, we found that the cells were divided into 12 clusters, and different clusters were roughly summarized as chondrcocytes, fibroblasts, and tissue stem cells. In cluster3, 7, 8, and 9, which we annotated as chondrocytes, we found high expression levels of genes related to skeletal development, ossification, and osteoblasts, such as “COL11A1’,”COMP”, and “POSTN”. By reading the original text of this dataset, they probably should have been more accurately annotated as mesenchymal fibroblasts with high expression of osteogenic genes (8). Finally, we marked the keloid and normal samples as red and blue respectively.

We conducted simulation analysis on the cell trajectory differentiation of all the fibroblasts, and finally found that as shown in Figure 3A, the darker the blue is, the earlier the cell differentiation is, indicating that with the differentiation of time, the fibroblasts differentiate from the left to the right, and the lightest blue is the latest differentiated cell. As shown in Figure 3B, in the process of differentiation of fibroblasts, there were altogether 5 differentiated states, which were marked with different colors, and the red one on the left was found to be the earliest differentiated type. We then investigated the differentiation process between keloid and normal fibroblasts and found that as shown in Figure 3C, keloid fibroblasts differentiated earlier, while normal fibroblasts differentiated later. Figure 3D indicates that all the cells analyzed are fibroblasts.

We found differential genes between keloid fibroblasts and normal fibroblasts were enriched in many biological processes. Figures 4A, C show the results of GO enrichment analysis, which are divided into the biological process (BP), cell component (CC), and molecular function (MF). First of all, in BP, differential genes were mainly enriched in the extracellular matrix organization, extracellular structure organization, ossification and collagen fibril organization. In CC, differential genes were mainly enriched in collagen−containing, endoplasmic reticulum lumen, and other cell components. In MF, differential genes were mainly concentrated in extracellular matrix structural constituent, glycosaminoglycan binding, and extracellular matrix structural constituent conferring tensile strength and other functions. Figures 4B, D show the results of the KEGG enrichment analysis. Differential genes were mainly enriched in protein digestion and absorption, PI3K−Akt signaling pathway, focal adhesion, and ECM−receptor interaction.

By clustering 32 samples in the dataset, we found that the 16 samples on the left and the 16 samples on the right were quite different and clustered into 2 clusters (Figure 5A). After setting the cutoff value to 55, we included the left 16 samples in the subsequent analysis. We found that as the threshold increased, the R^2 value increased and crossed 0.8 (Figure 5B). We chose the optimal threshold of 18. Then, we clustered similar genes into different modules. The results showed that the turquoise module accounted for a larger proportion of the genes (Figure 5C). To verify the accuracy of the WGCNA process and results, the dataset was randomly divided into a training set and a testing set, and conservative sequences were identified and retained. Finally, the values of yellow, red, and tan modules were found to be lower than 2, which were excluded, and the final findings were shown in Figure 5D. Then, through correlation analysis between different modules and phenotype files, it was finally found that purple modules were significantly positively correlated with normal skin fibroblast (correlation = 0.86, P <0.001, Figure 5E), but significantly negatively correlated with keloid fibroblast (correlation = -0.68, P <0.01, Figure 5E). At the same time, there is no significant correlation in non−lesional fibroblast, suggesting that the purple module plays an important role in keloid fibers and normal fibers, which may be related to the occurrence and development of keloid. The genes in the purple module were labeled as the WGCNA-hub genes.

We performed differential expression analysis on another chip data set, GSE145275. As shown in Figure 6A, we found that the expression of these genes differed greatly among different samples, and the redder the color was, the higher the expression level was. In Figure 6B, differentially expressed genes are shown in the form of a volcano map, where blue represents the genes that are low in keloid tissue and red represents the genes that are high in keloid tissue.

The differential genes obtained by single-cell analysis, the differential genes obtained by WGCNA analysis, and the differential genes in the dataset 145725 were collected by Venn graph (Figure 6C). Finally, a key gene TNC was obtained from the intersection, suggesting that this gene may play an important role in the development of keloid. Then, we plotted ROC curves in two chip data sets and found that the AUC of the ROC curve in data set GSE44270 was 0.9 (Figure 6E), and that in data set GSE145725, the AUC was 1.0 (Figure 6D). It is suggested that the TNC gene has a good effect on the diagnosis of keloid fibrous tissue and normal fibrous tissue. To further verify the above results, ROC curves were also drawn in the three pairs of keloid and normal samples sequenced by us, and the AUC value of the area under the curve was 1 (Figure 6F). Although the small sample size may affect the ROC value, it can be seen from the results that this gene also has a good effect on the diagnosis of keloid and normal tissue.

In the single-cell sequencing dataset, we found that the expression of TNC was different between keloid samples and normal samples (P < 2.2E-16, Figure 7A). Subsequently, we found that the expression of this gene was also different in keloid tissue and normal tissue in dataset GSE44270 (Figure 7B). Group 1 represented the keloid group and group 0 was the normal control. TNC expression is up-regulated in keloids. Then we used our RNA sequencing data to verify the expression of TNC in keloid and normal tissue and found that TNC was also up-regulated in keloids (Figure 7C). Raw RNA-sequencing data are summarized in Supplementary File 4. The pseudo-time series analysis mentioned above has found that fibroblasts in normal and keloid cells differ in cell differentiation trajectory, with a total of 5 differentiation states (Figures 3B, C). Next, we studied the significance of TNC in pseudo-time series analysis, as shown in Figure 7D. In the process of fibroblast differentiation, the expression of gene TNC first decreased, then tended to be stable, and finally increased again, suggesting that gene TNC may be related to fibroblast differentiation.

We selected 18 samples of all keloid tissues from two RNA microarray datasets (GSE44270 and GSE145275) and performed immune infiltration analysis. Figure 8B shows the composition of immune cells in each sample. It was found that the dominant immune cells in keloid tissue were T cells. We then performed correlation analysis for these immune cells (Figure 8A), with the darker red indicating a stronger negative correlation and the darker blue indicating the strongest positive correlation. We found that there was a strong correlation between NK cell activated and NK cell resting, and a strong correlation between Mast cell activated and NK cell activated in keloid tissue. Mast cells activated cells and NK cells resting cells had a strong negative correlation, while Macrophage M2 and T cells follicular helper had the strongest negative correlation. Finally, according to the expression of TNC, we divided them into TNC high expression group and TNC low expression group. The infiltration level of macrophage M0 was different between the two groups (Figure 8C). Macrophages M0 had a higher infiltration level in low TNC- group (P<0.05). By using the “ImmuneSubtypeClassifier” R package, keloid patients were divided into different immune subtypes. By evaluating each sample, the best corresponding subtype was determined, and two subtypes were found, namely immune subtype C1 and immune subtype C2 (Figure 8D). Then, we explored the expression of TNC among immune subtypes and found that there was a trend of high expression of TNC in subtype C1, but P =0.18 indicated that the difference was not statistically significant, which may be caused by the small sample size (Figure 8D).

GSEA analysis found that, in the group with high expression of TNC in keloids, the enrichment pathways included skin development, fibroblast, cancer progression, neutrophil at skin wound, and uterine fibroid (Figure 9A); In the low TNC expression group, the enriched pathways included YBX1 targets, ES-1, proliferation, transformed by RhoA, and defetive CFTR cause cystic fibrosis (Figure 9B).

We next conducted qRT-PCR experiment on 16 keloid patients, from whom the keloid tissue and peri-keloid normal tissue were taken for TNC mRNA quantification. Figure 10 showed that the TNC expression in keloid was significantly higher than the Peri-keloid normal tissue from the same patient. (***p<0.001).