Identification of Diagnostic Biomarkers Associated with Stromal and Immune Cell Infiltration in Fatty Infiltration After Rotator Cuff Tear by Integrating Bioinformatic Analysis and Machine-Learning

Introduction

Rotator cuff tear (RCT) is a common shoulder disorder causing shoulder pain and disability. The prevalence of full-thickness RCT is 20.7% in the general population, and increased with age.¹ Rotator cuff play essential roles in shoulder function and the treatment of proximal humeral fractures.^2,3 It is important to repair RCT. In the past decades, there was a fast-growing number of rotator cuff repair procedures.^4–6 However, massive RCT still has a high failure rate after rotator cuff repair.^7–9 The retear rate after repair of massive RCT could reach up to 36%.⁸ Fatty infiltration (FI) of the rotator cuff muscles was first proposed by Goutallier et al.¹⁰ It is associated with RCT and rotator cuff muscle atrophy, characterized as fat tissue fills in the interspaces of remodeled muscle fibers.¹¹ FI is an independent risk factor for poor outcomes and tendon retear after rotator cuff repair.^12–15 Animal studies showed that FI originated as early as two weeks after RCT and deteriorated six weeks later,¹⁶ and could persist after delayed rotator cuff repair.¹⁷ Magnetic resonance imaging studies illustrated that FI progressed even after successful rotator cuff repair.^15,18 In view of delayed diagnosis of FI may result in poor prognosis, the identification of diagnostic biomarkers may boost the development of early diagnostic method for FI.

To investigate the molecular mechanisms underlying FI will help identify diagnostic biomarkers and boost the development of therapeutic modalities. In recent years, several studies have showed that stromal and immune cells play important roles in the development of FI. Lakemeier demonstrated that vessel density was significantly correlated with FI after RCT.¹⁹ Liu et al showed that muscle fibrosis and FI coexisted after RCT, and muscle resident progenitors would differentiate into fibroblasts in fatty degenerative rotator cuff muscles.²⁰ Gumucio et al corroborated that macrophages recruitment and accretion of inflammatory cytokines such as IL-10 and IL-1β were observed in fatty infiltrated rotator cuff muscles.²¹ Mendias et al also observed macrophage infiltration in fatty degenerative muscles after RCT.²² Davies et al revealed that upregulated RhoA signaling with increased macrophage infiltration aggravated FI after RCT.²³ Thus, it is meaningful to assess stromal and immune cell infiltration in FI after RCT to clarify the molecular mechanisms and help identify molecular biomarkers and explore therapeutic modalities of FI. Murine microenvironment cell population counter (mMCP-counter) is a technique built on transcriptomic data and applied to quantify immune and stromal murine cell populations.²⁴ To date, few studies have conducted mMCP-counter to analysis stromal and immune cell infiltration in FI and explore candidate diagnostic biomarkers for FI.

High-throughput RNA sequencing (RNA-seq) analysis, together with bioinformatics analysis, has been widely used to screen feature genes related to multiple diseases that might serve as potential therapeutic targets or diagnostic biomarkers.^25–29 In recent years, machine learning algorithms have been proven to be effective ways for feature gene selection.^30–32 Random forest and multiple support vector machine recursive feature elimination (mSVM-RFE) are two popular machine learning algorithms, which can be used to make accurate predictions on high-dimension gene expression data.^33–35

In the present study, the murine RNA-seq datasets of FI after RCT in Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo) database were downloaded. Comprehensive bioinformatics analysis was carried out to identify differentially expressed genes (DEGs), and analyze gene set enrichment analysis (GSEA) of gene ontology (GO) and GSEA of Kyoto encyclopedia of genes and genomes (KEGG). Subsequently, random forest and mSVM-RFE algorithms were first combined to screen diagnostic biomarkers of FI. In addition, mMCP-counter method was first employed to further quantify stromal and immune cell infiltration in FI and normal samples. Moreover, the correlation between biomarkers and infiltrated immune and stromal cell subsets in FI samples were calculated to further analyze the molecular and immunologic mechanisms underlying FI.

Materials and Methods

Data Download

The RNA-seq count datasets of GSE130447³⁶ and GSE103266³⁷ were downloaded from GEO database, which were based on GPL17021 Illumina HiSeq 2500 (Mus musculus) platform and GPL18694 Illumina HiSeq 2500 (Rattus norvegicus) platform, respectively. GSE130447 dataset contains sixteen infraspinatus muscle samples collected from C57BL/6J (B6) and MRL/MpJ (MRL) mice, including eight fatty infiltrated samples established by a full-thickness infraspinatus tendon tear and eight control samples. GSE103266 dataset served as the validation cohort, contains sixteen supraspinatus muscle samples collected from Sprague-Dawley rats, including twelve fatty infiltrated samples made by a full-thickness supraspinatus tendon tear and four control samples.

RNA-Seq Data Processing and DEGs Identification

The regularized log (rlog) transformation depend on the DESeq2 package of R software (version 4.0.3, http://r-project.org/),³⁸ was applied to centralize the RNA-seq count datasets of GSE130447 and GSE103266. DESeq2, ggplot2,³⁹ gplots (https://CRAN.R-project.org/package=gplots) and RColorBrewer (https://CRAN.R-project.org/package=RColorBrewer) packages in R were implemented to analyze and visualize the quality of GSE130447 dataset. DEGs between rotator cuff tear (RCT) and control samples were winkled out using the DESeq2 package, and illustrated with a volcano plot graphed using ggplot2 and ggrepel (https://CRAN.R-project.org/package=ggrepel) packages. Adjusted p-value <0.05 and |log2FoldChange| >1 were set as the cut-off criteria.

GSEA of DEGs

GSEAs of GO and KEGG were conducted to identify the biological functions and pathways of DEGs using ClusterProfiler⁴⁰ and org.Mm.eg.db (https://bioconductor.org/packages/org.Mm.eg.db/) packages, and visualized using ggpubr (https://CRAN.R-project.org/package=ggpubr) and enrichplot (https://github.com/YuLab-SMU/enrichplot) packages. Adjusted p-value <0.05 was set as the cut-off criteria.

Diagnostic Biomarkers Selection and Validation

In order to screen the best underlying biomarkers for fatty infiltration of rotator cuff muscles, we used two varying machine learning algorithms. Random forest and mSVM-RFE algorithms were employed to fulfill feature selection. Feature selection built on random forest algorithm was carried out using the varSelRF package.⁴¹ Feature selection based on mSVM-RFE algorithm was performed by the e1071 package (https://CRAN.R-project.org/package=e1071). The intersecting variables between these two algorithms were obtained as the biomarkers. The diagnostic efficiencies of the screened biomarkers were further validated by the area under the ROC curve (AUC) in the GSE103266 dataset using pROC,⁴² ggplot2 and glmnet⁴³ packages.

Evaluation of Stromal and Immune Cell Infiltration

On account of the heterogeneity of stromal and immune cell infiltration in human and mouse, we adopted the mMCP-counter package,²⁴ a package specific for estimating the composition of immune and stromal cells in murine samples, to assess stromal and immune cell infiltration of the RNA-seq dataset of GSE130447. The ggstatsplot package (https://CRAN.R-project.org/package=ggstatsplot) was implemented to assess and visualize the correlation of immune and stromal cells, and to plot box/violin graphs of the differences of these cells between RCT and control groups.

Correlation Analysis Between Screened Biomarkers and Immune and Stromal Cells

The Spearman’s rank correlation tests between the screened biomarkers and the immune and stromal cells in RCT samples of GSE130447 were explored and visualized with scatterplots using the ggstatsplot package.

Results

RNA-Seq Data Processing and DEGs Identification

After the rlog transformation produce, the RNA-Seq values in GSE130447 ranged from 2.608 to 21.004. A density plot and a box plot displayed that the distributions of samples were homogeneous (Figure 1A and B). A hierarchical clustering heatmap built on the Euclidean distance between samples and a principal-components analysis (PCA) plot evinced that rotator cuff tear and control samples can easily be separated (Figure 1C and D). These results substantiated that the RNA-seq data was credible for further analyses. A total of 2123 DEGs, comprising 1417 upregulated genes and 706 downregulated genes, were identified in GSE130447 (Figure 2).

Figure 1 Distribution analyses of the samples. (A) The density plot of the GSE130447 data distribution. (B) The box plot of the GSE130447 data distribution. (C) The hierarchical clustering heatmap of the GSE130447 data. (D) The PCA plot of the GSE130447 data.

Figure 2 The volcano plot of DEGs in GSE130447.

GSEA of DEGs

GSEA of GO revealed that DEGs had enrichments for biological processes (BP), cellular component (CC) and molecular function (MF) (Figure 3). For BP, DEGs were predominantly enriched in immune system process, extracellular matrix organization, extracellular structure organization, immune response, regulation of immune system process, defense response, positive regulation of response to stimulus, positive regulation of immune system process, response to external stimulus and blood vessel morphogenesis (Figure 4). For CC, DEGs were principally enriched in extracellular region, extracellular space, extracellular matrix, collagen-containing extracellular matrix, basement membrane, cell surface, transferase complex, collagen trimer, cation channel complex and ion channel complex (Figure 5). For MF, DEGs were primarily enriched in extracellular matrix structural constituent, glycosaminoglycan binding, structural molecule activity, endopeptidase activity, heparin binding, antigen binding, ubiquitin-protein transferase activity, ubiquitin-like protein transferase activity, nucleic acid binding and fatty acid binding (Figure 6). In terms of GSEA of KEGG, DEGs were mainly enriched in ECM-receptor interaction, staphylococcus aureus infection, protein digestion and absorption, viral protein interaction with cytokine and cytokine receptor, phagosome, PPAR signaling pathway and osteoclast differentiation (Figure 7).

Figure 3 The histogram of GO analysis of DEGs in GSE130447.

Figure 4 The GSEA plot of BP.

Figure 5 The GSEA plot of CC.

Figure 6 The GSEA plot of MF.

Figure 7 The GSEA plot of KEGG.

Diagnostic Biomarkers Selection and Validation

The top 100 upregulated and top 100 downregulated DEGs were selected to perform feature selection. A subset of 64 features were screened using mSVM-RFE algorithm and 4 variables were identified using random forest algorithm (Figure 8A and B). Fatty acid binding protein 5 (FABP5) and matrix Gla protein (MGP) were ultimately selected from the intersection of features achieved by these two algorithms, and were considered as the diagnostic biomarkers. The diagnostic efficiencies of FABP5 and MGP were further validated using the validation set of GSE103266 RNA-seq dataset. The AUCs were 0.958 (95% CI 0.860–1) in FABP5 and 1 (95% CI 1–1) in MGP (Figure 9). These results illustrated that FABP5 and MGP had excellent diagnostic abilities.

Figure 8 Screening of diagnostic biomarkers using mSVM-RFE and Random forest. (A) mSVM-RFE algorithm for identifying biomarkers. When the number of features was sixty-four, the generalization error reached its minimum value. (B) Random forest algorithm to screen biomarkers. When the number of variables was four, the variables out-of-bag (OOB) error rate reached the minimum value.

Figure 9 The ROC curve of FABP5 and MGP.

Stromal and Immune Cell Infiltration

When conducted stromal and immune cell infiltration analysis in GSE130447 dataset, seven types of immune and stromal cells were detected. The correlation of seven types of immune and stromal cells was estimated (Figure 10). Monocytes were positively correlated with mast cells and eosinophils, whereas negatively correlated with vessels and endothelial cells. Mast cells were positively correlated with monocytes, but negatively correlated with vessels, endothelial cells and fibroblasts. Eosinophils were positively correlated with monocytes. Vessels were positively correlated with endothelial cells and fibroblasts, yet negatively correlated with monocytes and mast cells. Lymphatics were positively correlated with fibroblasts. Endothelial cells were positively correlated with vessels and fibroblasts, whereas negatively correlated with monocytes and mast cells. Fibroblasts were positively correlated with vessels, lymphatics and endothelial cells, but negatively correlated with mast cells. The exploration to the composition of immune and stromal cells in RCT samples versus control samples in GSE130447 dataset demonstrated that the proportions of monocytes (P = 0.005) and mast cells (P = 0.001) in RCT samples were lower than in control samples, and the proportions of vessels (P = 0.003), endothelial cells (P = 0.001) and fibroblasts (P = 0.003) in RCT samples were higher than in control samples (Figure 11).

Figure 10 The Correlation matrix plot of seven types of immune and stromal cells in GSE130447.

Figure 11 The box/violin plot of the differences of immune and stromal cells between RCT and control groups in GSE130447.

Correlation Analysis Between FABP5, MGP and Immune and Stromal Cells

FABP5 was positively correlated with vessels (ρ = 0.81, P = 0.015), whereas negatively correlated with monocytes (ρ = −0.71, P = 0.047) and mast cells (ρ = −0.86, P = 0.007) in RCT samples of GSE130447 (Figure 12). MGP was positively correlated with vessels (ρ = 0.86, P = 0.007), while negatively correlated with monocytes (ρ = −0.76, P = 0.028) and mast cells (ρ = −0.90, P = 0.002) in RCT samples of GSE130447 (Figure 13).

Figure 12 The scatterplots of correlation analysis between FABP5 and immune and stromal cells in GSE130447.

Figure 13 The scatterplots of correlation analysis between MGP and immune and stromal cells in GSE130447.

Discussion

FI of the rotator cuff muscles is a common complication after RCT. To explore molecular biomarkers and stromal and immune cell infiltration in FI may help the development of diagnostic and therapeutic strategies. To the best of our knowledge, the present study was the first to probe into diagnostic biomarkers in combination with stromal and immune cell infiltration in FI after RCT by mining RNA-seq datasets. We download two RNA-seq datasets of FI from GEO database and carried out an integrated bioinformatics analysis of the datasets.

A total of 2123 DEGs were identified, comprising 1417 upregulated genes and 706 downregulated genes. GSEA of GO revealed that DEGs were predominantly enriched in immune system process, extracellular region and extracellular matrix structural constituent, indicating that extracellular matrix (ECM) and immune response participated in the process of FI. In terms of GSEA of KEGG, DEGs were mainly involved in ECM-receptor interaction, staphylococcus aureus infection, protein digestion and absorption, viral protein interaction with cytokine and cytokine receptor, phagosome, PPAR signaling pathway and osteoclast differentiation. ECM is an elementary constituent of skeletal muscle; excessive cumulation of ECM components can result in muscle fibrosis.⁴⁴ FI and fibrosis of rotator cuff muscles could coexist after RCT.²⁰ Peroxisome proliferator-activated receptors (PPARs) with isotypes of PPARα, PPARβ/δ and PPARγ, are members of the ligand-inducible nuclear hormone receptor superfamily, and can be activated by fatty acids.⁴⁵ PPAR signaling pathway regulates a broad spectrum of cellular processes, such as fatty acid transport, oxidation and storage, apogenesis, inflammation, and can regulate immune function by triggering macrophages.^46,47 PPARs also play key roles in skeletal muscle lipid metabolism.⁴⁸ Itoigawa et al demonstrated that PPARγ was high expressed in fatty degenerative rotator cuff muscle after RCT.⁴⁹ Liu et al showed that PPAR signaling pathway was crucial for the regulation of muscle lipid accumulation in broiler chickens.⁵⁰

Based on random forest and mSVM-RFE algorithms, two diagnostic biomarkers were screened. These two biomarkers, FABP5 and MGP, exhibited excellent diagnostic abilities in validation dataset. FABP5 is an important member of fatty acid binding proteins (FABPs) family. FABPs are conservative family of cytoplasmic 14- to 15-kDa proteins that bind medium- and long-chain fatty acids; FABPs have 12 isoforms with highly conservative three-dimensional structures, and play essential roles in fatty acid uptake, transport, utilization, and are closely related with inflammatory responses and metabolic regulation.^51–53 Lee et al showed that RCT would trigger expression of fatty acid binding protein 4 (FABP4), another member of FABPs, causing intramuscular fat deposition in rotator cuff muscles after RCT.⁵⁴ FABP5 was expressed in many tissues, such as skin, lymph node, adipose tissue and skeletal muscle.^52,55,56 It was closely related to immune system. Field et al demonstrated that FABP5 could regulate regulatory T cells function, and its inhibition in regulate regulatory T cells jeopardized lipid metabolism, leading to mitochondrial DNA releasing and activation of the cyclic GMP-AMP synthase-stimulator of interferon genes signaling pathway, inducing high expression of the regulatory cytokine interleukin-10 and suppressing regulatory T cells activity.⁵⁷ Pan et al revealed that FABP5, together with FABP4, was essential for the survival and function of CD8⁺ tissue-resident memory T cells, while the shortage of FABP5 and FABP4 undermined the uptake of exogenous lipid in CD8⁺ tissue-resident memory T cells and lowered their longevity.⁵⁸ Abplanalp et al reported that in cardiac failure patients, FABP5 would be largely concentrated in classical monocyte and contribute to boost monocyte activation.⁵⁹ Moore et al found that FABP5 could regulate macrophage phenotypes, its deficiency would raise the production of anti-inflammatory cytokines.⁶⁰ Bogdan et al corroborated that suppression of FAPB5 would debilitate the induction of microsomal prostaglandin E synthase-1 and the synthesis of prostaglandin E₂ in macrophage during inflammatory response.⁶¹ Moreover, FABP5 was closely linked to skeletal muscle fat metabolism.⁶² Iso et al illustrated that FABP5 was highly expressed in capillary endothelial cells of skeletal muscle, and was crucial for skeletal muscle lipid uptake.⁵⁵ In addition, FABP5 expression was found to be positively correlated with PPARγ expression during adipogenesis.⁶³ MGP is a small 14.7 kDa vitamin K-dependent and γ-carboxyglutamic acid-rich protein, secreted and localized in the ECM of endothelial cells or chondrocytes.⁶⁴ It was widely expressed in many tissues, and played important roles in inhibiting tissue calcification and regulating the process of cell differentiation.^64,65 Velleman et al demonstrated that MGP controlled skeletal muscle cell differentiation and proliferation.⁶⁶ Mutch et al reported that MGP protein was abundantly secreted from adipocytes during adipogenesis, and might be a new adipokine in human.⁶⁷ Lanham et al showed that MGP knockout mice had a higher proportion of intrascapular brown adipose tissue volume than wild-type mice.⁶⁸ Recently, Li et al first reported the link between MGP and lipid metabolism in mouse model, and showed that MGP could regulate adipogenesis and lipolysis in vitro and lipid metabolism in vivo, and found that the level of serum desphosphorylated-uncarboxylated MGP was positively related to visceral adiposity and closely linked to lipoprotein level.⁶⁹ The aforementioned studies suggested that both FABP5 and MGP were closely involved in skeletal lipid metabolism, and might be potential biomarkers for FI.

The features of stromal and immune cell infiltration in FI were assessed using the mMCP-counter method. As a result, various immune and stromal cell subtypes were closely linked to FI. An increased infiltration of vessels, endothelial cells and fibroblasts, and a decreased infiltration of monocytes and mast cells were detected to be potentially linked to the occurrence and development of FI. Gibbons et al observed increased vessels in fatty degenerative rotator cuff muscles.⁷⁰ Endothelial cells were key components of the vessel wall. Muscle fibrosis and FI were common complications after RCT. Liu et al showed that Tie2⁺ progenitor cells and platelet-derived growth factor receptor alpha⁺ fibro/adipogenic progenitors would differentiate into fibroblasts in fatty infiltrated rotator cuff muscles.²⁰ Monocytes and mast cells both participated in the process of RCT. Krieger et al showed that multiple immune cells such as macrophages, monocytes and T cells were accumulated in rotator cuff muscles, while exhaustion of circulating monocytes expedited rotator cuff muscle atrophy after RCT.⁷¹ Sun et al found that preoperative lymphocyte to monocyte ratio less than 4.760 was a poor prognostic factor for arthroscopic rotator cuff repair.⁷² Moreover, circulating monocytes could recruit and switch into macrophages after skeletal muscle injury,⁷³ while macrophages would infiltrate in fatty degenerative rotator cuff muscles.^21–23 Mast cells could boost myoblast proliferation in injured skeletal muscle.⁷⁴ Matthews et al and Millar et al both reported that mast cells recruited in RCT, but the number of mast cells were lower in large and massive RCT than small and medium RCT.^75,76 However, in the present study, decreased infiltrations of monocytes and mast cells were observed in fatty infiltrated rotator cuff muscle. These may attribute to the own characteristics of stromal and immune cell infiltration in FI.

Correlation analysis between FABP5, MGP and immune and stromal cells showed that FABP5 and MGP were positively correlated with vessels, and negatively correlated with monocytes and mast cells. FABP5 could boost angiogenic responses in endothelial cells of blood vessels, and might regulate inflammation, lipid raft constitution of monocytes and metabolism in activated monocytes.^77,78 Moreover, FABP5 took part in the lipopolysaccharide-induced cytokine production of mast cells.⁷⁹ MGP was expressed in blood vessels, and would be up-regulated in lipopolysaccharide-stimulated THP-1 monocytes.^80,81 Thus, FABP5 and MGP may regulate vessels, monocytes and mast cells in the progress of FI.

The highlight of the present study was to employ two machine learning algorithms (random forest and mSVM-RFE) for date mining and adopt mMCP-counter method for the first time to analyze stromal and immune cell infiltration in FI. The present study has several limitations. This study was a bioinformatics research. The data was downloaded from GEO database for second mining, and further experiment should be conducted to validate the results. In addition, the sample size in the GSE103266 validation cohort was small.

Conclusions

In summary, FABP5 and MGP were identified as diagnostic biomarkers of FI after RCT. Vessels, endothelial cells, fibroblasts, monocytes and mast cells may be related to the process of FI. In addition, FABP5 and MGP were significantly associated with vessels, monocytes and mast cells. These immune and stromal cells may be potential targets of FI immunotherapy.