A Computational Study Identifies HIV Progression-Related Genes Using mRMR and Shortest Path Tracing

Since statistical relationships between HIV load and CD4+ T cell loss have been demonstrated to be weak, searching for host factors contributing to the pathogenesis of HIV infection becomes a key point for both understanding the disease pathology and developing treatments. We applied Maximum Relevance Minimum Redundancy (mRMR) algorithm to a set of microarray data generated from the CD4+ T cells of viremic non-progressors (VNPs) and rapid progressors (RPs) to identify host factors associated with the different responses to HIV infection. Using mRMR algorithm, 147 gene had been identified. Furthermore, we constructed a weighted molecular interaction network with the existing protein-protein interaction data from STRING database and identified 1331 genes on the shortest-paths among the genes identified with mRMR. Functional analysis shows that the functions relating to apoptosis play important roles during the pathogenesis of HIV infection. These results bring new insights of understanding HIV progression.


Introduction
Many efforts have been devoted to better understanding the mechanism governing disease progression and non-progression during HIV infection. Besides the direct cytotoxic effect against CD4+ T cells caused by HIV, immune activation is widely accepted as a good predictor for disease progression [1,2,3,4,5]. Furthermore, some clinical studies using immunity suppressive drugs to suppress immune activation slowed down the diseases progression [6,7]. However, the molecular mechanism underlying the immunopathogenesis remains obscure. To identify the host factors important for the HIV-1 pathogenesis and disease progression, high throughput techniques had been employed. Genome wide association studies revealed the protective effect against the virus of human leukocyte antigens (HLAs) including HLA-B*57:01, B*27:05 and risk alleles including HLA-B*35, Cw*07 [8,9]. These studies further led to the finding of the protective effect of HLA-C [10] Transcriptome studies also gained important insights, regarding to interferon stimulated genes (ISGs), immune activation, cell cycle and cell death during the infection. People had conducted transcriptome researches for identifying factors affect the viral control and the speed of CD4+ T cell loss [11,12]. A research with 137 HIV seroconverters, 16 elite controllers and 3 healthy blood donors attempted to identify some molecular factors associated with the viral control. More surprisingly, successful treatment made the transcriptome states of patients similar to the elite controllers and the HIV-negative donors [11]. Another study compared the transcriptoms of 6 viremic non-progressors (VNPs) and more than 20 rapid progressors (RPs). No significant result was found. Genes identified from the data of monkeys (CASP1, CD38, LAG3, SOCS1, EEIFD, and TNFSF13B) were deemed as the factors affecting the speed of disease progression [12].
Machine learning was recently proved to be an effective strategy for accurate classification of phenotypes based on transcriptome data (gene expression microarray) [13,14,15,16,17]. Among them, minimum redundancy -maximum relevance method (mRMR) is robust and represents a broad spectrum of characteristics [18,19]. It was also developed to identify disease-related genes from expression profiles [18,19].
Another useful informatics strategy for disease candidate gene identification is by known protein-protein interactions (PPIs). Since proteins not only function individually by themselves, but also co-function with their interaction partners; thus interaction partners of disease related genes are also important candidates for further disease casual studies. The STRING (Search Tool for the Retrieval of Interacting Genes) database is an online resource that provides PPI information by reporting from both prediction and experimental observations [20].
Here, we present a comprehensive informatics study based on transcriptional profiling of three different groups of HIV patientsrapid progressors (RPs), viremic controllers (ECs) and viremic nonprogressors (VNPs). We attempted to i) identify a gene set which can well classify the three groups, by using mRMR feature selection; ii) provide candidate casual genes for further experimental studies, by using shortest-path analysis of the above identified genes in a molecular interaction network contructed with the STRING data.

Gene expression profiling dataset of HIV patients
The dataset was from a research on HIV infection done by Rotger et al. [12]. In total, 78 chips were used in that research. We adopted data generated from CD4+ T cells, which contains 40 microarrays (8 elite controllers (ECs), 27 rapid progressors (RPs), 5 viremic non-progressors (VNPs)). Using the dataset alone, Rotger et al., didn't observe any differentially expressed genes. The data was downloaded from NCBI Gene Expression Omnibus (GEO) with the accession number of GSE28128. The expression profile was generated using the microarray Illumina HumanWG-6 v3.0 expression beadchip. Bead summary data was the output from Illumina's BeadStudio software without background correction. Genes declared as non-expressed (P.0.01) were excluded from further analysis. Data preprocessing, including quantile normalization and log2 transformation was completed in the Partek Genomics Suite package (Partek Inc.).

Minimum redundancy -maximum relevance algorithm
Minimum redundancy -maximum relevance algorithm for selecting features (genes) was developed based on the idea to balance features' 'relevance' to target (phenotype) and 'redundancy' between features [18]. Both relevance and redundancy are quantified using mutual information (MI). In this study, mRMR was realized using a R package 'mRMRe' [21], in which MI is estimated as, where I and r represent the MI and the correlation coefficient between variables x and y, respectively. Let y and X = {x 1 , …, x n } be the input variable (phenotype) and set of input features (genes), respectively. Given x i as the feature with highest MI with the output variable, so the set of features, denoted by S is then initialized with x i .
In the second step, the feature x j with the best balance between highest relevance and lowest redundancy was added to S. It is achieved by maximizing the score q as follows, The second step is repeat until a desired solution length has been reached.

Prediction engine
In our analysis, predicted phenotype of an individual was estimated in three ways: the phenotype of its nearest neighbor; the most phenotype of its 5 nearest neighbor; the phenotype of its nearest clustering center of each phenotype group (Table S1). Distance between two individuals were calculated according to Chou et al.'s studies [22,23], as follows, where d refers to the distance, i 1 and i 2 represent two samples, and e 1 and e 2 are vectors of selected features (expression levels of selected genes) of the corresponding individual.

Validation and Incremental Feature Selection
Jackknife validation was used to estimate the prediction accuracy of the selected features, considering its better performance and applicability to the dataset compared with other validation methods such as independent data set test and subsampling test [22,24,25]. Its idea is as follows. Given x samples with known outcome variable (phenotype) and n selected features (expression levels of genes), for each sample, we compared its known outcome variable with an estimated outcome variable, and the estimation was performed based on the rest x -1 samples. The accuracy is defined using the following equation, where TP, TN, FP and FN represents the numbers of true positives, true negatives, false positives and false negatives, separately. Incremental Feature Selection (IFS) was used to determine the number of prediction features. For n = 1 to 800 required number of features, each feature set was determined using mRMR and its prediction accuracy was estimated using Jackknife validation. The set with best prediction accuracy and smallest feature number was regard as final feature set. In this study, a set with 147 genes was chosen and its prediction accuracy is 0.8049.

Shortest path identification in STRING PPI network
The background weighted PPI network was constructed using data from STRING database (version 9.1) (http://string-db.org) [20]. Weights of edges in the PPI network present confidence of the PPI ranging from 1 to 999. The Dijkstra's algorithm implemented in a R package 'igraph' [26,27], was applied to identify shortest path between two pairs of proteins, and each of the protein corresponds to 86 protein coding genes of 147 mRMR-IFS identified genes. And the network of shortest paths was constructed using Cytoscape (version 3.0.1) [28].

Pathway enrichment analysis
KEGG pathway enrichment and GO functional enrichment analysis were carried out using web service of DAVID tools (version 6.7) [29,30]. Original enrichment P values as well as Benjamin multiple test corrected P values were estimated.

A set of 147 genes presents the best accuracy for prediction of RPs, ECs and VNPs
Based on the outputs of mRMR, we tested the predictor with one feature (probes of gene expression array), two features, three features, etc., and the IFS result is provided in Figure 1. In the IFS curve, the X-axis is the number of probes used for classification, and the Y-axis is the prediction accuracies of the nearest neighbor algorithm evaluated by the Jackknife validation. The accuracy reaches its maximum value when 147 features were included, corresponding to 147 different genes (Table 1 for 20 genes  identified with top mRMR scores and Table S2 for the full list). We discussed reported function of three genes in Table 1, SRP14P1, SLC45A2, and DNAJB1. Figure 2 A to C show their expression difference between VNP and RP. Although fold changes are not dramatic, there are significant difference between them (P value, 0.09957 to 1.636*10 25 ).

A PPI sub-network provides insights for the 147 genes on HIV infection
Furthermore, we constructed an undirected graph with the PPI data from STRING [20]. Then we picked all pairs of any two genes from 86 protein coding genes of the 147 genes identified with mRMR as described above, and revealed the shortest path between these two proteins using the Dijkstra's algorithm [26]. We eventually obtained a total of a sub-network of STRING PPIs based on the shortest paths ( Figure 3). There are a total of 4248 protein-protein interactions of 1331 proteins. Among the 1331 proteins, 1290 of their corresponding genes were annotated in the Ensemble Biomart database, and we ranked these genes according to their betweennesses ( Table 2 for the top 20 and Table S3 for the full list). Among these 1290 genes, UBE2K has the largest betweenness of 66, meaning that there are at least 33 shortest paths going through this gene. Accordingly, UBE2K may play an important role in connecting the 66 candidate genes and hence may be related to the loss of CD4+ T cells, although we didn't find any previous reports about the effect of UBE2K on the loss of CD4+ T cells.

Function and pathway enrichment of the identified genes
Using the functional annotation tool DAVID [29,30], the Gene Ontology (GO) and KEGG pathway enrichment analysis were carried out for the 147 genes identified using mRMR as well as the 1290 genes on the shortest paths (Table S4). The results showed that these genes were significantly enriched in the GO modules of protein localization in organelle and SMAD protein nuclear translocation GO modules. In KEGG pathway enrichment analysis, no pathway was significantly enriched. These results suggested that protein location may play critical role in the pathogenesis during the chronic HIV infection.
We also analyzed the genes on the shortest paths among the genes identified by mRMR algorithm with GO and KEGG enrichment. In the GO enrichment analysis, many genes were significantly enriched in GO modules relating to negative regulation of apoptosis (Table S5). Early studies had shown the importance of apoptosis during the pathogenesis of HIV infection [31,32]. Vulnerability to apoptosis is a hallmark character of CD4+ T cells due to immune activation and Tat protein [31,33]. Our results supported that resistance to apoptosis might be an important mechanism of the VNPs to maintain their CD4+ T cell levels. In the KEGG enrichment results, nothing strongly relating to the high CD4+ T cell level was discovered (Table S6).

Using mRMR to identify candidate gene set
We applied mRMR algorithm to study HIV progression related genes based on transcriptome data. The major difference between mRMR and another machine learning method is that mRMR considers redundancies among features as well as relevance between feature and target. Therefore, disease-related genes which shown less correlation in expression among themselves, are favored by the mRMR algorithm. This may be the theoretical reason that genes identified in our study represent various functional groups in GO annotations (Table S4).

Result validation
In this study, we used a jackknife method to validate our result. In jackknife validation, one sample x i is excluded from all the n samples X = {x i |i = 1,2,3 … n}. We predict its phenotype using the rest n-1 samples, and compared its predicted phenotype with the real one. The above step is repeated for all x i , i = 1,2,3 … n. In this way, prediction accuracy of X is considered as a parameter of validation. Previous studies also discussed the advantages of jackknife validation over another validation methods, such as subsampling test and independent data set test [22,24,25]. Jackknife validation can exclude ''memory'' effect as well as  arbitrariness problem existed using the other two validation methods.
Although jackknife validation was proved to be useful, we also considered using an independent dataset of general HIV infection as an additional validation to our results. The independent dataset is composed of expression profiles of CD4 T cells from 127 HIV untreated samples and 8 normal controls. Using student t test, we observed that 46 genes out of the 147 mRMR selected genes differentially expressed under P value 0.05 [11]. We included the related results in Table S2.
In fact, polymorphisms (e.g. SNPs) of human SLC transporter genes have been found to be closely related to HIV infection or responses to treatments of HIV-infected patients. For example, the genotypes of SLC11A1, a proton-coupled metal ion transporter, have determining effects on mortality in HIV infection [35]; SNPs of loci in or near SLC26A7, a multifunctional anion exchanger linked to energy metabolism and immunoregulations, are strongly associated with upper trunk and arm subcutaneous adipose tissue (SAT) distribution in antiretroviral therapy (ARV)-treated HIV-  infected patients [36]. Furthermore, literatures suggest that antiretroviral drugs such as efavirenz and nevirapine are substrates of a variety of SLCs and SNPs of them may influence highly active antiretroviral therapy (HAART) efficacy and AIDS-free survival [37,38]. Thus it is no surprise that the SLC genes included in our results of mRMR (4A8, 22A10, 22A15 -inorganic/organic ions transporters; 25A1 -mitochondrial carrier; 28A1 -nucleoside transporter; 45A2 -putative sugar transporter), which have similar functions to those in the literatures, may sustain similar functional roles and serve as good surrogates for prediction of HIV/AIDS outcomes.
Another gene among the 147 genes we noticed is Mouse Double Minute 2 Homolog (MDM2) (mRMR score = 0.1185). It is also known as E3 ubiquitin-protein ligase Mdm2, which is a negative regulator of p53 tumor suppressor [39]. Mdm2 inhibits p53 cellcycle arrest and apoptic functions and the interaction with Mdm2 can also result in a large reduction in p53 protein levels through enhanced proteasome-dependent degradation. Endogenous levels of Mdm2 are sufficient to regulate p53 stability, and overexpression of Mdm2 can reduce the amount of endogenous p53. Mdm2 is also found to be important for lymphopoiesis through the inhibition of p53 [40], which makes it a potential factor regulating the CD4+ T cell counts during chronic HIV infection.
The importance of gene DNAJB1 (DnaJ homolog subfamily B member 1, mRMR score = 0.1717,) to HIV/AIDS is self-evident, as experimental studies show that it is regulated by various subtypes of HIV (e.g. HIV-1 B, C and A/E) in dendritic cells (DCs), which are among the first targets of HIV infection and in turn play crucial roles in viral transmission to T cells and regulation of immune responses [41,42]. Meanwhile, its gene product DNAJB1 (Hsp40), is essential for HIV-1 Nef-mediated enhancement of viral gene expression and replication, thus playing a key role in the virus life cycle [43,44]. SRP14P1 (scored the highest 0.374 by mRMR) is named as Signal Recognition Particle 14 kDa (Homologous Alu RNA Binding Protein) Pseudogene 1, which is rarely studied. SRP14 is an Alu RNA binding protein. It regulates the expression of Alu elements, which is an important expression regulation element. This finding in our work suggested its potential role in affecting the pathogenesis of HIV infection.
Interferon response is of important function in HIV infections. But we didn't observe interferon stimulated genes. As suggested by Rotger et al. [11], during the course of HIV infection the interferon stimulated genes are stimulated by the virus. So it may explain that in our case, holding the viral load equal, the interferon stimulated genes are not identified.

Conclusion
In the previous study [12], M. Rotger et al. endeavored to explore some molecular factors relating to the maintenance of CD4+ T cell level using the transcriptome derived from CD4+ T cells of VNPs and RPs. The results were not encouraging. Therefore, we made the attempt to discover something new from a different perspective. Although most genes found in our results did not closely relate to the different CD4+ T cell levels between the two groups of people, functions associated with apoptosis regulation had been identified. These results reflected the complexity of the mechanism governing the decline of CD4+ T cell count during HIV infection. To some extent, this could also be explained by the design of the study. The data were generated from a retrospective study, thus the factors mediating the decline of CD4+ T cells might function early among the RPs and returned to normal since the CD4+ T cells had already dropped to a low level. Taking that into consideration, perspective studies for further discoveries could become feasible. Supporting Information Table S1 Predictions of phenotype of each individuals using the 147 genes identified using mRMR. (TXT)

Author Contributions
Conceived and designed the experiments: LL CM XD RL. Analyzed the data: CM XD. Contributed reagents/materials/analysis tools: CM. Wrote the paper: CM XD RL.