Online ISSN: 1945-4589
Learn about our FREE Post-Publication Promotion Services
Research Paper Volume 14, Issue 10 pp 4530—4555

Identification of hub biomarkers and immune cell infiltration in polymyositis and dermatomyositis

Si Chen1,2, , Haolong Li1, , Haoting Zhan1, , Xiaoli Zeng2, , Hui Yuan2, , Yongzhe Li1, ,

  • 1 Department of Clinical Laboratory, State Key Laboratory of Complex Severe and Rare Diseases, Peking Union Medical College Hospital, Chinese Academy of Medical Science and Peking Union Medical College, Beijing, China
  • 2 Department of Clinical Laboratory, Beijing Anzhen Hospital, Capital Medical University, Beijing, China

Received: December 23, 2021       Accepted: April 12, 2022       Published: May 24, 2022      

https://doi.org/10.18632/aging.204098
How to Cite

Copyright: © 2022 Chen et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Objective: Polymyositis (PM) and dermatomyositis (DM) are heterogeneous disorders. However, the etiology of PM/DM development has not been thoroughly clarified.

Methods: Gene expression data of PM/DM were obtained from Gene Expression Omnibus. We used robust rank aggregation (RRA) to identify differentially expressed genes (DEGs). Gene Ontology functional enrichment and pathway analyses were used to investigate potential functions of the DEGs. Weighted gene co-expression network analysis (WGCNA) was used to establish a gene co-expression network. CIBERSORT was utilized to analyze the pattern of immune cell infiltration in PM/DM. Protein–protein interaction (PPI) network, Venn, and association analyses between core genes and muscle injury were performed to identify hub genes. Receiver operating characteristic analyses were executed to investigate the value of hub genes in the diagnosis of PM/DM, and the results were verified using the microarray dataset GSE48280.

Results: Five datasets were included. The RRA integrated analysis identified 82 significant DEGs. Functional enrichment analysis revealed that immune function and the interferon signaling pathway were enriched in PM/DM. WGCNA outcomes identified MEblue and MEturquoise as key target modules in PM/DM. Immune cell infiltration analysis revealed greater macrophage infiltration and lower regulatory T-cell infiltration in PM/DM patients than in healthy controls. PPI network, Venn, and association analyses of muscle injury identified five putative hub genes: TRIM22, IFI6, IFITM1, IFI35, and IRF9.

Conclusions: Our bioinformatics analysis identified new genetic biomarkers of the pathogenesis of PM/DM. We demonstrated that immune cell infiltration plays a pivotal part in the occurrence of PM/DM.

Introduction

The idiopathic inflammatory myopathies (IIMs) are characterized by muscle injury, with an annual incidence of 7.98 per million people, according to data from 1966 to 2013 [1]. IIMs include polymyositis (PM), dermatomyositis (DM), inclusion body myositis (IBD), juvenile dermatomyositis (JDM), anti-synthetase syndrome (ASS), and immune-mediated necrotizing myopathies (IMNM); historically, the two most common IIMs are PM and DM [2]. The disease progression and muscle impairments of PM and DM are similar, irrespective of the presence of skin lesions [3]. In addition, clinical studies have shown similar serological characteristics between PM and DM, especially with certain myositis-specific autoantibodies [4, 5]. Despite several decades of research, the exact pathogenesis of PM and DM remains unknown; their etiology may be a combination of multiple factors, such as immune activation, genetic background, and environmental factors [6].

In general, elucidating specific molecular pathways associated with a given disease would be of great significance in identifying disease subgroups, monitoring disease activity, and selecting treatment approaches. Studies of PM/DM-associated genes in the immune system and inflammatory response signaling pathways have also been reported, such as interferon-alpha (IFN-α) [7, 8], nuclear factor-κB (NF-κB) [9], IFN-γ [7, 8], tumor necrosis factor α (TNF-α) [10], toll-like receptors (TLRs) [11], and retinoic acid-inducible gene 1 (RIG-1) [12]. According to the literature, gene microarray technology has been extensively utilized to analyze gene expression in muscle or skin tissues of PM/DM patients. However, there are some contradictions regarding these microarray data [13–17] in reflecting the importance of IFN1-induced genes in the pathogenesis of DM and PM. These contradictions may be attributed to factors such as diverse analytic methods and platforms, inconsistent specimen quantity, and different sample sources.

Bioinformatics analysis is an efficient means of deep detection and mining of transcriptomic data. Furthermore, the robust rank aggregation (RRA) method has been utilized to discriminate variation in mRNA profiles between several datasets in diverse disease types, such as in tumors [18–20]. In this study, we included five mRNA microarray datasets from Gene Expression Omnibus (GEO) and used the RRA method to select differentially expressed genes (DEGs) between PM/DM and healthy controls. We next analyzed the molecular mechanisms of PM/DM using Gene Ontology (GO) functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. In addition, other enrichment analyses such as gene set enrichment analysis (GSEA) and gene set variation analysis (GSVA) were used to study the molecular pathways of PM/DM. Weighted gene co-expression network analysis (WGCNA) was utilized to create a gene co-expression network and identify the most relevant modules in PM/DM. CIBERSORT analysis was utilized to estimate the various immune infiltrates in PM/DM. Protein–protein interaction (PPI) network analysis was performed to screen key genes, and Venn diagram analysis was used to assemble the DEGs of seven datasets. Subsequently, analysis of the correlation between core genes and muscle injury was conducted to indicate the potential functions of core genes in PM/DM. Finally, verification testing was carried out to identify five novel hub genes in the pathogenesis of PM/DM. We expect that the newly described DEGs and altered pathways between PM/DM and healthy controls found in this study may help to elucidate possible molecular mechanisms for their pathogenesis.

Materials and Methods

Study design and data collection

Gene expression datasets were filtered using the GEO (http://www.ncbi.nlm.nih.gov/geo) database [21]. We comprehensively retrieved microarray studies using the keywords “Polymyositis,” “Dermatomyositis,” “Myositis,” “Gene expression,” “Homo sapiens,” and “Microarray.” Datasets were filtered using the next enrollment criteria: (1) including more than 10 specimens; (2) original information or gene expression analysis data obtainable in GEO. In light of the above criteria, GSE1551 [13], GSE3112 [22], GSE39454 [23], GSE46239, and GSE128470 [24] were finally recruited. Subsequently, the normalization and quality control of these data were executed with the “affy” R package [25]. We used the "sva" package [26] to eliminate the batch effect. The probes were transformed into homologous gene symbols using the R package “Rsubread” [27]. If several probes matched an identical sign, their average value was obtained. Further, genes without a corresponding genetic symbol were deleted.

DEG screening

We established seven different groups for five GEO series: the GSE39454 series, which was split into GSE39454PM and GSE39454DM series; GSE128470, which was divided into GSE128470 PM and GSE128470 DM series. Data were analyzed using the “limma” (linear models for microarray data) R package [28] to detect all DEGs between PM/DM and healthy controls. The values for statistical significance were set as adjusted P≤0.05 and |log2 fold change (FC)|≥1, except for GSE39454 PM and GSE39454 DM, which were set as P≤0.05 and |log2FC|≥1. Volcano maps were drawn using the “ggplot2” [29] package. Principal component analysis (PCA) was used to extract two features from each group of genes. PCA score trajectory plots can be used to show whether there is overlap between PM/DM and the control group. When there was no substantial overlap, it suggested that there was a significant difference between PM/DM and the control group, which could be analyzed in the next step.

RRA analysis

To reduce the discrepancies and merge the outcomes of multiple microarray studies, we conducted RRA analysis to recognize substantial DEGs. RRA is an efficient tool to combine results from several arrays [30]. First, we acquired lists of up-ranked and down-ranked genes of each series that were produced by analyzing the FC of expression between PM/DM and controls. The “Robust Rank Aggregation” R package was utilized to aggregate all sorted gene lists in the datasets [30]. The Benjamini and Hochberg false discovery rate (FDR) method was utilized to generate the adjusted P-value ranked genes with adjusted P<0.05 and log2FC>0.5 were regarded as significant (set 1).

Functional and pathway enrichment analysis

To examine the effect of DEGs on the pathogenesis of PM/DM, we executed GO functional enrichment analysis and KEGG pathway analysis of the significant genes identified by RRA. The “clusterProfiler” R package automates the classification of biological terminology and gene cluster enrichment analysis. The analysis and visualization modules were amalgamated into a repeated workflow [31]. Furthermore, we used the “GOplot” R package (circle plot, chord plot, cluster plot) to append quantitative data about pathways to the GO terms of interest [32]. An adjusted P<0.05 and FDR <0.05 were considered the criteria for judgment.

Enrichment analysis by GSEA and GSVA

To identify the possible functional pathways associated with DEGs in PM/DM, we executed GSEA to investigate the biological processes and pathways of these genes and used the GSEA plot function of the “clusterProfiler” R package to carry out GSEA [31]. The annotated gene set “h.all.v7.0.symbols.gmt” in the MsigDB V6.2 database [33] was regarded as a reference gene set. FDR<0.25, normalized P<0.05 and |normalized enrichment score (NES)|>1 were deemed to show meaningful enrichment. Further, we used the gseaplot2 function of the “clusterProfiler” R package to visualize the results of GSEA [31]. The R package “GSVA” was utilized to calculate the signaling pathways in enrolled datasets [34]. Subsequently, the R package “limma” was used to research the meaningful distinguishing gene sets between PM/DM and healthy controls [28]. The gene set “h.all.v7.0.symbols.gmt” was chosen as the reference gene set. Volcano maps were drawn using the “ggplot2” [29] package. The “pheatmap” package was used to generate the heatmap plots to visualize the results of GSVA analysis [35]. Signaling pathways that were usually enriched by GSEA and GSVA analysis were regarded as possible PM/DM-related pathways.

WGCNA

We used the “WGCNA” R package to make a co-expression network on the basis of the above DEGs. We built a weighted adjoining matrix, identified outliers by sample clustering, and removed the outliers. In addition, we defined the association degree (soft threshold parameter) to indicate a strong intergenic association while carrying out scale-free network verification. Further, we converted the adjoining matrix to a topological overlap matrix (TOM) to estimate the degree of connectedness between genes. To avoid bias and error, the minimum number of genes in each module was set to 30. Using the average linkage hierarchical clustering method on account of the TOM discrepancy measure, we categorized the genes with similar expression spectra using gene modules, displaying them on branches and with different colors of clustering trees to display the relationships between modules. The correlation cutoff was set at 0.8. To verify the interactions occurring significantly above the expected probability due to chance, a control network analysis was also performed using a Z-test, and q-value analysis was also performed to the correction analysis of multiple tests. The relationship between gene module and phenotype was computed and the modules related to clinical characteristics were selected.

Assessment of immune cell infiltration by CIBERSORT

Here, we adopted the CIBERSORT method to analyze the expression of 22 immune cell categories across seven gene expression matrix datasets [36]. We selected the characteristic gene matrix file LM22, which is a leukocyte gene characteristic matrix for identifying 22 human hemopoietic cell phenotypes [36]. We then set the run mode to batch mode, utilized the quantile standardization of expression matrix, and set 500 permutations for penetration analysis. Results with P<0.05 were considered statistically significant. The Mann-Whitney U test was utilized to discover the meaningfully distinct infiltrating immune cell sorts between PM/DM and healthy controls. In addition, we used “ggplot2” package [29] to design boxplots to show the differentiation in immune cell infiltration.

PPI network analysis

The STRING database is an online database for retrieving known proteins and forecasting PPI, including the direct physical interrelations between proteins and their connections to indirect functions [37]. We provided the meaningful genes from the above-mentioned RRA analysis to the STRING database to construct a PPI network (set 2). We set the medium confidence score as>0.4 when outputting the network analysis results and exported the TSV format data to Cytoscape3.7.2 [38]. The CytoHubba (version 0.1) plug-in for Cytoscape was utilized to distinguish the association degrees of genes in the PPI network and visualize the network. The Molecular Complex Detection (MCODE version 1.6.1) software (http://apps.cytoscape.org/apps/mcode) was adopted to choose the key modules from the PPI network in Cytoscape with MCODE scores >5 (set 3 and set 4).

Venn diagram analysis

Venn diagram analysis was carried out using the Venn Diagram R package (version 2.12.1) [39] for the DEGs of seven datasets. Overlapping DEGs (set 5) were maintained for further analysis. The comprehensive analysis of sets 1- 5 indicated the core genes. The core genes were used for the next analysis.

Correlation between core genes and muscle injury

We searched for muscle injury-associated genes using the GeneCards website (https://www.genecards.org/) with the term “muscle injury.” A relevance score based on a scale of 0 to 100 was used to demonstrate the strength of the association between genes and muscle injury. The higher the score, the more relevant it was. The relevant scores were sorted in descending order, and the top 16 genes were considered the uppermost muscle injury-related genes. Furthermore, we used the “corrplot” package [40] to calculate the correlation between core genes and muscle injury-related genes. In addition, we used “ggcorrplot” package [41] to visualize the correlation. The final core genes were determined according to the correlation score.

Diagnostic effectiveness evaluation

For diagnostic analysis, we selected GSE39454 and GSE128470, which contain both PM and DM samples. In addition, we carried out confirmation studies using data from a microarray dataset (GSE48280, including 5 DM patients, 5 PM patients, and five healthy controls) [42]. The receiver operator characteristic (ROC) curves were diagramed and area under curve (AUC) was counted respectively to appraise the performance of each dataset (GSE39454 PM, GSE39454 DM, GSE128470 PM, GSE128470 DM, GSE48280 PM and GSE48280 DM) utilizing the R packages “pROC” [43]. We defined a guideline to distinguish different diagnostic criteria, including excellent accuracy (0.9≤AUC<1), decent accuracy (0.8≤AUC<0.9), fair accuracy (0.7≤AUC<0.8), poor accuracy (0.6≤AUC<0.7), and insufficient accuracy (0.5≤AUC<0.6) [44]. If the AUC value of a hub gene was >0.8, it was regarded as having excellent specificity and sensitivity for identifying PM/DM.

Results

Essential information of selected microarrays

In light of the above inclusion criteria, GSE1551, GSE3112, GSE39454, GSE46239, and GSE128470 were finally selected. A total of 145 samples (81 DM samples, 22 PM samples, and 42 control samples) were evaluated in our study. Analyses of GSE1551, GSE3112, and GSE128470 series were undertaken on the GPL96 platform (Affymetrix Human Genome U133A Array), and GSE39454 and GSE46239 were performed on the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array). The details of these datasets are presented in Table 1. Our analysis of microarray data is based on the basic workflow (Figure 1).

Table 1. Characteristics of the enrolled microarray datasets.

GSE IDDMPMControlTissuesAnalysis typePlatformYear
GSE15511310MuscleArrayGPL962005
GSE3112611MuscleArrayGPL962005
GSE39454895Skeletal muscleArrayGPL5702012
GSE46239484SkinArrayGPL5702013
GSE12847012712MuscleArrayGPL962019
Data analysis and processing flow. The data processing in this study is divided into five steps.

Figure 1. Data analysis and processing flow. The data processing in this study is divided into five steps.

Evaluation of DEGs in PM/DM

First, to exclude individual variations between participants, all five microarray datasets were normalized using the “affy” R package. The standardized boxplots are presented in Supplementary Figure 1; all participants in each dataset attained acceptable homogeneity. Secondly, the PCA plots of all data series are shown in Supplementary Figure 2. In view of the gene expression of all participants, PCA showed differing distribution patterns between PM/DM and control groups. The ranges between the participants in the control group were similar, as were the ranges between the participants in the DM or PM groups. Furthermore, we utilized the “limma” R package to screen out the DEGs on account of the above screening criteria, and the volcano plots of the seven microarrays are displayed in Figure 2.

Volcano diagrams of the seven microarrays and heatmap of the robust rank aggregation (RRA) analysis. Red points indicate upregulated genes, while blue points indicate downregulated genes. Gray points indicate genes with no meaningful difference. (A) GSE1551 (dermatomyositis, DM); (B) GSE3112 (polymyositis, PM); (C) GSE39454 (DM); (D) GSE39454 (PM); (E) GSE46239 (DM); (F) GSE128470 (DM); (G) GSE128470 (PM); (H) Heatmap of the top 20 upregulated and five downregulated genes in the RRA method. Red and blue indicate high and low expression of genes in patients with PM/DM, respectively.

Figure 2. Volcano diagrams of the seven microarrays and heatmap of the robust rank aggregation (RRA) analysis. Red points indicate upregulated genes, while blue points indicate downregulated genes. Gray points indicate genes with no meaningful difference. (A) GSE1551 (dermatomyositis, DM); (B) GSE3112 (polymyositis, PM); (C) GSE39454 (DM); (D) GSE39454 (PM); (E) GSE46239 (DM); (F) GSE128470 (DM); (G) GSE128470 (PM); (H) Heatmap of the top 20 upregulated and five downregulated genes in the RRA method. Red and blue indicate high and low expression of genes in patients with PM/DM, respectively.

RRA integrated analysis of DEGs

The RRA analysis posits that each gene is irregularly arranged in each dataset. The lower the P-value in RRA outcomes, the higher the gene grade and the reliability of distinguishing gene expression. After the integrated analysis, 82 significant DEGs (70 upregulated and 12 downregulated) were identified (Supplementary Table 1; set 1). The heatmap of the top 20 upregulated and five downregulated genes is displayed in Figure 2H. The top 10 meaningfully upregulated genes found in PM/DM included ISG15 (P = 4.77E-08), C1QB (P = 7.73E-08), HLA-A (P = 9.33E-07), HLA-C (P = 1.53E-06), HLA-B (P = 2.15E-06), PSMB8 (P = 7.57E-06), C1QA (P = 9.28E-06), GBP1 (P = 1.05E-05), UBE2L6 (P = 1.07E-05), PARP12 (P = 3.04E-05).

Functional annotation

According to the above experimental results, 82 DEGs were subjected to GO [including three main function modules: biological process (BP), molecular function (MF), and cellular component (CC)] and KEGG analyses. The analyses demonstrated that the type I interferon signaling pathway (GO: 0060337; adjusted P = 2.14E-31) was the most meaningfully enriched BF, followed by cellular response to type I interferon (GO: 0071357; adjusted P = 2.14E-31) and response to type I interferon (GO: 0034340; adjusted P = 3.93E-31). Different analysis methods showed differing results; thus, we used the GO bar plot (Figure 3A), cluster plot (Figure 3B), chord plot (Figure 2C), circle plot (Figure 3D), and cneplot (Figure 3E) to visualize the DEGs and GO terms. Furthermore, KEGG pathway enrichment analysis showed that Epstein-Barr virus infection (hsa05169; adjusted P = 7.91E-09) and coronavirus disease - COVID-19 (hsa05171; adjusted P = 2.98E-07) were significantly enriched (Figure 3F).

Gene ontology (GO) functional enrichment analysis and Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis of differentially expressed genes. (A) GO bar plot (Figure 3A); (B) GO cluster plot; (C) GO chord plot; (D) GO circle plot; (E) GO cneplot; (F) KEGG bar plot.

Figure 3. Gene ontology (GO) functional enrichment analysis and Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis of differentially expressed genes. (A) GO bar plot (Figure 3A); (B) GO cluster plot; (C) GO chord plot; (D) GO circle plot; (E) GO cneplot; (F) KEGG bar plot.

GSEA and GSVA

We comprehensively analyzed the results of all GSEA and GSVA datasets. We screened out four commonly enriched pathways: allograft rejection, inflammatory response, interferon-alpha response, and interferon-gamma response (Figures 4–6A and Supplementary Figure 3). Analyses of the present research revealed that interferon response and inflammatory response were the biological pathways most relevant to the pathogenesis of PM/DM.

Gene set enrichment analysis (GSEA) results of six microarrays. Normalized enrichment score (NES) demonstrates the analysis outcomes across gene sets. False discovery rate (FDR) indicates if a set was meaningfully enriched. (A) GSE1551 (dermatomyositis, DM); (B) GSE3112 (polymyositis, PM); (C) GSE39454 (DM); (D) GSE46239 (DM); (E) GSE128470 (DM); (F) GSE128470 (PM).

Figure 4. Gene set enrichment analysis (GSEA) results of six microarrays. Normalized enrichment score (NES) demonstrates the analysis outcomes across gene sets. False discovery rate (FDR) indicates if a set was meaningfully enriched. (A) GSE1551 (dermatomyositis, DM); (B) GSE3112 (polymyositis, PM); (C) GSE39454 (DM); (D) GSE46239 (DM); (E) GSE128470 (DM); (F) GSE128470 (PM).

Gene set variation analysis (GSVA) results of six microarrays. (A) GSE1551 (dermatomyositis, DM); (B) GSE3112 (polymyositis, PM); (C) GSE39454 (DM); (D) GSE39454 (PM); (E) GSE46239 (DM); (F) GSE128470 (DM).

Figure 5. Gene set variation analysis (GSVA) results of six microarrays. (A) GSE1551 (dermatomyositis, DM); (B) GSE3112 (polymyositis, PM); (C) GSE39454 (DM); (D) GSE39454 (PM); (E) GSE46239 (DM); (F) GSE128470 (DM).

Gene set variation analysis (GSVA) of GSE128470 (polymyositis, PM) and weighted gene co-expression network analysis (WGCNA) results of five datasets. (A) GSVA result of GSE128470 (PM); (B) WGCNA result of GSE1551; (C) WGCNA result of GSE3112; (D) WGCNA result of GSE39454; (E) WGCNA result of GSE46239; (F) WGCNA result of GSE128470.

Figure 6. Gene set variation analysis (GSVA) of GSE128470 (polymyositis, PM) and weighted gene co-expression network analysis (WGCNA) results of five datasets. (A) GSVA result of GSE128470 (PM); (B) WGCNA result of GSE1551; (C) WGCNA result of GSE3112; (D) WGCNA result of GSE39454; (E) WGCNA result of GSE46239; (F) WGCNA result of GSE128470.

WGCNA

In this study, the WGCNA R package was utilized for co-expression network construction. Before analysis, we determined suitable soft thresholds to establish a scale-free network (Supplementary Figure 4). To ensure the scale-free topology model curve fit the plateau smoothly, different soft thresholds were selected for each dataset (Supplementary Figure 4). Hereafter, the stepwise method was used for dynamic cluster analysis. First, dynamic hybrid cutting was used to generate the dendrogram (Supplementary Figure 5). Each leaf represents a single gene that has a close expression profile when placed together to form a branch, and the branch represented a gene module. Subsequently, we also computed the characteristic genes of each module and clustered the modules in parallel, especially when the correlation was >0.8; then, these modules were merged. Thus, we obtained significant modules excluding the nonsense one (gray; Supplementary Figure 5). The adjacency heatmap of characteristic genes showed that the red and blue modules were the most positively and negatively correlated with the occurrence of PM/DM, respectively (Figure 6B–6F). Taken together, these results indicated that MEblue and MEturquoise are key target modules in PM/DM.

Immune cell infiltration outcomes

The boxplots of the immune cell infiltration results demonstrated that compared with those in the healthy control participants (immune cell infiltration results of healthy controls are shown in Supplementary Figure 6), M1 (P = 0.03) and M2 (P<0.0001) macrophage infiltration in GSE1551 was higher and that of regulatory T cells (Tregs) (P< 0.0001) in GSE1551 was less (Figure 7A). M2 macrophages (P< 0.0001) in GSE3112 infiltrated more, and memory B cells (P = 0.007), eosinophils (P = 0.003) and follicular helper T cells (P = 0.015) in GSE3112 infiltrated less (Figure 7B). M1 macrophages (P = 0.002) in GSE39454 (DM) infiltrated more (Figure 7C), and M1 macrophages (P = 0.019) in GSE39454 (PM) infiltrated more (Figure 7D). M1 macrophages (P = 0.01) and M2 macrophages (P = 0.001) in GSE46239 infiltrated more and resting dendritic cells (P = 0.001) and Tregs (P = 0.041) in GSE46239 infiltrated less (Figure 7E). M0 macrophages (P< 0.0001) and M1 macrophages (P = 0.001) in GSE128470 (DM) infiltrated more, and CD8 T cells (P = 0.004) in GSE46239 infiltrated less (Figure 7F). M1 macrophages (P = 0.001) in GSE128470 (PM) infiltrated more, and neutrophils (P = 0.003), plasma cells (P = 0.01) and Tregs (P = 0.013) in GSE128470 (PM) infiltrated less (Figure 7G). Overall, our results indicated that M0, M1, and M2 macrophages in PM/DM patients infiltrated more, and Tregs in PM/DM patients infiltrated less.

Boxplots of the proportion of 22 immune cell sorts in polymyositis (PM) and dermatomyositis (DM) and the outcomes of protein–protein interaction (PPI) network analysis. (A) GSE1551 (DM); (B) GSE3112 (PM); (C) GSE39454 (DM); (D) GSE39454 (PM); (E) GSE46239 (DM); (F) GSE128470 (DM); (G) GSE128470 (PM); (H) Results of Cytoscape plug-in CytoHubba; (I) the top significant module of Cytoscape plug-in MCODE.

Figure 7. Boxplots of the proportion of 22 immune cell sorts in polymyositis (PM) and dermatomyositis (DM) and the outcomes of protein–protein interaction (PPI) network analysis. (A) GSE1551 (DM); (B) GSE3112 (PM); (C) GSE39454 (DM); (D) GSE39454 (PM); (E) GSE46239 (DM); (F) GSE128470 (DM); (G) GSE128470 (PM); (H) Results of Cytoscape plug-in CytoHubba; (I) the top significant module of Cytoscape plug-in MCODE.

PPI network analysis

Using the STRING online database, we input significant genes from the RRA method for PPI network analysis, and we used Cytoscape software to show the results. In the PPI network, the connectivity between nodes was highlighted to determine the interrelations between the proteins encoded by genes in PM/DM, and a PPI network of 65 object genes was established. The genes situated in the centric node were recognized as key genes that may play crucial physiological regulatory actions in PM/DM (Figure 7H, orange node genes) (set 2). Among the 65 nodes, 16 proteins were picked by degree through Cytoscape CytoHubba plug-in; namely, central genes (ISG15, MX1, STAT1, CXCL10, OAS2, TRIM22, IFI6, IFITM1, IFI35, IFI27, IRF9, GBP2, OAS1, RSAD2, GBP1, and IFIT3) (Figure 7H, orange node genes) (set 2). In addition, another Cytoscape plug-in, MCODE, was used to further categorize the PPI to recognize sub-network. Here, we utilized MCODE for cluster analysis. The top meaningful module was selected after clustering (Figure 7I, orange node genes) (set 3). Cluster 1 included 28 genes and the 12 core genes in the central node contained ISG15, MX1, IFI6, IFITM1, IFI35, IFI27, IRF9, OAS1, OAS2, RSAD2, GBP1, and IFIT3. The top two significant module was obtained and included nine genes (C1QB, VSIG4, CD163, CCL2, MS4A4A, CCL8, CCR1, TYROBP, and C1QA) (Figure 8A) (set 4).

Protein–protein interaction (PPI) network analysis, Venn diagram analysis and correlation between core genes and muscle injury. (A) The top two significant modules of Cytoscape plug-in MCODE. (B) Results of Venn diagram analysis; (C) Association result of GSE1551 (DM); (D) Association result of GSE3112 (PM); (E) Association result of GSE39454 (DM); (F) Association result of GSE39454 (PM); (G) Association result of GSE46239 (DM); (H) Association result of GSE128470 (DM); (I) Association result of GSE128470 (PM).

Figure 8. Protein–protein interaction (PPI) network analysis, Venn diagram analysis and correlation between core genes and muscle injury. (A) The top two significant modules of Cytoscape plug-in MCODE. (B) Results of Venn diagram analysis; (C) Association result of GSE1551 (DM); (D) Association result of GSE3112 (PM); (E) Association result of GSE39454 (DM); (F) Association result of GSE39454 (PM); (G) Association result of GSE46239 (DM); (H) Association result of GSE128470 (DM); (I) Association result of GSE128470 (PM).

Venn diagram analysis and identification of core genes

We used the Venn diagram analysis for DEGs of seven datasets and identified nine overlapping genes (Figure 8B) (set 5), including ISG15, IFIT3, IRF9, LY96, C1R, CD163, C1QB, C1QA, and LGALS3BP (Figure 8B). The comprehensive analysis of set 1, set 2, set 3, set 4, and set 5 recognized core genes. The core genes, including, IFIT3, ISG15, MX1, STAT1, CXCL10, OAS2, TRIM22, IFI6, IFITM1, IFI35, IFI27, IRF9, GBP2, OAS1, RSAD2, and GBP1 were selected for further analysis.

Correlation between core genes and muscle injury

We searched genes associated with the muscle injury in the GeneCards database. The muscle injury-related genes in the top 16 of relevance scores contained RYR1, TTN, DMD, IL6, TNF, CAV3, DES, MSTN, DYSF, INS, MYOD1, ACTA1, MYH7, TGFB1, SOD1, and SCN4A. The correlation between muscle injury-related genes in the top 16 of relevance scores and above core genes was analyzed (Figure 8). According to the results of correlation analysis, we identified five final hub genes (TRIM22, IFI6, IFITM1, IFI35, and IRF9) associated with PM/DM, which were used for subsequent investigation. Hub genes are generally regarded as key functional genes and are highly associated with other genes.

Diagnostic value of the five hub genes

To validate the diagnostic value of five hub genes (TRIM22, IFI6, IFITM1, IFI35, and IRF9) in PM/DM patients, the next phase of our study was to execute ROC analyses to investigate the sensitivity and specificity of hub genes for PM/DM diagnosis. The ROC outcomes verified that five hub genes could differentiate between PM/DM patients and healthy controls in GSE39454 DM, GSE39454 PM, GSE48280 DM and GSE128470 DM (all P<0.05), and the AUCs were between 0.799 and 1 (Figure 9A–9D). However, the diagnostic value of the five hub genes (especially IFI35 and IFI6) in GSE1551DM, GSE3112PM, and GSE46239 DM was uncertain (Supplementary Figure 7A–7E). This result may be attributed to the large difference in sample size between DM patients and healthy controls in GSE46239 and the small sample sizes of GSE1551 and GSE3112 that lead to some bias. Our results indicated that expression of TRIM22, IFI6, IFITM1, IFI35, and IRF9 was related to disease diagnosis and that these genes could act as biomarkers to verify the diagnosis of PM/DM and validate the effectiveness of PM/DM treatment.

Receiver operating characteristics of five hub genes. (A) GSE39454 (DM); (B) GSE39454 (PM); (C) GSE48280 (DM); (D) GSE128470 (DM).

Figure 9. Receiver operating characteristics of five hub genes. (A) GSE39454 (DM); (B) GSE39454 (PM); (C) GSE48280 (DM); (D) GSE128470 (DM).

Discussion

PM and DM are common autoimmune disorders clinically presenting with skeletal muscle injury [45]. At present, there is no specific diagnostic antibody for PM and DM. The diagnosis mainly depends on medical experience and invasive muscle biopsy [45]. Hence, there is an imperative need to better understand the pathogenesis of PM/DM to generate novel strategies for the diagnosis and treatment of PM/DM.

In the current research, a large number of bioinformatic analysis tools have been used to identify five hub genes (TRIM22, IFI6, IFITM1, IFI35, and IRF9) between PM/DM and health control subjects on the basis of gene expression profiles attained from GSE1551 (DM), GSE3112 (PM), GSE39454 (DM), GSE39454 (PM), GSE46239 (DM), GSE128470 (DM), and GSE128470 (PM) datasets. We explored the biological functions of these DEGs utilizing GO enrichment and KEGG pathway analyses; these results demonstrated that DEGs are significantly associated with the IFN response pathway. Moreover, deep enrichment analysis by GSEA and GSVA revealed that immune function and the IFN signaling pathway are pivotal features implicated in PM/DM, which was in accordance with the findings of previous studies [13–17, 46]. Most notably, WGCNA was conducted to identify meaningful modules correlated with PM/DM, indicating that MEblue and MEturquoise could be identified as key target modules in PM/DM. Furthermore, this study is the first to use a bioinformatics analysis with an immune cell infiltration analysis in PM/DM. The results indicated that M0 macrophages, M1 macrophages, and M2 macrophages in PM/DM patients infiltrated more, and Tregs in PM/DM patients infiltrated less. We performed PPI network analysis and Venn diagram analysis, and we assessed the associations between core genes and muscle injury to identify hub genes. Next, we performed ROC analyses to investigate the sensitivity and specificity of five hub genes for the diagnosis of PM/DM, and results indicated that expression of TRIM22, IFI6, IFITM1, IFI35, and IRF9 are related to the diagnosis of PM/DM.

The five hub genes (TRIM22, IFI6, IFITM1, IFI35, and IRF9) are all IFN pathway genes. In all signaling pathways of PM/DM, the one centered on IFNs has been the most studied and IFNs have been identified as playing vital roles in PM/DM. Initial investigation of cytokine expression revealed the upregulation of IFN-γ in PM/DM muscle [47], which leads to the localized overexpression of IFN-γ-related genes [48]. These findings suggested that transcriptomic changes are involved in the pathogenesis of PM/DM. According to the cell surface receptor binding ligand family, the IFN pathway is divided into three categories: type 1 IFNs (IFN1; including IFN-α and IFN-β), type 2 IFNs (IFN2 or IFN-γ), and type 3 IFNs (IFN3 or IFN-λ), which share overlapping signaling pathways [49, 50]. Analysis has suggested that IFN1 is potentially relevant to the pathogenesis of DM [13, 14]. In particular, overexpression of IFN1-related genes has been described in the muscle [13], skin [14], and circulating leukocytes [15, 16] of DM patients. Furthermore, the expression of IFN1-induced genes is associated with the disease activity of DM [15, 16]. In addition, a previous study on the skin of DM patients demonstrated that both IFN-β and IFN-γ are highly correlated with the degree of the IFN transcriptional response, whereas IFN-α are not [14]. However, related investigation on the IFN pathway in the pathogenesis of PM is lacking. Although both PM and DM were implicated in some studies, the IFN pathway was found to be associated with DM but not with PM [13, 17]. Nevertheless, in other studies, PM and DM samples showed similar gene expression profiles, although there were differences in several other details [15, 46]. Thus, the role of the IFN pathway in the pathogenesis of PM/DM remains controversial.

TRIM22 encodes an IFN-induced member of the tripartite motif (TRIM) family, which is localized to the cytoplasm. Previous studies indicated that this protein might participate in the antiviral effects of IFN [51]. Previous reports on TRIMs in autoimmunity have been inconclusive. TRIMs, as active modulators of IFN production and inflammasome activity, may counteract autoimmune or inflammatory disorders if their activity or expression is improved. Conversely, TRIMs that negatively regulate IFN production or inflammasome activity may be beneficial in some diseases, as they can prevent the excessive production of pathogenic cytokines, thus avoiding the progression of autoimmune and autoinflammatory diseases [52]. In addition, the researchers found that variants in TRIM22 influenced the activation of NOD2-dependent IFN-β signaling and the upregulation of NF-κB pathways in early-onset inflammatory bowel disease, and the TRIM22-NOD2 network affected antiviral pathways leading to inflammation [53]. Meanwhile, only one early-onset Crohn’s disease (CD) patient had a TRIM22 variant that was consistent with the previous result [54]. These results suggested the TRIM22 variant affects pediatric patients with inflammatory bowel disease. Another study showed that psoriasis was associated with the overexpression of antiviral genes (such as ISG15 and TRIM22) in the skin, but not in the blood [55]. However, the role of TRIM22 in the development of multiple sclerosis (MS) has been debatable [56–58]. Therefore, the role of TRIM22 in connective tissue disorders needs further elucidation and interpretation.

IFI6, also known as G1P3, is one of the IFN-stimulating genes (ISGs). IFI6 may be crucial in mediating apoptosis in many cancers [59], and it may have antiviral activity against the hepatitis C virus [60]. Furthermore, the expression of IFI6 is elevated in the blood or platelet samples of patients with systemic lupus erythematosus (SLE), rheumatoid arthritis (RA), primary antiphospholipid syndrome, and MS [61–63]. Meanwhile, bioinformatic analysis of psoriasis and SLE cases suggested IFI6 is a primary IFN-inducible gene [64, 65]. IFITM1, a member of the IFN-induced transmembrane protein family, is a known regulator of immunity and antiviral activity [66]. It is interesting to note that mounting evidence has shown that IFITM1 is highly expressed in many tumor tissues [67]. However, the definite mechanism of IFITM1 in autoimmune diseases remains unclear. IFI35 encodes IFN-induced 35 kDa protein, an IFN-induced protein, which mainly participates in, and regulates the response of, the innate immune system [68–70]. IFI35 forms complexes with N-myc and STAT interactor (NMI), regulating various immune responses, such as restricting virus-triggered IFN-β production [70] and passively adjusting NF-κB signaling, leading to the restraint of endothelial cell proliferation and migration [70]. IFI35 plays different roles in connective tissue disease. For example, IFI35 is deemed a biomarker of neuroinflammation and therapeutic reaction in MS [71]. Additionally, IFI35 promotes the proliferation of mesangial cells, which is modulated by methyl CpG-binding domains in lupus nephritis [72].

IRF9 encodes a transcription factor that plays a principal role in antiviral immunity. It participates in the IFN response and modulates cell proliferation [73] and immune system activity [74–76]. In general, TRIM22, IFI6, IFITM1, IFI35, and IRF9 are IFN-inducible genes that play significant functions in antiviral, cell proliferation, differentiation, apoptosis, and immune regulation functions. Previous studies have generally studied IFN-inducible genes in DM muscle tissue, but fewer have been conducted on PM tissues [13, 14, 17]. Recently, RNA sequence analysis of muscle biopsies in patients with DM and other forms of myositis (not including PM) suggested IFN1- and IFN2-inducible genes are differentially activated in various kinds of myositis, especially in DM [8]. Moreover, the role of IFN-inducible genes in PM and DM blood samples is controversial [15, 77, 78]. Therefore, bioinformatics analyses are needed to aggregate the findings of previous studies. This investigation showed that IFN-inducible genes are related to PM and DM by analyzing the data of previously performed transcriptome microarrays.

Another important finding of our study is the role of immune cell infiltration in PM/DM. While immune infiltration has been well studied in other autoimmune diseases, it had not yet been studied in PM/DM. We utilized CIBERSORT to perform a widespread assessment of the immune microenvironment in PM/DM, discovering elevated levels of macrophage infiltration, while resting Tregs infiltrated less; these findings may be correlated to the progression of PM/DM. Previous studies on macrophage differentiation patterns in PM/DM have demonstrated that DM and PM have different activation patterns and macrophage distribution. Specifically, in PM patients, macrophages activated in the early stage were mainly located in the myointima, which was the dominant histological manifestation. Whereas in DM, macrophages activated in late stages mainly appeared in the perifascicular area [79]. Other studies on macrophages and PM/DM have mainly focused on the activated receptors on the macrophage surface, such as CD206 and CD163. These studies have shown that activated macrophages are strongly correlated with the poor prognosis and mortality of PM/DM [80–82]. Another concern is that macrophage activation syndrome is caused by phagocytosis of hematopoietic components by activated macrophages. Macrophage activation syndrome is more common in RA and SLE; however, it has also been reported in some cases of DM [83, 84]. The above reports combined with our results have illustrated that macrophage infiltration is vital in the occurrence and progression of PM/DM.

The purpose of this investigation was to appraise biomarkers of PM/DM and to further investigate the function of immune cell infiltration in PM/DM. Our research has some limitations. First, we did not conduct in vivo tests to verify these outcomes. Second, we need to further study the definite mechanism of immune response caused by the five hub genes. Third, the CIBERSORT method is based upon finite genetic data, which might diverge from the heterogeneous interrelations of cells, characteristics of diseases, or plasticity of phenotypes. Finally, Mariampillai et al. [85] used unsupervised multiple correspondence analysis and hierarchical clustering analysis to aggregate patients by a database of the French myositis network according to myositis-specific antibodies, transcriptomic signatures, and clinicopathological correlations. They reclassified IIMs into four subgroups: DM, IMNM, ASS, and IBM [85, 86] and put forward the view that patients with PM were mainly present in IMNM and ASS clusters, and the use of PM should probably be discontinued. According to the latest classification, the PM patients we selected may be mixed with IMNM and ASS patients. However, as the original data of the five GEO series were not serologically grouped and discussed, it is not clear whether PM patients had a mixture of IMMN and ASS. Moreover, we did not explore the association of the five hub genes with the serological phenotypes (autoantibody profiles) of PM/DM patients. Although bioinformatics can reveal the internal mechanism, the results of our study need to be further validated by in vivo and in vitro tests and medical analysis.

We have comprehensively supplied a profound understanding of the overall molecular changes in the pathological mechanism of PM/DM and recognized five hub genes as potential therapeutic targets, including TRIM22, IFI6, IFITM1, IFI35, and IRF9. Moreover, through functional enrichment and pathway analysis, we discovered that these DEGs were generally enriched in immune function and the IFN signaling pathway. The WGCNA method identified MEblue and MEturquoise as key target modules in PM/DM. To our knowledge, this is the first study to focus on the role of immune cell infiltration in PM/DM. Additionally, we found that macrophage infiltration is influential in the occurrence and progression of PM/DM.

Author Contributions

SC conceived and designed the research. HL and HZ extracted data and conducted quality assessment. SC, XZ, HY, and YL analyzed the data. SC wrote the paper. All authors are accountable for all aspects of the study, and attest to the accuracy and integrity of the results. All authors contributed to the article and approved the submitted version.

Conflicts of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Funding

This work was supported by the National Natural Science Foundation of China Grants (No. 81770466, 81671618, 81871302, and 81800435); the Youth Plan of Beijing Hospital Management Center (QML20190602); the Beijing Outstanding Young Talents backbone individual project (2018000021469G242); the National Key Research and Development Program of China (2018YFE0207300); the CAMS Innovation Fund for Medical Sciences (CIFMS 2017-I2M-3-001and CIFMS 017-I2M-B&R-01); and the Beijing Key Clinical Specialty for Laboratory Medicine - Excellent Project (No. ZK201000).

References

  • 1. Meyer A, Meyer N, Schaeffer M, Gottenberg JE, Geny B, Sibilia J. Incidence and prevalence of inflammatory myopathies: a systematic review. Rheumatology (Oxford). 2015; 54:50–63. https://doi.org/10.1093/rheumatology/keu289 [PubMed]
  • 2. Rider LG, Miller FW. Deciphering the clinical presentations, pathogenesis, and treatment of the idiopathic inflammatory myopathies. JAMA. 2011; 305:183–90. https://doi.org/10.1001/jama.2010.1977 [PubMed]
  • 3. Behan WM, Behan PO. Immunological features of polymyositis/dermatomyositis. Springer Semin Immunopathol. 1985; 8:267–93. https://doi.org/10.1007/BF00197300 [PubMed]
  • 4. Yang WM, Chen JJ. Advances in biomarkers for dermatomyositis. Clin Chim Acta. 2018; 482:172–7. https://doi.org/10.1016/j.cca.2018.03.042 [PubMed]
  • 5. Chen Z, Cao M, Plana MN, Liang J, Cai H, Kuwana M, Sun L. Utility of anti-melanoma differentiation-associated gene 5 antibody measurement in identifying patients with dermatomyositis and a high risk for developing rapidly progressive interstitial lung disease: a review of the literature and a meta-analysis. Arthritis Care Res (Hoboken). 2013; 65:1316–24. https://doi.org/10.1002/acr.21985 [PubMed]
  • 6. Ceribelli A, De Santis M, Isailovic N, Gershwin ME, Selmi C. The Immune Response and the Pathogenesis of Idiopathic Inflammatory Myositis: a Critical Review. Clin Rev Allergy Immunol. 2017; 52:58–70. https://doi.org/10.1007/s12016-016-8527-x [PubMed]
  • 7. Moneta GM, Pires Marafon D, Marasco E, Rosina S, Verardo M, Fiorillo C, Minetti C, Bracci-Laudiero L, Ravelli A, De Benedetti F, Nicolai R. Muscle Expression of Type I and Type II Interferons Is Increased in Juvenile Dermatomyositis and Related to Clinical and Histologic Features. Arthritis Rheumatol. 2019; 71:1011–21. https://doi.org/10.1002/art.40800 [PubMed]
  • 8. Pinal-Fernandez I, Casal-Dominguez M, Derfoul A, Pak K, Plotz P, Miller FW, Milisenda JC, Grau-Junyent JM, Selva-O’Callaghan A, Paik J, Albayda J, Christopher-Stine L, Lloyd TE, et al. Identification of distinctive interferon gene signatures in different types of myositis. Neurology. 2019; 93:e1193–204. https://doi.org/10.1212/WNL.0000000000008128 [PubMed]
  • 9. Chinoy H, Li CK, Platt H, Fertig N, Varsani H, Gunawardena H, Betteridge Z, Oddis CV, McHugh NJ, Wedderburn LR, Ollier WE, Cooper RG, and UK Adult Onset Myositis Immunogenetic Consortium and UK Juvenile Dermatomyositis Research Group. Genetic association study of NF-κB genes in UK Caucasian adult and juvenile onset idiopathic inflammatory myopathy. Rheumatology (Oxford). 2012; 51:794–9. https://doi.org/10.1093/rheumatology/ker379 [PubMed]
  • 10. Chinoy H, Salway F, John S, Fertig N, Tait BD, Oddis CV, Ollier WE, Cooper RG, and UK Adult Onset Myositis Immunogenetic Collaboration (AOMIC). Tumour necrosis factor-alpha single nucleotide polymorphisms are not independent of HLA class I in UK Caucasians with adult onset idiopathic inflammatory myopathies. Rheumatology (Oxford). 2007; 46:1411–6. https://doi.org/10.1093/rheumatology/kem145 [PubMed]
  • 11. De Paepe B. Progressive Skeletal Muscle Atrophy in Muscular Dystrophies: A Role for Toll-like Receptor-Signaling in Disease Pathogenesis. Int J Mol Sci. 2020; 21:4440. https://doi.org/10.3390/ijms21124440 [PubMed]
  • 12. Zhang L, Xia Q, Li W, Peng Q, Yang H, Lu X, Wang G. The RIG-I pathway is involved in peripheral T cell lymphopenia in patients with dermatomyositis. Arthritis Res Ther. 2019; 21:131. https://doi.org/10.1186/s13075-019-1905-z [PubMed]
  • 13. Greenberg SA, Pinkus JL, Pinkus GS, Burleson T, Sanoudou D, Tawil R, Barohn RJ, Saperstein DS, Briemberg HR, Ericsson M, Park P, Amato AA. Interferon-alpha/beta-mediated innate immune mechanisms in dermatomyositis. Ann Neurol. 2005; 57:664–78. https://doi.org/10.1002/ana.20464 [PubMed]
  • 14. Wong D, Kea B, Pesich R, Higgs BW, Zhu W, Brown P, Yao Y, Fiorentino D. Interferon and biologic signatures in dermatomyositis skin: specificity and heterogeneity across diseases. PLoS One. 2012; 7:e29161. https://doi.org/10.1371/journal.pone.0029161 [PubMed]
  • 15. Walsh RJ, Kong SW, Yao Y, Jallal B, Kiener PA, Pinkus JL, Beggs AH, Amato AA, Greenberg SA. Type I interferon-inducible gene expression in blood is present and reflects disease activity in dermatomyositis and polymyositis. Arthritis Rheum. 2007; 56:3784–92. https://doi.org/10.1002/art.22928 [PubMed]
  • 16. Baechler EC, Bauer JW, Slattery CA, Ortmann WA, Espe KJ, Novitzke J, Ytterberg SR, Gregersen PK, Behrens TW, Reed AM. An interferon signature in the peripheral blood of dermatomyositis patients is associated with disease activity. Mol Med. 2007; 13:59–68. https://doi.org/10.2119/2006-00085.Baechler [PubMed]
  • 17. Salajegheh M, Kong SW, Pinkus JL, Walsh RJ, Liao A, Nazareno R, Amato AA, Krastins B, Morehouse C, Higgs BW, Jallal B, Yao Y, Sarracino DA, et al. Interferon-stimulated gene 15 (ISG15) conjugates proteins in dermatomyositis muscle with perifascicular atrophy. Ann Neurol. 2010; 67:53–63. https://doi.org/10.1002/ana.21805 [PubMed]
  • 18. Kim S, Rhee JK, Yoo HJ, Lee HJ, Lee EJ, Lee JW, Yu JH, Son BH, Gong G, Kim SB, Singh SR, Ahn SH, Chang S. Bioinformatic and metabolomic analysis reveals miR-155 regulates thiamine level in breast cancer. Cancer Lett. 2015; 357:488–97. https://doi.org/10.1016/j.canlet.2014.11.058 [PubMed]
  • 19. Liu J, Li Y, Gan Y, Xiao Q, Tian R, Shu G, Yin G. Identification of ZNF26 as a Prognostic Biomarker in Colorectal Cancer by an Integrated Bioinformatic Analysis. Front Cell Dev Biol. 2021; 9:671211. https://doi.org/10.3389/fcell.2021.671211 [PubMed]
  • 20. Topno R, Singh I, Kumar M, Agarwal P. Integrated bioinformatic analysis identifies UBE2Q1 as a potential prognostic marker for high grade serous ovarian cancer. BMC Cancer. 2021; 21:220. https://doi.org/10.1186/s12885-021-07928-z [PubMed]
  • 21. Barrett T, Wilhite SE, Ledoux P, Evangelista C, Kim IF, Tomashevsky M, Marshall KA, Phillippy KH, Sherman PM, Holko M, Yefanov A, Lee H, Zhang N, et al. NCBI GEO: archive for functional genomics data sets--update. Nucleic Acids Res. 2013; 41:D991–5. https://doi.org/10.1093/nar/gks1193 [PubMed]
  • 22. Greenberg SA, Bradshaw EM, Pinkus JL, Pinkus GS, Burleson T, Due B, Bregoli L, O’Connor KC, Amato AA. Plasma cells in muscle in inclusion body myositis and polymyositis. Neurology. 2005; 65:1782–7. https://doi.org/10.1212/01.wnl.0000187124.92826.20 [PubMed]
  • 23. Zhu W, Streicher K, Shen N, Higgs BW, Morehouse C, Greenlees L, Amato AA, Ranade K, Richman L, Fiorentino D, Jallal B, Greenberg SA, Yao Y. Genomic signatures characterize leukocyte infiltration in myositis muscles. BMC Med Genomics. 2012; 5:53. https://doi.org/10.1186/1755-8794-5-53 [PubMed]
  • 24. Greenberg SA, Pinkus JL, Kong SW, Baecher-Allan C, Amato AA, Dorfman DM. Highly differentiated cytotoxic T cells in inclusion body myositis. Brain. 2019; 142:2590–604. https://doi.org/10.1093/brain/awz207 [PubMed]
  • 25. Gautier L, Cope L, Bolstad BM, Irizarry RA. affy--analysis of Affymetrix GeneChip data at the probe level. Bioinformatics. 2004; 20:307–15. https://doi.org/10.1093/bioinformatics/btg405 [PubMed]
  • 26. Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. 2012; 28:882–3. https://doi.org/10.1093/bioinformatics/bts034 [PubMed]
  • 27. Liao Y, Smyth GK, Shi W. The R package Rsubread is easier, faster, cheaper and better for alignment and quantification of RNA sequencing reads. Nucleic Acids Res. 2019; 47:e47. https://doi.org/10.1093/nar/gkz114 [PubMed]
  • 28. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015; 43:e47. https://doi.org/10.1093/nar/gkv007 [PubMed]
  • 29. Wickham H. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. 2016. https://doi.org/10.1007/978-3-319-24277-4
  • 30. Kolde R, Laur S, Adler P, Vilo J. Robust rank aggregation for gene list integration and meta-analysis. Bioinformatics. 2012; 28:573–80. https://doi.org/10.1093/bioinformatics/btr709 [PubMed]
  • 31. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–7. https://doi.org/10.1089/omi.2011.0118 [PubMed]
  • 32. Walter W, Sánchez-Cabo F, Ricote M. GOplot: an R package for visually combining expression data with functional analysis. Bioinformatics. 2015; 31:2912–4. https://doi.org/10.1093/bioinformatics/btv300 [PubMed]
  • 33. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 2015; 1:417–25. https://doi.org/10.1016/j.cels.2015.12.004 [PubMed]
  • 34. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013; 14:7. https://doi.org/10.1186/1471-2105-14-7 [PubMed]
  • 35. Kolde R, Kolde MR. Package ‘pheatmap’. R package, 2015; 1:790. https://CRAN.R-project.org/package=pheatmap.
  • 36. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015; 12:453–7. https://doi.org/10.1038/nmeth.3337 [PubMed]
  • 37. Szklarczyk D, Gable AL, Nastou KC, Lyon D, Kirsch R, Pyysalo S, Doncheva NT, Legeay M, Fang T, Bork P, Jensen LJ, von Mering C. The STRING database in 2021: customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res. 2021; 49:D605–12. https://doi.org/10.1093/nar/gkaa1074 [PubMed]
  • 38. Bauer-Mehren A. Integration of genomic information with biological networks using Cytoscape. Methods Mol Biol. 2013; 1021:37–61. https://doi.org/10.1007/978-1-62703-450-0_3 [PubMed]
  • 39. Chen H, Boutros PC. VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinformatics. 2011; 12:35. https://doi.org/10.1186/1471-2105-12-35 [PubMed]
  • 40. Taiyun W, Viliam S. R package “corrplot”: Visualization of a Correlation Matrix. R package Version 0.84. 2017. https://github.com/taiyun/corrplot.
  • 41. Alboukadel K. ggcorrplot: Visualization of a Correlation Matrix using 'ggplot2'. R package version 0.1.3. 2019. https://CRAN.R-project.org/package=ggcorrplot.
  • 42. Suárez-Calvet X, Gallardo E, Nogales-Gadea G, Querol L, Navas M, Díaz-Manera J, Rojas-Garcia R, Illa I. Altered RIG-I/DDX58-mediated innate immunity in dermatomyositis. J Pathol. 2014; 233:258–68. https://doi.org/10.1002/path.4346 [PubMed]
  • 43. Robin X, Turck N, Hainard A, Tiberti N, Lisacek F, Sanchez JC, Müller M. pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics. 2011; 12:77. https://doi.org/10.1186/1471-2105-12-77 [PubMed]
  • 44. El Khouli RH, Macura KJ, Barker PB, Habba MR, Jacobs MA, Bluemke DA. Relationship of temporal resolution to diagnostic performance for dynamic contrast enhanced MRI of the breast. J Magn Reson Imaging. 2009; 30:999–1004. https://doi.org/10.1002/jmri.21947 [PubMed]
  • 45. Dalakas MC, Hohlfeld R. Polymyositis and dermatomyositis. Lancet. 2003; 362:971–82. https://doi.org/10.1016/S0140-6736(03)14368-1 [PubMed]
  • 46. Zhou X, Dimachkie MM, Xiong M, Tan FK, Arnett FC. cDNA microarrays reveal distinct gene expression clusters in idiopathic inflammatory myopathies. Med Sci Monit. 2004; 10:BR191–7. [PubMed]
  • 47. Lundberg I, Brengman JM, Engel AG. Analysis of cytokine expression in muscle in inflammatory myopathies, Duchenne dystrophy, and non-weak controls. J Neuroimmunol. 1995; 63:9–16. https://doi.org/10.1016/0165-5728(95)00122-0 [PubMed]
  • 48. Illa I, Gallardo E, Gimeno R, Serrano C, Ferrer I, Juárez C. Signal transducer and activator of transcription 1 in human muscle: implications in inflammatory myopathies. Am J Pathol. 1997; 151:81–8. [PubMed]
  • 49. Ivashkiv LB, Donlin LT. Regulation of type I interferon responses. Nat Rev Immunol. 2014; 14:36–49. https://doi.org/10.1038/nri3581 [PubMed]
  • 50. Pestka S, Krause CD, Walter MR. Interferons, interferon-like cytokines, and their receptors. Immunol Rev. 2004; 202:8–32. https://doi.org/10.1111/j.0105-2896.2004.00204.x [PubMed]
  • 51. Barr SD, Smiley JR, Bushman FD. The interferon response inhibits HIV particle production by induction of TRIM22. PLoS Pathog. 2008; 4:e1000007. https://doi.org/10.1371/journal.ppat.1000007 [PubMed]
  • 52. Jefferies C, Wynne C, Higgs R. Antiviral TRIMs: friend or foe in autoimmune and autoinflammatory disease? Nat Rev Immunol. 2011; 11:617–25. https://doi.org/10.1038/nri3043 [PubMed]
  • 53. Li Q, Lee CH, Peters LA, Mastropaolo LA, Thoeni C, Elkadri A, Schwerd T, Zhu J, Zhang B, Zhao Y, Hao K, Dinarzo A, Hoffman G, et al. Variants in TRIM22 That Affect NOD2 Signaling Are Associated With Very-Early-Onset Inflammatory Bowel Disease. Gastroenterology. 2016; 150:1196–207. https://doi.org/10.1053/j.gastro.2016.01.031 [PubMed]
  • 54. Ashton JJ, Mossotto E, Stafford IS, Haggarty R, Coelho TAF, Batra A, Afzal NA, Mort M, Bunyan D, Beattie RM, Ennis S. Genetic Sequencing of Pediatric Patients Identifies Mutations in Monogenic Inflammatory Bowel Disease Genes that Translate to Distinct Clinical Phenotypes. Clin Transl Gastroenterol. 2020; 11:e00129. https://doi.org/10.14309/ctg.0000000000000129 [PubMed]
  • 55. Raposo RA, Gupta R, Abdel-Mohsen M, Dimon M, Debbaneh M, Jiang W, York VA, Leadabrand KS, Brown G, Malakouti M, Arron S, Kuebler PJ, Wu JJ, et al. Antiviral gene expression in psoriasis. J Eur Acad Dermatol Venereol. 2015; 29:1951–7. https://doi.org/10.1111/jdv.13091 [PubMed]
  • 56. Morris G, Maes M, Murdjeva M, Puri BK. Do Human Endogenous Retroviruses Contribute to Multiple Sclerosis, and if So, How? Mol Neurobiol. 2019; 56:2590–605. https://doi.org/10.1007/s12035-018-1255-x [PubMed]
  • 57. Martire S, Navone ND, Montarolo F, Perga S, Bertolotto A. A gene expression study denies the ability of 25 candidate biomarkers to predict the interferon-beta treatment response in multiple sclerosis patients. J Neuroimmunol. 2016; 292:34–9. https://doi.org/10.1016/j.jneuroim.2016.01.010 [PubMed]
  • 58. Nexø BA, Hansen B, Nissen KK, Gundestrup L, Terkelsen T, Villesen P, Bahrami S, Petersen T, Pedersen FS, Laska MJ. Restriction genes for retroviruses influence the risk of multiple sclerosis. PLoS One. 2013; 8:e74063. https://doi.org/10.1371/journal.pone.0074063 [PubMed]
  • 59. Cheriyath V, Kaur J, Davenport A, Khalel A, Chowdhury N, Gaddipati L. G1P3 (IFI6), a mitochondrial localised antiapoptotic protein, promotes metastatic potential of breast cancer cells through mtROS. Br J Cancer. 2018; 119:52–64. https://doi.org/10.1038/s41416-018-0137-3 [PubMed]
  • 60. Meyer K, Kwon YC, Liu S, Hagedorn CH, Ray RB, Ray R. Interferon-α inducible protein 6 impairs EGFR activation by CD81 and inhibits hepatitis C virus infection. Sci Rep. 2015; 5:9012. https://doi.org/10.1038/srep09012 [PubMed]
  • 61. Sol N, Leurs CE, Veld SGI', Strijbis EM, Vancura A, Schweiger MW, Teunissen CE, Mateen FJ, Tannous BA, Best MG, Würdinger T, Killestein J. Blood platelet RNA enables the detection of multiple sclerosis. Mult Scler J Exp Transl Clin. 2020; 6:2055217320946784. https://doi.org/10.1177/2055217320946784 [PubMed]
  • 62. Rodríguez-Carrio J, López P, Alperi-López M, Caminal-Montero L, Ballina-García FJ, Suárez A. IRF4 and IRGs Delineate Clinically Relevant Gene Expression Signatures in Systemic Lupus Erythematosus and Rheumatoid Arthritis. Front Immunol. 2019; 9:3085. https://doi.org/10.3389/fimmu.2018.03085 [PubMed]
  • 63. Ugolini-Lopes MR, Torrezan GT, Gândara APR, Olivieri EHR, Nascimento IS, Okazaki E, Bonfá E, Carraro DM, de Andrade DCO. Enhanced type I interferon gene signature in primary antiphospholipid syndrome: Association with earlier disease onset and preeclampsia. Autoimmun Rev. 2019; 18:393–8. https://doi.org/10.1016/j.autrev.2018.11.004 [PubMed]
  • 64. Zhao X, Zhang L, Wang J, Zhang M, Song Z, Ni B, You Y. Identification of key biomarkers and immune infiltration in systemic lupus erythematosus by integrated bioinformatics analysis. J Transl Med. 2021; 19:35. https://doi.org/10.1186/s12967-020-02698-x [PubMed]
  • 65. Zhang YJ, Sun YZ, Gao XH, Qi RQ. Integrated bioinformatic analysis of differentially expressed genes and signaling pathways in plaque psoriasis. Mol Med Rep. 2019; 20:225–35. https://doi.org/10.3892/mmr.2019.10241 [PubMed]
  • 66. Wang J, Wang CF, Ming SL, Li GL, Zeng L, Wang MD, Su BQ, Wang Q, Yang GY, Chu BB. Porcine IFITM1 is a host restriction factor that inhibits pseudorabies virus infection. Int J Biol Macromol. 2020; 151:1181–93. https://doi.org/10.1016/j.ijbiomac.2019.10.162 [PubMed]
  • 67. Liang R, Li X, Zhu X. Deciphering the Roles of IFITM1 in Tumors. Mol Diagn Ther. 2020; 24:433–41. https://doi.org/10.1007/s40291-020-00469-4 [PubMed]
  • 68. Das A, Dinh PX, Pattnaik AK. Trim21 regulates Nmi-IFI35 complex-mediated inhibition of innate antiviral response. Virology. 2015; 485:383–92. https://doi.org/10.1016/j.virol.2015.08.013 [PubMed]
  • 69. Xiahou Z, Wang X, Shen J, Zhu X, Xu F, Hu R, Guo D, Li H, Tian Y, Liu Y, Liang H. NMI and IFP35 serve as proinflammatory DAMPs during cellular infection and injury. Nat Commun. 2017; 8:950. https://doi.org/10.1038/s41467-017-00930-9 [PubMed]
  • 70. Jian D, Wang W, Zhou X, Jia Z, Wang J, Yang M, Zhao W, Jiang Z, Hu X, Zhu J. Interferon-induced protein 35 inhibits endothelial cell proliferation, migration and re-endothelialization of injured arteries by inhibiting the nuclear factor-kappa B pathway. Acta Physiol (Oxf). 2018; 223:e13037. https://doi.org/10.1111/apha.13037 [PubMed]
  • 71. De Masi R, Orlando S. IFI35 as a biomolecular marker of neuroinflammation and treatment response in multiple sclerosis. Life Sci. 2020; 259:118233. https://doi.org/10.1016/j.lfs.2020.118233 [PubMed]
  • 72. Zhang L, Zhu H, Li Y, Dai X, Zhou B, Li Q, Zuo X, Luo H. The role of IFI35 in lupus nephritis and related mechanisms. Mod Rheumatol. 2017; 27:1010–8. https://doi.org/10.1080/14397595.2016.1270387 [PubMed]
  • 73. Weihua X, Lindner DJ, Kalvakolanu DV. The interferon-inducible murine p48 (ISGF3gamma) gene is regulated by protooncogene c-myc. Proc Natl Acad Sci USA. 1997; 94:7227–32. https://doi.org/10.1073/pnas.94.14.7227 [PubMed]
  • 74. Rauch I, Rosebrock F, Hainzl E, Heider S, Majoros A, Wienerroither S, Strobl B, Stockinger S, Kenner L, Müller M, Decker T. Noncanonical Effects of IRF9 in Intestinal Inflammation: More than Type I and Type III Interferons. Mol Cell Biol. 2015; 35:2332–43. https://doi.org/10.1128/MCB.01498-14 [PubMed]
  • 75. Smith S, Fernando T, Wu PW, Seo J, Ní Gabhann J, Piskareva O, McCarthy E, Howard D, O’Connell P, Conway R, Gallagher P, Molloy E, Stallings RL, et al. MicroRNA-302d targets IRF9 to regulate the IFN-induced gene expression in SLE. J Autoimmun. 2017; 79:105–11. https://doi.org/10.1016/j.jaut.2017.03.003 [PubMed]
  • 76. Huber M, Suprunenko T, Ashhurst T, Marbach F, Raifer H, Wolff S, Strecker T, Viengkhou B, Jung SR, Obermann HL, Bauer S, Xu HC, Lang PA, et al. IRF9 Prevents CD8+ T Cell Exhaustion in an Extrinsic Manner during Acute Lymphocytic Choriomeningitis Virus Infection. J Virol. 2017; 91:e01219–7. https://doi.org/10.1128/JVI.01219-17 [PubMed]
  • 77. Król P, Kryštůfková O, Polanská M, Mann H, Klein M, Beran O, Vencovský J. Serum levels of interferon α do not correlate with disease activity in patients with dermatomyositis/polymyositis. Ann Rheum Dis. 2011; 70:879–80. https://doi.org/10.1136/ard.2010.141051 [PubMed]
  • 78. Greenberg SA, Higgs BW, Morehouse C, Walsh RJ, Kong SW, Brohawn P, Zhu W, Amato A, Salajegheh M, White B, Kiener PA, Jallal B, Yao Y. Relationship between disease activity and type 1 interferon- and other cytokine-inducible gene expression in blood in dermatomyositis and polymyositis. Genes Immun. 2012; 13:207–13. https://doi.org/10.1038/gene.2011.61 [PubMed]
  • 79. Rostasy KM, Piepkorn M, Goebel HH, Menck S, Hanefeld F, Schulz-Schaeffer WJ. Monocyte/macrophage differentiation in dermatomyositis and polymyositis. Muscle Nerve. 2004; 30:225–30. https://doi.org/10.1002/mus.20088 [PubMed]
  • 80. Horiike Y, Suzuki Y, Fujisawa T, Yasui H, Karayama M, Hozumi H, Furuhashi K, Enomoto N, Nakamura Y, Inui N, Ogawa N, Suda T. Successful classification of macrophage-mannose receptor CD206 in severity of anti-MDA5 antibody positive dermatomyositis associated ILD. Rheumatology (Oxford). 2019; 58:2143–52. https://doi.org/10.1093/rheumatology/kez185 [PubMed]
  • 81. Zuo Y, Ye L, Liu M, Li S, Liu W, Chen F, Lu X, Gordon P, Wang G, Shu X. Clinical significance of radiological patterns of HRCT and their association with macrophage activation in dermatomyositis. Rheumatology (Oxford). 2020; 59:2829–37. https://doi.org/10.1093/rheumatology/keaa034 [PubMed]
  • 82. Peng QL, Zhang YL, Shu XM, Yang HB, Zhang L, Chen F, Lu X, Wang GC. Elevated Serum Levels of Soluble CD163 in Polymyositis and Dermatomyositis: Associated with Macrophage Infiltration in Muscle Tissue. J Rheumatol. 2015; 42:979–87. https://doi.org/10.3899/jrheum.141307 [PubMed]
  • 83. Paul N, Avalos C, Estifan E, Swyden S. Interstitial lung disease in dermatomyositis complicated by right ventricular thrombus secondary to macrophage activation syndrome- a case report. AME Case Rep. 2020; 4:18. https://doi.org/10.21037/acr.2020.03.06 [PubMed]
  • 84. Zhu DX, Qiao JJ, Fang H. Macrophage activation syndrome as a complication of dermatomyositis: A case report. World J Clin Cases. 2020; 8:2339–44. https://doi.org/10.12998/wjcc.v8.i11.2339 [PubMed]
  • 85. Mariampillai K, Granger B, Amelin D, Guiguet M, Hachulla E, Maurier F, Meyer A, Tohmé A, Charuel JL, Musset L, Allenbach Y, Benveniste O. Development of a New Classification System for Idiopathic Inflammatory Myopathies Based on Clinical Manifestations and Myositis-Specific Autoantibodies. JAMA Neurol. 2018; 75:1528–37. https://doi.org/10.1001/jamaneurol.2018.2598 [PubMed]
  • 86. Tanboon J, Uruha A, Stenzel W, Nishino I. Where are we moving in the classification of idiopathic inflammatory myopathies? Curr Opin Neurol. 2020; 33:590–603. https://doi.org/10.1097/WCO.0000000000000855 [PubMed]
Download PDF
Corresponding Author

Yongzhe Li
yongzhelipumch@126.com
https://orcid.org/0000-0002-8267-0985

Keywords
idiopathic inflammatory myopathy macrophage infiltration interferon signaling autoimmune disease gene set expression analysis
Paper Sections