Identifying candidate diagnostic markers for tuberculosis: A critical role of co-expression and pathway analysis

1.   Introduction

Tuberculosis, caused by Mycobacterium tuberculosis, has been posing a profound burden on global public health ^[1,2]. Despite successful discoveries of several anti-mycobacterial drugs, the results of the applied chemotherapy are far from satisfactory. Treatment requires many months of taking several drugs, many of which have side effects. This, coupled with the complicated infectious processes and high transmissibility of the disease in the population, has resulted in widespread emergence of drug-resistant and multidrug-resistant strains of Mycobacterium tuberculosis ^[3].

Most individuals with latent tuberculosis infection do not present symptoms, so detection of tuberculosis infection has been a real challenge. Some diagnosis methods currently in use have limitations including lack of high degree of sensitivity, false positive tendency, and high cost ^[4,5]. As tuberculosis infection is the outcome of a complicated host-pathogen interaction, identifying host factors involved in this interaction is key to better infection control and effective diagnostic biomarkers. In particular, exploration of key genes associated with tuberculosis is important for diagnosis development and pathogenesis understanding. There are, however, only a handful of host genes confirmed to participate in the aforementioned interaction ^[6].

Gene expression database GEO based on gene chip technology provides a large number of information related to the disease. Here we focus on a particular set of gene expression data associated with W-Beijing families of Mycobacterium tuberculosis strains. We are interested in this data set of 40 samples because the data set contains detailed information about test time after infection by tuberculosis strains. This information, as our study will show, is critical for the examination of the relationship between tuberculosis development and infection time. Using a multidimensional scaling analysis, we will observe the significantly different gene expressions at different time stages. Applying the linear models in limma package and empirical Bayes moderation in Bioconductor of R to contrasting the gene expression in different infection time stages and performing gene set testing, we will discover more host genes associated with tuberculosis. To verify these key genes from the underlying biological mechanism and to explore the corresponding pathogenicity process, we will generate a gene co-expression network map based on the mutual regulation between the discovered significant genes and the gene regulatory network associated with tuberculosis for the interaction of the significant genes based on gene regulations. We will also explore the gene pathway figures based on our discovered significant genes, providing further insights about the pathogenesis of tuberculosis.

2.   Materials and methods

2.1.   Data

This data was obtained through NCBI (National Center for Biotechnology: https://www.ncbi.nlm.nih.gov/geo/download/?acc=GDS4258). This includes the gene expression of more than 50,000 probes for 39 THP1 cell samples infected by 12 different W-Beijing strains of tuberculosis and one control sample. The most significant feature of the data set is the important information about infection times since the samples were tested after 4 hours, 18 hours and 48 hours of infection for each strain (see Table. 1).

Table 1.  Data description.

Strains	Testing Time			Row Sum
	After 4 hours	After 18 hours	After 48 hours
R1.4	1	1	1	3
R17.1	1	1	1	3
ZA9.2	1	1	1	3
ZA9.4	1	1	1	3
R19.4	1	1	1	3
CHN50.1	1	1	1	3
MAD2.1	1	1	1	3
CHN50.2	1	1	1	3
R17.3	1	1	1	3
R19.5	1	1	1	3
H37Rv	1	1	1	3
MAD2.2	1	1	1	3
Column Sum	13	13	13	39

---

 | Show Table

DownLoad: CSV

2.2.   Data pre-processing

The raw counts data and normalized data are both available in GDS4258. The data after RPKM (reads per kilobase of transcript per Million mapped reads) normalization was used in this study. Altogether, there were 54675 gene probes. Combination of the duplicate gene symbols was done by adopting the average. The mean-variance relationship was plotted to check the level of (biological) variation in the experiments and the level of filtering performed upstream. Usually, experiments with high variation result in flatter trends, where variance values plateau at high expression values. Experiments with low variation tend to result in sharp decreasing trends. If filtering of lowly-expressed genes is insufficient, a drop in variance levels can be observed at the low end of the expression scale due to very small counts.

2.3.   Unsupervised clustering of samples

For this study, it is important to see which genes are expressed at different levels between the different sample classifications profiled. Therefore, classification of the samples is an essential step in our bioinformatics analysis. In particular, only on the basis of good classification, contrast could be done among different classifications and differentially expressed genes could be discovered effectively. We used the multidimensional scaling (MDS) unsupervised clustering method, a means of visualizing the level of similarity of individual cases of a dataset. MDS is a form of non-linear dimensionality reduction, and the MDS plot shows similarities and dissimilarities between samples in an unsupervised manner so that one can have an idea of the extent to which differential expressions can be detected before carrying out formal tests. Ideally, samples would cluster well within the primary condition of interest, and any sample straying far from its group could be identified and followed up for sources of error or extra variation. In our dataset, time tested and virus strains are two possible classification criterions and were therefore examined.

2.4.   Differential expression analysis

Our analysis followed the workflow in Bioconductor which carried out linear modelling in limma package and fitted to the data with the assumption that the underlying data was normally distributed. To initiate, we set up a design matrix with the classification information, i.e., we set up the contrasts for pairwise comparisons between the classifications in limma using the makeContrasts function. We then carried out Empirical Bayes moderation by borrowing information across all the genes to obtain more precise estimates of gene-wise variability ^[7].

Differential expression levels were defined using an adjusted p-value cutoff, that was set at $5\%$ by default. The number of significantly up- and down-regulated genes can be summarised. In this study, for a stricter definition on significance, log-fold-changes (log-FCs) to be above a minimum value was applied where the treat method ^[8] was used to calculate p-values from empirical Bayes moderated t-statistics with a minimum log-FC requirement.

We generated mean-difference plots for the contrasts of different classifications which display log-FCs from the linear model fitting against the average log-CPM values with the differentially expressed genes highlighted. We then created a heat map for the top differentially expressed genes (as ranked by adjusted p-value) from one of the contrast and this allowed us to look at the expression of a subset of genes. This process can give useful insights into the expression of individual groups and samples without losing perspective of the overall study when focusing on individual genes, or losing resolution when examining patterns averaged over thousands of genes at the same time. Our heat map correctly clustered samples into sample type and rearranged the order of genes to form blocks of similar expressions. Finally, genes that were differentially expressed (up-or down-regulated) in multiple comparisons were extracted as the most significant genes.

2.5.   Gene co-expression analysis

Based on the mutual regulation between significant genes selected, we generated a gene network map with the gene regulatory network associated with tuberculosis through Genevestigator (https://genevestigator.com/gv/), a powerful tools search engine for gene expression with advanced analysis possibilities. Genevestigator includes the search for genes that are specifically expressed under certain conditions and the search for groups of genes sharing similar expression patterns by means of clustering and biclustering algorithms. In this study, the co-expression function of Genevestigator was used to find co-regulated genes with a set of differential expression genes of interest. P-value 0.05 and fold-change 2 or 0.5 were selected as criteria. We also used the co-expression networks database Coexpedia (http://www.coexpedia.org/) to search the tuberculosis related co-expression genes.

2.6.   Gene pathway analysis

We explored the gene pathway analysis based on the human gene database GeneCards (http://www.genecards.org/) through the software Pathway Builder Tool 2.0. This exploration detected some relationships between significant genes and the pathogenesis of tuberculosis.

3.   Results

3.1.   Unsupervised clustering of samples

After combing the duplicated genes by adopting average, there were 29822 gene symbols left. The model's residual variances analysis against average expression (Figure 1 in supplementary file) showed that the variance was not dependent on the mean expression level, excluding the high biological variation of the experiment.

Figure 1.  Classification of samples.

DownLoad: Full-Size Img PowerPoint

In this dataset, 39 samples clustered well within test time over dimension 1 and 2 (Figure 1). The first dimension of a data point represents the leading-fold-change that best separates samples, and explains that the largest proportion of variation in the data with subsequent dimensions, orthogonal to each other, has smaller effects. Dimensions 3 and 4 were also examined using the grouping defined by strains and the classification is not good.

Whilst all samples cluster by time, the largest transcriptional difference was observed between test time after 4 hours (t4) and test time after 18 hours (t18), and t4 and t48 (test time after 18 hours) over dimension 1. For this reason, it was expected that pairwise comparisons between samples will result in a greater number of differentially expressed genes for comparisons involving t4.

3.2.   Differential expression analysis

To carry out the linear models in limma packag and discover the differentially expressed genes, we set up a design matrix with the classification information. Here we chose the model that removes the intercept from the first factor and test time, but keeps the intercept in the second factor strains for our analysis, since setting up model contrasts became more straightforward in the absence of an intercept for test time. Contrasts for pairwise comparisons between the classifications were set up in limma using the makeContrasts function and the main contrasts were between the different test time groups. Based on the empirical Bayes moderation, the number of significantly up- and down-regulated genes was obtained. The number of differentially expressed was found to be 10581 (including both up- and down-regulated genes) for t4 versus t18, 13379 genes for t4 versus t48, and 9614 genes for t18 versus t48 (Table 2, left).However, after the treat method was used and a minimum log-FC was considered, the number of differentially expressed genes greatly diminished to a total of 493 genes for t4 versus t18,738 genes for t4 versus t48, and 117 genes for t18 versus t48 when testing required genes to have a log-FC that was significantly greater than 1(Table 2, right). The integration of the differentially expressed genes in different contrasts allowed us to extract a total of 287 genes which are differentially expressed in both t4 versus t18 and t4 versus t48 (Figure 2). These were taken as most significant genes associated with tuberculosis (Table 1 in supplementary file).

Table 2.  Number of up-and down-regulated genes for empirical Bayes (left) and treat method (right).

	t4 vs t18	t4 vs t48	t18 vs t48			t4 vs t18	t4 vs t48	t18 vs t48
Down	4408	6351	4977		Down	393	431	62
Notsig	19241	16443	20208		Notsig	29329	29084	29705
Up	100	307	55		Up	6173	7028	4637

---

 | Show Table

DownLoad: CSV

Figure 2.  Differentially expressed genes in both t4 vs t18 and t4 vs t48.

DownLoad: Full-Size Img PowerPoint

3.3.   Gene co-expression analysis

We built a gene network map based on the 287 differentially expressed genes where 212 genes formed the network. All co-expression links were evaluated for functional association by statistical assessment (left panel of Figure 3). Moreover, we noted a sub-network (the yellow part of the left panel of Figure 3) with 27 genes that links according to Medical Subject Headings (MeSH, the National Library of Medicine's controlled vocabulary thesaurus) terms and associated with tuberculosis (right panel of Figure 3). The links among the sub-network have direct evidence of anatomical or disease context information (Table 3).

Figure 3.  Gene network map.

DownLoad: Full-Size Img PowerPoint

Table 3.  The rank of the genes in the sub-network.

Rank	Gene Symbol	Gene Name	Score
1	RSAD2	radical S-adenosyl methionine domain containing 2	41.5
2	IFI44	interferon induced protein 44	35.46
3	CMPK2	cytidine monophosphate (UMP-CMP) kinase 2, mitochondrial	30.94
4	IFI44L	interferon induced protein 44 like	30.51
5	IFIT3	interferon induced protein with tetratricopeptide repeats 3	27.41
6	XAF1	XIAP associated factor 1	27
7	IFIT1	interferon induced protein with tetratricopeptide repeats 1	24.7
8	IFIH1	interferon induced, with helicase C domain 1	21.99
9	GBP1	guanylate binding protein 1, interferon-inducible	16.91
10	IFIT5	interferon induced protein with tetratricopeptide repeats 5	14.34
11	HERC5	HECT and RLD domain containing E3 ubiquitin protein ligase 5	14.17
12	SAMD9L	sterile alpha motif domain containing 9-like	11.3
13	STAT1	signal transducer and activator of transcription 1	10.16
14	ISG15	ISG15 ubiquitin-like modifier	9.15
15	OAS1	2'-5'-oligoadenylate synthetase 1	7.29
16	GBP5	guanylate binding protein 5	6.57
17	IFIT2	interferon induced protein with tetratricopeptide repeats 2	5.29
18	RTP4	receptor (chemosensory) transporter protein 4	5.28
19	RPLP2	ribosomal protein, large, P2	3.29
20	RPL31	ribosomal protein L31	3.29
21	UBE2L6	ubiquitin conjugating enzyme E2L 6	2.7
22	CASP1	caspase 1	2.05
23	TNFSF10	tumor necrosis factor superfamily member 10	2.05
24	MX1	MX dynamin-like GTPase 1	1.91
25	SERPING1	serpin peptidase inhibitor, clade G (C1 inhibitor), member 1	1.4
26	IFITM2	interferon induced transmembrane protein 2	1.25
27	IFITM1	interferon induced transmembrane protein 1	1.25

---

 | Show Table

DownLoad: CSV

3.4.   Gene pathway analysis

To explore the pathogenesis of tuberculosis, we explored the gene pathway figure based on our identified significant genes. Four gene pathway figures were obtained and 24 genes were found participated in the signalling pathways, among which there are 17 genes participated in the NF-$\kappa$B signalling ^{[9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31]}, 2 genes participated in the JAK-STAT pathway ^[32,33,34], 2 genes participated in the Autophagy pathway ^[35,36] and 8 genes participated in the Apoptosis pathway ^{[37,38,39,40,41,42,43]} (http://pathwaymaps.com/maps/721/) (Figure 4). For the 16 genes in the NF-$\kappa$B signalling, 6 of them were selected in the gene network map associated with tuberculosis. The NF-$\kappa$B signalling and JAK-STAT pathway can accelerate apoptosis and autophagy of host cell ^[44,45], and apoptosis and autophagy are two main ways for macrophages to resistant and remove Mycobacterium tuberculosis so as to maintain the stability of the environment in the body ^[46,47].

Figure 4.  NF-$\kappa$B signalling, JAK-STAT pathway, Autophagy pathway and Apoptosis pathway and their relationship.

DownLoad: Full-Size Img PowerPoint

4.   Discussion

Tuberculosis caused by Mycobacterium tuberculosis has been a serious threat to the public health worldwide. Identifying key genes which express differentially with respect to infection times and construction of co-expression network to reveal the process of gene regulations are important for the detection of tuberculosis.

A major finding in our bioinformatics study of an important data set from NCBI is the identification of some key genes express differentially with respect to different infection times. Our gene co-expression analysis showed that some of these differentially expressed genes have played an important role. Based on differentially expressed genes, our constructed co-expression network revealed the process of gene regulations. We identified hubs, which are topologically central in the co-expression network which have maximal informational connections with other genes. We also identified a sub-network relevant to tuberculosis. We further refined the network to have obtained a more complete network with deferent expressed genes. With this refinement, genes in the original network but beyond the sub-network could be used to assist in predicting new pathway on gene regulation and expression.

We also found that some of the key genes participated in the important pathways such as the NF-$\kappa$B signalling and JAK-STAT pathway which accelerate apoptosis and autophagy of host cells. As Macrophage is well known to be a key cellular immune response to the host's resistance to Mycobacterium tuberculosis, once the organism is infected by Mycobacterium tuberculosis the surface receptors of the macrophage in the host recognize this immediately and use some special ways to fight against tuberculosis (such as apoptosis and autophagy). Therefore, the search for molecules interacting with the host of Mycobacterium tuberculosis is essential to diagnose and control Mycobacterium tuberculosis ^[48], important for the development of new drugs, and useful for the search of new drug targets.

The gene co-expression network might vary as time goes by, examining this variation is an interesting topic for future research. Some genes whose expressions do not vary much over time might also play significant roles in the co-expression network, exploring those representative genes remains to be an important next step of investigation.

Acknowledgments

We are grateful to the three reviewers for their valuable comments and suggestions that have helped us improving the presentation of this paper. This study was supported by the National Science Foundation of China (No.11701471 and 11671327), Basic Science and Frontier Technology Research Project of Chongqing (cstc2017jcyjAX0476), the Canada Research Chair Program and the NSERC/Sanofi Industrial Research Chair program.

Conflict of interest

All authors declare no conflicts of interest in this paper.