THD-Module Extractor: An Application for CEN Module Extraction and Interesting Gene Identification for Alzheimer’s Disease

There exist many tools and methods for construction of co-expression network from gene expression data and for extraction of densely connected gene modules. In this paper, a method is introduced to construct co-expression network and to extract co-expressed modules having high biological significance. The proposed method has been validated on several well known microarray datasets extracted from a diverse set of species, using statistical measures, such as p and q values. The modules obtained in these studies are found to be biologically significant based on Gene Ontology enrichment analysis, pathway analysis, and KEGG enrichment analysis. Further, the method was applied on an Alzheimer’s disease dataset and some interesting genes are found, which have high semantic similarity among them, but are not significantly correlated in terms of expression similarity. Some of these interesting genes, such as MAPT, CASP2, and PSEN2, are linked with important aspects of Alzheimer’s disease, such as dementia, increase cell death, and deposition of amyloid-beta proteins in Alzheimer’s disease brains. The biological pathways associated with Alzheimer’s disease, such as, Wnt signaling, Apoptosis, p53 signaling, and Notch signaling, incorporate these interesting genes. The proposed method is evaluated in regard to existing literature.

Researchers have also introduced several tools and methods to support extraction of network modules from co-expression networks [2][3][4][5] .
Medical applications of analysis of co-expressed modules are many. One of such applications is to find interesting genes related to a disease. Most of the network based approaches extract co-expressed network modules from a CEN. Analysis of such network modules facilitates the study of molecular differences across different stages of a disease 6 . Methods presented in refs 7-11 suggest the application of CEN module extraction techniques in the analysis of genes related to diseases, such as Alzheimer's Disease (AD). These studies consider the co-expressed genes in network modules based on expression similarity.
Gene expression data assess only experimental information, but not the functional associations among genes. However, for biological relevance of results, Gene Ontology (GO) derived information must be incorporated, which encode additional information, such as, functional similarity and co-membership of genes or proteins in a biological pathway 12 . Functional similarity between two genes with annotations can be measured using their semantic similarity. Semantic similarity is based on GO database and describes the topological distance between two terms in the hierarchical taxonomy. Semantic similarity between two genes can be measured on the basis of information content (Ic). Ic is defined as shown in Equation (1).
where, = P t number of annotations involving a GO term t number of genes ( ) ' ' This paper considers the important issue of analysis of CEN using both gene expression similarity and semantic similarity. The proposed method analyzes CEN considering not only the highly co-expressed genes, but also the genes with less expression similarity, yet high semantic similarity, which are termed as border genes in this paper. The border genes obtained are found to be involved in biological pathways, and therefore, are highly functionally related.
The proposed method is applied to find the potential biomarkers related to the progression of AD. Alzheimer' disease, mostly seen in elderly persons, is related to degeneration of neurons, leading towards dementia. With the growth of data from genomics, proteomics, and interactomics, it has become easy for researchers to understand the dysfunctionalities at molecular level during progression of a disease. In case of AD, pathways, such as Wnt signaling, Alzheimer's disease, Apoptosis signaling, and Glycolysis, are related to the progression of the disease. Therefore, understanding of the genes participating in a network module with Pathway analysis, KEGG enrichment analysis, and GO terms helps to identify the biomarkers related to the disease.
Existing methods 9,11,[13][14][15] , use CEN module extraction techniques that prioritize the genes and study them at molecular level during progression of AD and other neurodegenerative diseases. Alzbase 11 is an integrated approach, which reveals links between upstream genetic variations and downstream endo-features and prioritizes genes in relation to AD. Method proposed by Yue et al. 9 combines gene pair scores of four existing methods to provide robust and useful results enriched with pathways, such as AD, Parkinson's disease, and Huntington's disease.
In this paper, a method named THD-Module Extractor is introduced, which can extract biologically relevant modules from CEN and identify border genes related to AD with less expression similarity and high semantic similarity. The method accepts a microarray gene expression data M and two thresholds, namely, expression similarity threshold (δ) and minimum neighborhood threshold (ρ). The method constructs CEN from gene expression data and extracts highly co-expressed modules.
The proposed method uses SSSim 16 as the expression similarity measure and Lin 17 as semantic similarity measure. SSSim finds similarity between genes exhibiting shifting, scaling, and shifting-and-scaling patterns. The SSSim measure is chosen over other correlation measures since SSSim is robust to noisy values. For semantic similarity, the Lin measure is chosen as it is a normalized relative measure that finds the differences in information content of two GO terms being compared 12 .
The effectiveness of the proposed method is validated using various datasets, including AD dataset. In case of AD dataset, the method identifies the border genes involved in pathways related to the progression of the disease. It is observed that the highly differentially expressed border genes share some important pathways, GO terms, and KEGG pathways, known to be related to AD. The border genes are assessed to find their significance in progression of AD. The work flow of THD-Module Extractor, is illustrated in Fig. 1.
This work provides a comprehensive solution towards identifying the biomarkers of diseases, such as AD, based on expression and semantic similarity of genes, obtained from module extraction method. The following are the major contributions of this work.
1. In this paper, a CEN module extraction method named THD-Module Extractor is introduced, which identifies not only the highly co-expressed genes, but also extracts the genes with less expression similarity and high semantic similarity. In contrast, these genes are not considered by any of the existing related methods. 2. The method is validated using datasets pertaining to various species, in terms of both statistical and biological significance measures. For the AD dataset, the genes with less expression similarity and high semantic similarity are assessed using GO, pathway analysis, and KEGG enrichment analysis to prioritize the genes related to the progression of the disease. 3. Correlation analysis of each pair of genes in a large dataset is computationally expensive. This work leverages the parallel computing capabilities of the Graphics Processing Units (GPU) to find the SSSim 16 correlation matrix, implemented using the NVIDIA CUDA library.

Construction of CEN and extraction of network modules. THD-Module
Extractor accepts a microarray gene expression dataset and two thresholds, the expression similarity threshold (δ) and the minimum neighborhood threshold (ρ), to construct CEN and to extract modules with high expression similarity using a similarity measure, SSSim. The method identifies border genes with low expression similarity and high semantic similarity from the highly co-expressed network modules. These border genes are analyzed to identify disease related genes from AD dataset. The threshold (δ) of the module extraction process is gradually updated from 0.9 to 0.5 by a factor α in each iteration. The minimum neighborhood threshold (ρ) is set to 3 in the conducted experiments. The generated co-expressed modules are validated in terms of p and q values.
For a list of genes, the p value represents the association of query genes with different GO terms. A lower value of p implies more biological significance of the network module. Similarly, q value is a statistical measure, which gives the minimum false discovery rate (FDR). In Tables 1 and 2, it is shown that the network modules obtained from Dataset 1 and Dataset 2 have less p and q values, respectively, which implies the biological significance of our results. Table 1 shows that the modules are biologically enriched with GO terms namely, cell cycle process, mitotic cell cycle process, chromosomal part, cell cycle, nuclear chromosome part with 3.801e-24, 5.760e-22, 2.351e-21, 2.188e-20, and 3.413e-18 p values, respectively. Table 2 shows that the modules extracted from Dataset 2 are enriched with functions like gluthathione peroxidate activity, response to wounding, monocarboxylic acid metabolic After preprocessing, it accepts two thresholds, namely, expression similarity threshold (δ) and minimum neighborhood threshold (ρ) to construct CEN and to extract modules with high expression similarity using similarity measure, SSSim. THD-Module Extractor updates δ by a factor α in each iteration of the module extraction process. The method identifies disease related genes by analyzing the border gene set extracted from the highly co-expressed modules and validates them using GO, KEGG Enrichment, pathway analysis, and existing literature. process, jasmonic acid biosynthetic process, perioxidase activity with q values 2.87e-16, 4.38e-15, 2.10e-12, 5.47e-10, 1.22e-9, respectively. Table 3 gives a comparison between FUMET 4 and the proposed method in terms of p values obtained for Dataset 3. Additionally, Tables 4 and 5 depict that GO terms obtained using the proposed method for Dataset 1 have lower p values and q values than that of FUMET and Qcut, respectively. This signifies the ability of the method to detect biologically enriched modules with low p and q values. As shown in Tables 3, 4

Application of THD-Module Extractor on AD dataset. Application of THD-Module Extractor on AD
dataset shows that the method can extract biologically and statistically enriched network modules. As shown in Table 6 and 7, the genes in network modules obtained from AD dataset have lower p and q values, respectively. For example, in case of Module 1 in Table 6 Analysis of modules to find genes related to AD. Each network module are extracted from a CEN to find those genes that have less expression similarity with the core gene, yet high semantic similarity with all genes. A gene or a node in CEN v i ∈ V is a core gene iff it is the starting node of a module m j and satisfies the following conditions.
1. v i ∉ any other module m k where j ≠ k.
where A is the Adjacency matrix, p and q denote the rows and columns in A, and n is the total number of genes.    On the other hand, a gene v b within a module is a border gene iff From each biologically enriched modules extracted from the CEN for AD dataset, the border genes are found to have low expression similarity with respective core genes and high semantic similarity among each other. These border genes are then analyzed to find the differentially expressed genes across normal as well as disease samples. The 2000 differentially expressed genes with highest variance across normal and disease samples are considered. These differentially expressed genes are then analyzed using queryMany() of mygene package in R library to find the entrez IDs of the corresponding genes. The semantic correlations among these entrez IDs are found using mgeneSim of GOSemSim package in R using Lin's semantic correlation measure 17 and Biological Process (BP) structure. In GO-BP structure, genes with annotations are involved with multiple functions 12 . Figure 2 shows an evidence of genes connected in p53 pathway. From the semantic correlation matrix, the genes for which semantic correlations are greater or equal to β are extracted, where β is the mean semantic correlation score. From this analysis, 912 genes are found to be interesting, with semantic correlation ≥ β.
These 912 interesting genes are analyzed using PANTHER and validated using GeneCards. Table 8 provides description of some of the interesting border genes, namely, APBB2, CASP2, CSNK1D, CDK5, HSD17B10, MAPT, PSEN2, and RCAN1. Some of the pathways associated with these border genes are depicted in Table 9 and their description is given in Table 10. A few of these pathways, namely, Wnt signaling pathway, Alzheimer disease-amyloid secretase pathway, Apoptosis signaling pathway, and Glycolysis are found to be related to AD, which are further validated using GeneCards and existing literature 18,19 .
Further, these 912 interesting genes are analyzed to find GO terms related to AD. The GO terms that are found to be related to AD are: membrane-bounded organelle, negative regulation of apoptotic process, and antigen processing and presentation of peptide antigen via MHC class, which are validated using GeneCards.
Moreover, a Web-based tool called TopoGSA 20 is used to study the topological structures. TopoGSA compares the interesting genes against KEGG database and finds KEGG enrichment scores for the given genes. Table 11 depicts that the interesting border genes extracted from the modules, obtained using THD-Module Extractor, are enriched with KEGG pathways. The KEGG pathways Wnt signaling, Apoptosis, p53 signaling pathway, Notch signaling pathway, Alzheimer's disease signaling, associated with the interesting border genes and validated using GeneCards and existing literature 18,19 , are found to be related to AD. A brief description of the pathways related to AD and the interesting border genes are given in Table 9.

Discussion
In the past two decades, a good number of methods, such as, FCM 21 Table 7.

p values and q values of the network modules for Dataset 4 (continued).
genes based on both expression and semantic similarity. The experiments conducted show that the border genes with low expression similarity and high semantic similarity are differentially expressed between normal and disease set of samples with high variance.
The key consideration is that genes or proteins function in co-ordination and not in isolation. The genes or gene products interact in biomolecular network or pathways, which are ultimately associated with the functions of living organisms and other various traits. Therefore, analysis of genes involved in molecular pathway is important for identifying disease-related genes or biomarkers, responsible for dysregulations of functions in disease patients.
AD is a common cause of dementia. In AD patients, it is found that there are some differences in brain tissues in the form of misfolded proteins, known as senile plagues and neurofibrillary tangles. The senile plagues or amyloid plagues are some clumps of β amyloid proteins, which block communication between neurons, and implicate neuronal death. The formation of β amyloid clumps at early onset of AD is initiated by mutations in proteins, known as Amyloid Beta (A4) Precursor Protein (APP) 27 , Presenilin 1 (PSEN1), and Presenilin 2 (PSEN2) (nlm.nih.gov). On the other hand, neurofibrillary tangles are formed due to mutations in a protein called Tau, which is responsible for ensuring clear passage for food molecules inside microtubule. Mutation in the protein Tau causes disintegration of microtubule and formation of tangles, which obstructs food passage to the cell neurons, resulting in death of neurons 28 . Involvement of genes or proteins in molecular pathways related to AD, such as, Wnt signaling, p53 signaling, Alzheimer disease-amyloid secretase, Apoptosis signaling, and Glycolysis, show a close relation with pathogenesis of AD.  In the proposed method, the interesting border genes are extracted from network modules using PANTHER and validated using GeneCards. The border genes are semantically similar as they share some common functionalities in the GO-BP structure, and therefore are found to be involved in molecular pathways. The pathways associated with the border genes related to AD are Wnt signaling, p53 signaling, Alzheimer disease-amyloid secretase, Apoptosis signaling, and Glycolysis. Some of the interesting border genes that are found to be related to AD and validated against the existing literature 1,29 are AATF, APBA2, APBB1, APOD, BACE2, CASP2, CAST, CDK5, CSNK1D, GAP43, HSD17B10, KCNIP3, KLK6, MAPT, NQO1, NRGN, OGT, PADI2, PSEN2, and RCAN1. Microtubule-Associated Protein Tau (MAPT) is linked with frontal lobe dementia 29 . PSEN2 is associated with deposition of longer form of amyloid-beta in AD brains, and Caspase 2 (CASP2) is involved in activation of cascade of caspases responsible for apoptosis, which increases cell death 1 . Gene CDK5 is associated with p53 signaling pathway. The p53 signaling pathway is a suppressor of cell growth, i.e., death of neurons in AD patients. Genes, such as, CSNK1D and PSEN2 related to Wnt signaling pathways are found to play an important role in neurodegenerative diseases 30 .

Method and datasets
The proposed THD-Module Extractor method is implemented using Matlab R2015a on a machine with an Intel(R) core I3-2120 CPU@3.30Ghz processor, a 2.00 GB RAM, and a 64-bit Windows 7 operating system. The algorithm is evaluated on four datasets. The description of the datasets is given in Table 12. Figure 1 depicts the work flow of the method. The various steps of execution are detailed below.   Preprocessing. Preprocessing of datasets is important before implementation of any scientific algorithm.
Many datasets contain noisy values and NULL values. The following preprocessing tasks are performed on the considered datasets using Matlab, before storing them in tab delimited form. The AD dataset GSE4226 consists of blood mononuclear cells collected from AD patients and age-gender matched normal individuals with 9600 probes and 28 samples. This work considers perturbation of gene expression values during progression of the disease, from normal to disease stage. Therefore, expression values of probes for AD patients and normal individuals are considered and rest are discarded.
After preprocessing of the datasets, the SSSim score of each gene pair is computed. However, computation of SSSim score for each gene pair takes lots of execution time. The datasets considered are large, often containing several thousands of genes, each with many expression values for multiple varying conditions. In our case, the preprocessed AD dataset contains expression values of more than 9000 genes under 28 different conditions, even after removing the unnecessary entries. This requires computation of more than 81 million getSSSim values. Traditional implementation of the measure on a standard CPU requires a huge amount of processing time. Therefore, this work leverages the parallel processing capabilities of GPU to speed up the computation of the large getSSSim correlation matrix.
NVIDIA CUDA implementation for SSSim computation. The NVIDIA CUDA 7.5 API and its C/C+ + based compiler (NVCC) are used to implement the getSSSim correlation matrix computation on GPU. An NVIDIA GeForce 980 GPU with compute capability 5.2 31 is used in this implementation. The computation of the getSSSim correlation matrix is distributed into a logical 2-D grid of thread blocks, where each thread processes one getSSSim value. The dimension of each thread block is 4 × 4 × 28, which specifies the number of threads per  block. The dimension of the grid is 2271 × 2271, which specifies the number of blocks in a grid. Figure 3 shows the distribution of threads into the grid. The GPU used in our implementation contains 2048 cores and can process equal number of threads simultaneously. In case of the AD dataset, the speedup obtained using the GPU for the correlation matrix computation is more than 1200 times, in comparison to a Matlab based implementation using a standard Intel-i7 CPU.
Once the score matrix has been computed, the CEN is constructed from the adjacency matrix, A. From the CEN, co-expressed modules are extracted for different values of δ and a static ρ. Each module contains a core gene that has maximum number of neighbor genes, satisfying the expression similarity threshold δ within the module. These co-expressed modules are then validated in terms of external validation metrics. A Web-based tool FuncAssociate 2.0 is used to find p values of each GO term associated with the genes of the co-expressed modules. Another Web-based tool called GeneMania 1 is used to find the q values of each GO term. From each co-expressed module extracted from AD dataset, the border genes that satisfy the expression similarity threshold δ with their neighbors, but not the minimum neighborhood threshold ρ, are found. THD-Module Extractor in finding interesting genes for AD. In addition to biologically relevant network module extraction, the proposed THD-Module Extractor allows identification of disease related genes by analyzing the border genes, extracted from the highly co-expressed modules. These border genes are analyzed using R library. Among the border genes, the k differentially expressed genes with highest variance between normal and disease samples are considered. Further, the entrez ID of each differentially expressed gene is found using mygene package of R. Thereafter, the semantic similarity of each of these border genes among each other is computed using GOSemSim package of R. On the basis of a static threshold, the border genes with high semantic similarity are filtered out. These interesting border genes are found to have low expression similarity with the core genes, yet high semantic similarity among each other. These interesting border genes are then validated using pathway analysis, KEGG enrichment analysis, GO terms enrichment analysis, and existing literature. The organization of the grid is hierarchical. The grid is a 2D array of blocks and each block is a 3D array of threads. For simplicity, the blocks are shown as 2D arrays of threads.