Comparative host transcriptomics as a tool to identify candidate biomarkers for immune reactions in leprosy using meta-analysis

Introduction

Leprosy reactions are acute inflammatory episodes that occur during infection with Mycobacterium leprae.¹ These reactions pose a high risk to the infected individuals, resulting in physical deformity followed by nerve damage. Leprosy reactions are classified into two types: Type-I or reversal reactions and Type-II or commonly presenting as erythema nodosum leprosum (ENL).² Both types of reactions may occur in leprosy patients before and during treatment, but may also occur after the treatment is completed.³ Type-I reactions occur due to an increase in cell-mediated immunity which results in nerve inflammation at the point of mycobacterial infection. The skin lesions become tender and nerve involvement produces pain, abnormal sensation and functional impairment.⁴ Prompt intervention is required to prevent permanent nerve impairment since the expected rate of recovery for nerve function is 60–70%.⁵ Type-II reactions occur in multibacillary patients and are immune complex mediated. Patients with Type-II reactions develop tender, erythematous nodules on the face, trunk or legs.⁶ The treatment regimen for ENL is quite long and delayed treatment is an important cause of neuropathy and consequent disability.⁷ Thus, timely management of leprosy reactions is a big challenge and there is a need to understand the underlying complex immunological processes and gene expression.

Transcriptome profiling can be viewed as a vital approach in disseminating valuable information related to gene expression and associated host response to infection, assessing disease severity and discovering novel gene candidates for diagnosis and prognosis. Due to differences in experimental design and multiple hypotheses in individual studies, these differentially expressed genes (DEGs) should be interpreted accordingly.⁸ However, their combined analysis can unravel some common pathological mechanisms. In this study, the meta-analysis of individual datasets paved the way to integrate gene sets from multiple individual studies and work out a common list of DEGs for further biological investigation.⁹ In this study, leprosy reaction–specific datasets were analysed and integrated using meta-analysis to identify DEGs that could be considered as candidate biomarkers in leprosy pathogenesis. Further, the list of DEGs was subject to pathway analysis revealing vital processes and their likely association.

Type-I or reversal reaction represents an elevation in cellular immune response which progresses towards the tuberculoid pole of the disease. This type of reaction is common in borderline forms of leprosy and does not occur in the polar forms. The clinical expression of Type-I or reversal reaction includes an increase in the levels of activated lymphocytes which causes the existing skin lesions to become erythematous and oedematous. If it is very severe, then ulceration may occur at times. Peripheral nerves become enlarged and tender due to acute neuritis. It results in rapid loss of sensory and motor functions. The duration of the reversal reaction is often several weeks, but it may extend for many months.

Type-II reaction can be chronic and may occur in half of the patients with lepromatous or borderline lepromatous type of leprosy during the disease. The entire immunology of ENL is not yet completely understood but the condition is aggravated by immune complexes formed with Mycobacterium leprae antigens. The clinical expression includes neuritis, iridocyclitis, and painful skin lesions which are erythematous papules and nodules. Chronic ENL is associated with amyloidosis and glomerulonephritis and may cause mortality from renal insufficiency.

Methods

In the present study, we selected four high throughput sequencing and microarray datasets from a whole blood sample and skin biopsy of patients with adverse leprosy reactions from the GEO (Gene Expression Omnibus) database¹⁰ and investigated clinical biomarker candidates for early diagnosis and prognosis of immune reactions of leprosy based on functional and molecular pathway analyses of DEGs. We selected four datasets of whole blood samples from leprosy patients with immune reactions, GSE16844, GSE125943, GSE129033, and GSE74481, using the following keywords: ‘reversal reaction’, ‘leprosy’, ‘Homo sapiens,’ and ‘erythema nodosum leprosum’. These datasets were downloaded from the GEO database. The search results were examined and datasets were excluded if they consisted of (i) only the pathogen gene expression, (ii) only a single sample that does not allow comparison, (iii) treatment by any compound, and (iv) studies targeting only non-coding RNAs. Raw data were preferentially downloaded for the selected studies wherever possible; otherwise, the gene expression matrix was used. The datasets were included when (i) they consist of only host gene expression, (ii) they contain samples from either whole blood or skin biopsy and (iii) they target the mRNA. No studies were discarded due to differences in the sequencing platform, publication or experiment date, or unavailability of raw data as it would have further lowered the number of datasets available. A summary of the individual datasets is shown in Table 1. The workflow of the analysis is depicted in Figure 1.

Close

Figure 1a:

Flowchart depicting the meta-analysis procedure. (SR - Sum of ranks, REM - Random effects model.)

Export to PPT

Close

Figure 1b:

Differentially expressed genes selected after the meta-analysis procedure for individual comparisons using multiple methods (sdef, maxP, SR, fisher), LL - Lepromatous leprosy, ENL - Erythema nodosum leprosum, SR - Sum of Ranks, RR - Reversal reaction, HC - Healthy control.

Export to PPT

Meta-analysis

As proposed by the developers of the sdef R package, it aims for fewer false positives when compared to other methods. Thus, a comparative approach was incorporated to compare genes with other standard meta-analysis methods with those obtained from the sdef method. MetaDE R package¹⁸ was used for differential gene expression. MetaDE package incorporates major meta-analysis methods such as Fisher, maximum P-value (maxP), the sum of ranks (SR), rth ordered P-value (roP), and random effects model (REM) which is based on effect sizes and are tested on the datasets used in this study. MetaDE also provides options for gene matching across studies and gene filtering before meta-analysis. In this study, P-value-based methods such as maxP, SR, and roP were used. Combining P-values from different studies has two major advantages: (i) simplicity and extensibility to different outcome variables and (ii) when the outcome variable is not binary then effect size may not be computed; however, P-value association among studies can still be calculated. In REM, the estimated effect size in each study is assumed to come from an underlying true effect size which also includes measurement errors like experimental or sampling error. Further, the studies also contain a random effect that can incorporate unknown cross-study heterogeneities in the model.¹⁹

Results

Identification of common differentially expressed genes using sdef and meta-analysis methods

Individual datasets were processed to identify all the expressed genes. Then, similar comparisons from different datasets were considered to identify DEGs. Four similar comparisons from datasets were considered for ratio association analysis, namely, HC vs ENL, HC vs RR, LL vs ENL, and LL vs RR. Genes common to each study within the specific comparison were used in the ratio association test (sdef). Additional gene expression meta-analysis tools were used as described in the MetaDE R package, such as- rth ordered P-value (roP), sum of ranks (SR), and maximum P-value (maxP), among others. Figure 2 shows the upset plots representing the number of genes selected by each method along with their intersection with a false discovery rate (FDR) < 0.01. maxP and Fisher seems to be a more generous method in the selection of DEGs across comparisons. sdef analysis showed a significant association between the expressed genes in the datasets within each comparison category. Overall, the sdef method seems to be moderate in selection of genes. Other meta-analysis methods are quite conservative as they select fewer genes than the ratio association method. Table 2 represents the list of genes selected with |log₂FC| > 3 and FDR < 0.001 in the HC vs ENL comparison. A total of 41 genes were selected by all four meta-analyses methods. maxP, fisher, and SR methods selected 29 genes. Top DEGs which are upregulated in HC individuals are ATF3, BCL2A1, C11orf96, C9orf47, CHST15, PLAUR and S100A9 while those more expressed in ENL are – AGFG2, ANK3, CD1E, CTSF, DOCK9, KRT73, NAV2, PCDH1, among others.

Close

Figure 2a:

Upset plots depicting the number of differentially expressed genes (DEGs) from each meta-analysis method and their intersection. The vertical bar represents the total number of expressed genes from different methods (intersection size) and the horizontal bar depicts the total gene size (set size) expressed by the meta-analysis method. DEGs from HC vs ENL comparison from two independent datasets (GSE129033 and GSE74481), 41 genes were identified by all the methods.

Export to PPT

Close

Figure 2b:

Upset plots depicting the number of differentially expressed genes from each meta-analysis method and their intersection. Th vertical bar represents the total number of expressed genes from different methods (intersection size) and the horizontal bar depicts the total gene size (set size) expressed by the meta-analysis method. Genes from HC vs RR from two datasets (GSE129033 and GSE74481), 13 genes were found common by meta-analysis methods.

Export to PPT

Close

Figure 2c:

Upset plots depicting the number of differentially expressed genes from each meta-analysis method and their intersection. The vertical bar represents the total number of expressed genes from different methods (intersection size) and the horizontal bar depicts the total gene size (set size) expressed by the meta-analysis method. Differentially expressed genes resulted from LL vs ENL comparison from two datasets (GSE16844 and GSE74481), 41 DEGs were found.

Export to PPT

Close

Figure 2d:

Upset plots depicting the number of differentially expressed genes (DEGs) from each meta-analysis method and their intersection. The vertical bar represents the total number of expressed genes from different methods (intersection size) and the horizontal bar depicts the total gene size (set size) expressed by the meta-analysis method. DEGs from LL vs RR comparison (GSE125943 and GSE74481), a total of 44 DEGs were selected by all meta-analysis methods.

Export to PPT

Table 3 represents the DEGs selected. A cut-off of |log₂FC| > 2 and FDR < 0.001 in LL vs ENL is considered. A total of 334 genes were selected by the random effects model (REM) method excluding SR; however, 41 genes were found to be significant by all four methods along with REM. Genes that showed more expression in LL are ATP1B1, CD9, ITGBL1, PRKACB and PSD3 while genes more expressed in ENL are ADAMTS5, ADAMTS9, IFITM2, IFITM3, KIRREL, KRT16, KRT6A, PTX3, SERPINA3, SLC22A4, among others.

DEGs selected by all four methods (Fisher, SR, maxP, and sdef) in HC vs RR and LL vs RR comparison groups are shown in Supplementary Table 1 and Supplementary Table 2, respectively.

Comparison among groups considered in this study

The number of DEGs from the comparison groups such as HC vs ENL, HC vs RR, LL vs ENL, and LL vs RR were then further compared to reveal genes specific to an individual group such as healthy, LL, ENL, and RR. The DEGs that belong to HC, LL, ENL, and RR were found to be 2536, 2861, 2716 and 2585, respectively. The Venn diagram depicting the genes specific to individual groups and common to different groups is shown in Figure 3.

Close

Figure 3:

Venn diagram depicting the number of differentially expressed genes specific to each group and common to these groups (healthy control, LL - Lepromatous, ENL - Erythema nodosum leprosum, RR - Reversal reaction).

Export to PPT

Gene ontology and functional enrichment analysis

The list of DEGs from different studies across multiple comparisons is used to understand the biological role of these genes. In the HC vs ENL comparison, 41 genes were selected by all four methods with an FDR < 0.01, and standardised |log₂FC| > 3 was used in the analysis. The top 20 GO-enriched BPs are shown in Figure 4 [Supplementary Table 3]. The expressed genes correspond to processes like ‘response to stress’, ‘negative regulation of apoptotic process’, ‘cytokine production’, ‘defense response,’ and ‘myeloid leukocyte activation’. However, analysis of genes more expressed in HC than ENL showed GO enriched terms such as – ‘cellular metabolic process’, ‘ion binding’, ‘membrane-bound organelle, and ‘intracellular’. Similarly, the ontologies enriched for genes more expressed in ENL than HC are the’cellular macromolecule metabolic process’, ‘protein metabolic process,’ and ‘heterocyclic compound binding’. The other enriched terms for genes in HC vs ENL comparison are shown in Supplementary Table 3.

Close

Figure 4a:

Dot plots depicting the enriched Gene Ontology (GO) terms. GO terms denote the biological process, cellular component, and molecular function of the expressed genes. The X-axis shows the gene ratio, Y-axis denotes the enriched term. Adjusted P-value < 0.1 is set as a cut-off. The colour scale represents adjusted P-values and the size of the dot represents the gene number. Enriched terms for genes expressed in HC vs ENL comparison.

Export to PPT

Close

Figure 4b:

Dot plots depicting the enriched Gene Ontology (GO) terms. GO terms denote the biological process, cellular component, and molecular function of the expressed genes. The X-axis shows the gene ratio, Y-axis denotes the enriched term. Adjusted P-value < 0.1 is set as a cut-off. The colour scale represents adjusted P-values and the size of the dot represents the gene number. Enriched terms from LL vs ENL comparison, ATP - Adenosine triphosphate.

Export to PPT

Supplementary Table 3

Download

For LL vs ENL comparison, overall 334 genes were expressed at a cut-off of standardised |log₂FC| > 2 and FDR < 0.001. However, only 41 genes were found to be significant by all meta-analysis methods. Therefore, these 41 genes were used to understand their BPs. Functional enrichment analysis reveals the GO enriched BPs such as ‘cell cycle’, ‘regulation of ion transport’, ‘inflammatory response,’ and ‘response to endogenous stimulus’ which are indicative of host response against a pathogen. The top 20 enriched ontologies for genes differentially expressed in LL vs ENL are shown in Figure 4 [Supplementary Table 4]. Genes specific to a particular GO enriched term for HC vs ENL and LL vs ENL comparison sets are shown in Figure 4. Heat plots showing the enriched terms for identified genes based on standardised log₂FC for comparison HC vs ENL and LL vs ENL are shown in Figure 5. Figure 6 represents an overview of the meta-analysis of the datasets representing the two types of leprosy reactions: Type-I (Reversal Reactions) and Type-II (Erythema Nodosum Leprosum).

Supplementary Table 4

Download

Close

Figure 5a:

Heat plots showing enriched genes with specific GO terms according to the standardised log₂FC for HC vs ENL comparison; positive log₂FC denotes genes more expressed in HC individuals while negative log₂FC denotes genes more expressed in ENL patients.

Export to PPT

Close

Figure 5b:

Heat plots showing enriched genes with specific GO terms according to the standardised log₂FC for LL vs ENL comparison; positive log₂FC denotes genes more expressed in LL while negative log₂FC denotes genes more expressed in ENL patients.

Export to PPT

Close

Figure 6:

The overview of the meta-analysis of the datasets representing Type-I and Type-II leprosy reactions. RR- Reversal reaction, DEGs - Differentially expressed genes, ENL - Erythema nodosum leprosum.

Export to PPT

We have identified several DEGs more expressed in RR than in HC such as ABLIM1, ANK3, FAM110C, and MORC4 while those less expressed in RR as compared to HC like BCL2A1, IFI30, and PLAUR and they correspond to GO enriched terms like cellular aldehyde metabolic process and cell surface receptor signalling pathway [Supplementary Table 5]. The genes that resulted from LL vs RR, namely, APOC1, CSTB, FTL, SLC11A1, and TSPAN10 are more expressed in RR than LL and not described in the literature earlier. ADAMTS4, CCL17, CD1B, CD1E, and G0S2 among others are more expressed in LL than RR [Supplementary Table 6]. A major class of genes expressed in LL belongs to chemokines and plays a significant role in recruiting neutrophils and lymphocytes to the mycobacterium-infected site.²³ The GO terms for LL vs RR comparison are immune system process cell communication, and signal transduction which are representative of the need for the cell to generate an immune response and would be essential to understand the immunopathogenesis of the disease. Gene ontology- enriched terms for DEGs in LL vs RR and HC vs RR are shown in Supplementary Table 5 and Supplementary Table 6, respectively. CCL17 and CYP27B1 seem to be more expressed in LL patients and are also reported to have no significant correlation with Type-I reaction ²⁴ while MERTK showed higher expression in RR patients [Supplementary Figure 1].

Supplementary Table 5

Download

Supplementary Table 6

Download

Supplementary Figure 1

Download

Discussion

Leprosy reactions play a major role in morbidities associated with the disease. These immune reactions lead to nerve damage, and weakness and pose a great risk of injury and deformity.² These immune-mediated complications occur in more than 50% of leprosy patients.⁶ The pathological mechanisms of these immune reactions are complex and have not yet been completely understood. Comparative host transcriptomics has immense utility in elucidating screening profiles and is capable of efficiently narrowing down the candidate genes associated with multifactorial diseases like leprosy and investigates the underlying mechanisms of immune reactions for their diagnosis and prognosis. In this study, publicly available datasets from GEO were considered for analysis. The datasets specific to leprosy Type-I and Type-II reactions were identified based on search terms. Individual analysis of these datasets revealed statistically significant DEGs. Integration of these individual studies based on a common comparison pattern was carried out using meta-analysis statistical methods to find common expression patterns across studies. Datasets that resulted from both the microarray and RNA sequencing platforms were included to ensure that no specific dataset is left out merely because of variation in the experimental platform being used. At last, we identified a list of DEGs in each comparison set across different studies which are indicative of their possible role in progression and immune exacerbations in leprosy. Therefore, such common genes can be prioritised as their functional validation would be useful because they are consistently identified across multiple studies.

The criteria used for identifying the differential expression are based on adjusted P-value and standardised fold-change or effect size.²⁵ However, the hypothesis testing depends on the sample size, and studies with less sample size do not or had fewer differentially expressed genes, and thus chances of getting false-positives and false-negatives are expected. When a gene is identified as differentially expressed in one study, it might not necessarily be statistically significant in other unless that gene is biologically replicated. If a gene is selected by the major meta-analysis methods, such as maxP, SR, Fisher, and sdef, then the chances of it being a false-positive are less likely. Also, if a gene is statistically significant in one study, it might not be necessarily picked up by all methods in a meta-analysis by virtue of differences in P-values or effect size as it is not consistent across studies. This corresponds to the fact that only a small proportion of the genes were selected by these methods over those that only differentially expressed in one or few studies.

Only studying the final set of DEGs from individual studies might limit the discoverability of other potential genes as certain factors contribute to making it biased. For example, the statistical framework used in an individual study might differ; all genes from one study might not be expressed in another on account of different platforms in use, etc. Thus, analysis of individual studies removes these biases and helps to incorporate a common statistical method for them to reach a balanced comparable conclusion.

The meta-analysis methods used in this study are based on the P-value or effect size of expressed genes and thus help to reveal DEGs not identified earlier.²⁶ The potential role of these genes can be identified using functional analysis or previously defined pathways relevant to the host response to leprosy. Some of the identified genes have been previously reported in studies relevant to leprosy like PTX3 which showed a higher expression in ENL,²⁷ S100A9 which is more expressed in HC than in the ENL cases and generates an ineffective response to leprosy.²⁸ GO enrichment analysis of lists of DEGs from HC vs ENL reveals the BPs involved in the inflammatory response, cell population proliferation, and cytokine production. Moreover, a similar profile is seen in GO-enriched BPs from LL vs ENL comparison. It represents terms like inflammatory response, membrane organisation, and regulation of transmembrane transport which connotes genes such as SLC11A1 and showed higher expression in ENL patients.²⁹

This study unleashes several potential biomarker candidate genes that may serve the purpose of early detection of reactions in leprosy patients. Over the years, several genes were associated with ENL, for example, Pentraxin-3 (PTX3) which is a protein present in secondary neutrophilic granules and is a potential biomarker of ENL. According to this study, it shows a downregulation in expression from its normal levels. PTX3 can be used to differentiate between episodes of ENL from a reversal reaction. Another gene SLC11A1, a member of the solute carrier family, encodes a protein called natural resistance–associated macrophage protein 1(NRAMP1 or SNP 274C/T) which is associated with leprosy reaction. The presence of the ‘C’ allele becomes a risk factor for Type-I reaction; however, it becomes protective for Type-II reaction. A class of interferons was also identified as a part of candidate biomarkers which constitute a range of cytokines performing multifaceted roles. In this study, IFITM2 and IFITM3 are shown to be upregulated in ENL as also confirmed by other studies.

The gene expression can be validated by techniques such as the qPCR (quantitative polymerase chain reaction). It is quite interesting to note that the analysis of several microarray and RNA sequencing datasets by different people generating analogous results is indicative of a true differential expression. However, proper replication of these findings would be required to rule out any bias.³⁰

Conclusion

This study analyses immune reactions in leprosy-specific datasets publicly available in GEO using a standardised approach, thereby revealing the DEGs common among different studies. Integration of DEGs from individual studies using meta-analysis reveals several genes that are vital in understanding the mechanisms of immune reactions and their possible correlation with the host response. We have identified several genes ADAMTS5, ADAMTS9, IFITM2, IFITM3, KIRREL, ANK3, CD1E, CTSF, DOCK9, and KRT73 to name a few which can serve as potential biomarkers for Type-II reaction. For Type-I reaction, DEGs are identified as ACP5, APOC1, CCL17, S100B, and SLC11A1 among others, and may likely serve as biomarkers. Further, the analysed studies imply some important BPs which help in getting insights into the pathogenesis of reactions in leprosy. A combination of biomarkers would be necessary for its early diagnosis or prognosis to be used for clinical applications.