GSEA enrichment analysis on DEGs of blood samples between healthy controls and RA patients
R package limma(21) was used to analyze the DEGs between RA and healthy controls. In order to include more DEGs, FDR <0.05 was selected as the threshold and the volcano map as shown in Fig.1A. Moreover, the fold differences of gene expression were used as the rank ordering, and GSEA was used for enrichment analysis of DEGs. As shown in Fig1B and C, compared with the healthy control sample, the up-regulated DEGs were mainly concentrated in the region with large difference folds, while the down-regulated DEGs were mainly concentrated in the region with small difference folds.
Identification of infliximab therapeutic response-related gene co-expression modules for RA
To further explore the biomarkers associated with the infliximab response, the expression profile data of GSE78068 was obtained (containing 140 samples). Afterwards, 7692 DEGsbetween samples with unresponsiveness and responsiveness against infliximab treatment were extracted and R software package WGCNA was subsequently used to screen the co-expression module. As a result, the co-expression network was shown to conform to the scale-free network, that was to say, the logarithm log(k) of the node with the connection degree k was negatively correlated with the logarithm log(P(k)) of the probability of the node, and the correlation coefficient should be greater than 0.8. To ensure that the network was a scale-free network, β=12 was set (Fig. 2A and B). A total of 19 modules were obtained (Fig. 2C), where the grey module was a collection of genes that cannot be aggregated into other modules. The transcript statistics of each module was shown in Table 1, showing that 4198 transcripts were assigned to 18 co-expression modules. We subsequently calculated the correlation between the eigenvector of the 19 modules and the therapeutic response of infliximab. As shown in Fig. 2D, the tan module exhibited the most significantly negative correlation with the therapeutic response of infliximab, suggesting that the genes in the tan module might be closely related to the therapeutic response of infliximab.
Table1:Statistics on the number of genes corresponding to each module
Module
|
Count
|
black
|
219
|
blue
|
444
|
brown
|
442
|
cyan
|
69
|
green
|
380
|
greenyellow
|
120
|
grey
|
3494
|
grey60
|
47
|
lightcyan
|
62
|
lightgreen
|
43
|
magenta
|
181
|
midnightblue
|
67
|
pink
|
202
|
purple
|
130
|
red
|
317
|
salmon
|
75
|
tan
|
102
|
turquoise
|
869
|
yellow
|
429
|
To further investigate the correlation between the genes in the tan module and the therapeutic response of infliximab, we analyzed the difference of the eigenvector of the tan module between infliximab-responsive samples and infliximab-unresponsive samples. As shown in Fig.3A, the eigenvector of the tan module was significantly higher in the infliximab-unresponsive samples than that in the infliximab-responsive samples, suggesting the expression levels of genes in the tan module may be more active in infliximab- unresponsive samples. We subsequently analyzed the difference of the levels of 102 genes of the tan module between infliximab-responsive samples and infliximab-unresponsive samples. As shown in Fig.3B, the gene expression levels were significantly higher in the infliximab-unresponsive samples than that in the infliximab-responsive samples, consistent with our expectation. Finally, we analyzed the differential expression of these genes in RA patients and healthy control samples. As shown in Fig.3C, the expression levels of these genes in the blood samples of RA patients were significantly higher than those of healthy controls, suggesting that the 102 genes in the tan module were not only related to the therapeutic response of infliximab in RA patients, but also critically involved in the occurrence and development of RA.
Functional enrichment analysis of the 102 genes in tan module
To investigate the gene functions in the tan module, the genes in the tan module were mapped to the g:profiler online analysis platform for enrichment analysis, with threshold FDR<0.05. As a result, these genes could be enriched into multiple GO terms, including 91 biological processes, mainly related to immune responses and stimuli defense (as shown in Fig.4A); five cellular components, mainly associated with cytoplasma and inflammasome (as shown in Fig.4B); five molecular functions, mainly related to protein and RNA binding (as shown in Fig.4C). Subsequently, these 102 genes were mainly enriched into seven KEGG pathways, mainly associated with inflammatory infection, NOD-like receptor signaling pathway (as shown in Fig.4D); and 13 reactom pathways, mainly correlated with immunity and the interferon signaling pathway (as shown in Fig.4E). Finally, the miRNAs enrichment analysis of these genes showed that 16 out of the 102 genes could be enriched into hsa-miR-146a-5p (as shown in Fig.4F). hsa-miR-146a-5p has been validated to be closely related to the occurrence and development of various complex diseases(22, 23). The serum level of hsa-miR-146a-5p has been reported as a diagnostic marker for RA and predictive marker for clinical response to TNF inhibitor(24, 25). As we all known, RA is a chronic systemic autoimmune disease characterized by joint lesions. In our study, the genes in the gene co-expression module tan, which was most significantly correlated with the therapeutic response of infliximab, are involved in the multiple immunological processes and relevant signaling transduction.
Mining and validation of key genes associated with the occurrence, development and infliximab response of RA diseases
In order to identify key gene markers, genes in the tan module were mapped to the String database, and combine score>0.4 was selected, which finally gave rise to 330 interactions of 63 genes (as shown in last.ppi.cm0.4.txt), followed by Cytoscope visualization (as shown in Fig.5A). The statistics of the degree of genes in the network were shown in Fig.5B, revealing relatively high degree of most genes, with an average of 10.47, which suggested the relatively close interaction among these genes. The degree distribution of the network was shown in Fig.5C, showing that the degree of most genes was relatively small, while the degree of a few genes was relatively large, which indicated that the close gene-gene interactions in the network may be dominated by a small number of genes with large degree. The Cytoscope plug-in cytohubba was subsequently used to analyze the topological properties of each node, and genes simultaneously within the top 20 of the degree, top 20 of the Closeness, and top 20 of the Betweenness were selected (as shown in lst.gene.txt) as key genes correlated with the occurrence, development and infliximab response in RA.
Afterwards, to confirm the relationship between these 14 genes and the therapeutic response of infliximab, the differential expression of these 14 genes between the infliximab-responsive samples and the infliximab-unresponsive samples of RA patients was analyzed. As shown in Fig.6, there were significant difference of 10 genes between the two groups, and the expression levels of the 10 genes (OAS1, STAT1, IFIT3, IFIT2, IFI6, RSAD2, OASL, GBP2, IRF7 and IFIT1)were significantly lower in the infliximab-responsive samples than those in the infliximab-unresponsive samples. This finding indicated that these genes were associated with the therapeutic response of infliximab, which could be used as a blood biomarker for individualized treatment of infliximab (even anti-TNF therapy) in RA patients.
Finally, in order to investigate the relationship of these 14 genes with the occurrence and development of RA, the expression difference of the 14 genes were analyzed between RA samples and healthy control samples. As shown in Fig. 7, there were significant difference of the 14 genes between the two groups, and the expression levels of these genes were significantly higher in RA samples than those in healthy controls, suggesting that these 14 genes were associated with the occurrence and development of RA, which could be used as a blood diagnostic marker or new therapeutic targets for RA.
As a summary (Fig.8), these 14 genes play an important role in the occurrence and development of RA, and further high-expression of 10 out of them will lead to a decreased responsiveness of infliximab in RA patients.