Extracting unrecognized gene relationships from the biomedical literature via matrix factorizations

Abstract

Background

The construction of literature-based networks of gene-gene interactions is one of the most important applications of text mining in bioinformatics. Extracting potential gene relationships from the biomedical literature may be helpful in building biological hypotheses that can be explored further experimentally. Recently, latent semantic indexing based on the singular value decomposition (LSI/SVD) has been applied to gene retrieval. However, the determination of the number of factors k used in the reduced rank matrix is still an open problem.

Results

In this paper, we introduce a way to incorporate a priori knowledge of gene relationships into LSI/SVD to determine the number of factors. We also explore the utility of the non-negative matrix factorization (NMF) to extract unrecognized gene relationships from the biomedical literature by taking advantage of known gene relationships. A gene retrieval method based on NMF (GR/NMF) showed comparable performance with LSI/SVD.

Conclusion

Using known gene relationships of a given gene, we can determine the number of factors used in the reduced rank matrix and retrieve unrecognized genes related with the given gene by LSI/SVD or GR/NMF.

Background

Latent semantic indexing based on the singular value decomposition (LSI/SVD) [1, 2] uses the truncated singular value decomposition as a low-rank approximation of a term-by-document matrix. Recently, LSI/SVD has been applied to gene clustering so as to retrieve genes directly and indirectly associated with the Reelin signaling pathway [3]. This approach may provide us with a powerful tool for the functional relationship analysis of discovery-based genomic experiments. However, this work did not utilize a priori knowledge of gene-gene relationships that are generally available. Moreover, the determination of the number of factors k used in the reduced rank matrix is still an open problem even though it is an important parameter that determines a concept space in which gene-documents are projected. In this paper, we suggest a method to estimate the reduced rank k in LSI/SVD by taking advantage of known gene relationships. In addition, we propose a gene retrieval method based on the non-negative matrix factorization (GR/NMF), which is a new framework for extracting unrecognized gene relationships from the biomedical literature.

Given a non-negative matrix A of size m × n and a desired reduced dimension k < min{m, n}, NMF finds two non-negative matrices W ∈ ℝ^m×kand H ∈ ℝ^k×nso that A ≈ WH. A solution to the NMF problem can be obtained by solving the following optimization problem:

(1)

where W ∈ ℝ^m×kis a basis matrix, H ∈ ℝ^k×nis a coefficient matrix, ||·|| _F is the Frobenius norm, and W, H ≥ 0 means that all elements of W and H are non-negative. Due to k <m, dimension reduction is achieved and the lower-dimensional representation is given by H. The NMF gives a direct interpretation due to non-subtractive combinations of non-negative basis vectors. In addition, some practical problems require non-negative basis vectors. For example, pixels in digital images, term frequencies in text mining, and chemical concentrations in bioinformatics are typically non-negative [4]. The NMF has been successfully applied to many problems including text data mining [5, 6] and gene expression data analysis [7, 8]. Non-negative dimension reduction is desirable for handling the massive quantity of high-dimensional data that require non-negative constraints. The determination of the reduced dimension k and the initialization of W and H are open problems. Some NMF algorithms [5, 9, 10] require that both W and H be initialized, while NMF based alternating non-negativity-constrained least squares that we describe in this paper only requires the initialization of H. We initialized a part of the matrix H by incorporating a known cluster structure and determined the reduced dimension k that can capture well the known gene relationships.

Results and discussion

For evaluation of our methods, we study two biological pathways: (1) Reelin signaling pathway and (2) Alzheimer's disease pathway. We try to extract unrecognized gene-gene relationships from the biomedical literature by taking advantage of known gene relationships. Table 1 shows 50 genes used for our experiments.

Table 1 The genes considered in the data set. The letters 'A', 'C' and 'D' in brackets show the relation with Alzheimer's disease, cancer, development, respectively.

Reelin signaling pathway

Reelin is a large extracellular protein that controls neuronal positioning, formation of laminated structures (including the cerebellum) and synapse structure in the developing central nervous system [11, 12]. Reelin binds directly to lipoprotein receptors, the very low-density lipoprotein receptor (VLDLR) and the apolipoprotein E receptor-2 (ApoER2), and induces tyrosine phosphorylation of the cytoplasmic adapter protein Disabled-1 (Dab1) by fyn tyrosine kinase. APOER2 is a gene alias name of LRP8. By using these knowledge, we chose five genes directly associated with Reelin signaling pathway, i.e. {RELN, DAB1, LRP8, VLDLR, FYN}.

We will examine if the following indirect gene relationships can be found by using the above knowledge. Dab1 is phosphorylated on serine residues by cyclin-dependent kinase 5 (Cdk5) [13]. The proteins encoded by CDK5R1 (p35) and CDK5R2 (p39) are neuron-specific activators of Cdk5. They associate with Cdk5 to form an active kinase. Apolipoprotein E (ApoE) is a small lipophilic plasma protein and a component of lipoproteins such as chylomicron remnants, very low density lipoprotein (VLDL), and high density lipoprotein (HDL). The ApoER2 is involved in cellular recognition and internalization of these lipoproteins. ApoE blocks the interaction of Reelin with its receptors. The Src related family member fyn tyrosine kinase mediates the effect of Reelin on Dab1 [14, 15]. MAPT encodes the microtubule-associated protein tau. Cdk5 is one of the major kinases that phosphorylates tau [16]. MAPT gene mutations have been associated with several neurodegenerative disorders such as Alzheimer's disease, Pick's disease, frontotemporal dementia, cortico-basal degeneration and progressive supranuclear palsy. The six genes indirectly associated with the Reelin signaling pathway are CDK5, CDK5R1, CDK5R2, APOE, SRC, and MAPT.

Alzheimer's disease pathway

We obtained the Alzheimer's disease pathway data from the KEGG pathway database [17]. From this pathway, we can obtain an overview of the general picture of the Alzheimer's disease pathway. Amyloid beta precursor protein (APP) encodes a cell surface receptor and transmembrane precursor protein that is cleaved by secretases to form a number of peptides. The pathway includes {APP, APBB1, LRP1, APOE, A2M, PSEN1, PSEN2, MAPT} among our 50 genes. These eight genes are known genes associated with the Alzheimer's disease pathway.

However, we cannot guarantee that the pathway contains all information regarding the Alzheimer's disease. We will determine whether we can find the following unrecognized knowledge from the above known knowledge. Amyloid beta precursor-like protein 1 (APLP1) affects the endocytosis of APP and makes more APP available for α-secretase cleavage [18]. Site-specific proteolysis of the amyloid-beta precursor protein (APP) by BACE 1 and γ-secretase, a central event in Alzheimer disease, releases a large secreted extracellular fragment (called APP(S)), peptides of 40–43 residues derived from extracellular and transmembrane sequences (Abeta), and a short intracellular fragment (APP intracellular domain) that may function as a transcriptional activator in a complex with the adaptor protein Fe65 and the nuclear protein Tip60. APP is closely related to APP-like protein (APLP) 1 and APLP2, and similar to APP, APLP1 and APLP2 are also cleaved by BACE 1 [19]. Amyloid beta precursor protein-binding, family A, member 1 (APBA1) stabilizes APP and inhibits production of proteolytic APP fragments including the Abeta peptide that is deposited in the brains of Alzheimer's disease patients. Some of the knowledge about genes is from the gene summary entries in the Entrez Gene database. The three unrecognized genes associated with the Alzheimer's disease pathway are APLP1, APLP2, and APBA1.

Performance comparison

For performance comparison, we tested two additional reference methods to identify unrecognized gene-gene relationships. The first method counts the number of shared PubMed citations cross-referenced in the Entrez Gene IDs for each gene. For the Reelin signaling pathway, we counted the number of PubMed co-citations between RELN and other genes. For the Alzheimer's disease pathway, we counted the number of PubMed co-citations between APP and other genes. If a paper is cross-referenced in two genes, the two genes are likely to have a direct or indirect association. The larger number of co-citations provides us with the more probable relationship between two genes, which may be a direct or indirect association. Provided we already knew genes directly associated with a pathway, we can find genes indirectly associated with the pathway. In Table 2, the number of PubMed co-citations between RELN and DAB1 was 9. Even though this method could not find some indirect relationships, i.e. (RELN – APOE) and (RELN – MAPT), it could find most of the direct and indirect relationships in the Reelin signaling pathway. However, it found only a known relationship (APP – PSEN1) in the Alzheimer's disease pathway (see Table 3). It cannot suggest potential gene relationships if they do not have co-citations.

Table 2 Genes directly and indirectly associated with the Reelin signal pathway. The cosine similarities between RELN and genes in the full space and the reduced dimensional space obtained from NMF are also presented. (n/a: not applicable)

Table 3 Known and unrecognized genes associated with the Alzheimer's disease pathway. The cosine similarities between APP and genes in the full space and the reduced dimensional space obtained from NMF are also presented. (n/a: not applicable)

The second method counts the frequency of gene symbol 'gene-B' in the gene-document of gene-A, and the frequency of gene symbol 'gene-A' in the gene-document of gene-B to find the relationship between gene-A and gene-B. For example, we searched for a symbol 'DAB1' in the RELN gene-document and a symbol 'RELN' in the DAB1 gene-document in order to find a relationship of (RELN – DAB1). The total frequency of symbol-match for (RELN – DAB1) was 47. Though this method could find all direct relationships, it could not find most of the indirect relationships except (RELN – CDK5) and (RELN – SRC) in the Reelin signaling pathway. This low recall problem is primarily due to inconsistencies in gene symbol usage in the literature. It could not find a known relationship (APP – APOE) in the Alzheimer's disease pathway. In contrast to these two reference methods, our gene retrieval method based on NMF could extract most of the direct and indirect gene-gene relationships by using the cosine similarity measure in the reduced dimensional space. We ranked genes by cosine similarity with a query gene RELN for the Reelin signaling pathway. A higher cosine similarity indicates a more probable relationship between two genes, which may be a direct or indirect association. In the full dimensional space, the ranks of two indirect relationships ((RELN – APOE) and (RELN – MAPT)) were 28 and 29, which were larger than those obtained from GR/NMF. APP was used as a query gene for the Alzheimer's disease pathway. In the full dimensional space, the ranks of two known relationships ((APP – APOE) and (APP – A2M)) were 13 and 18. In the reduced dimensional space obtained from NMF, the ranks were reduced so that we could capture the relationships. NMF (k = 3) can also be used to visualize genes in three dimensional space when it can retrieve most of the known relationships. Figure 1 shows gene relationships in the Reelin signaling pathway. Figure 2 illustrates gene relationships in the Alzheimer's disease pathway. By using different initializations of H, we were able to focus on the specific gene relationships in the different pathways. The proposed NMF initialization scheme was evaluated by retrieving genes associated with the Alzheimer's disease pathway in the reduced three-dimensional space obtained from NMF with k = 3. After computing the NMF with different initializations, we obtained F-measure values when 10 genes were retrieved. From 50 different random initializations, the NMF produced the maximal F-measure value (0.9524) only 36 times. Typically, the NMF is sensitive to random initialization since it converges only to a critical point.

Figure 1

Visualization of genes associated with the Reelin signaling pathway. Three-dimensional representations of 50 genes, which were obtained from NMF (k = 3) using an initial matrix H built from genes directly associated with Reelin signaling pathway. The j-th gene is located at (h_1j, h_2j, h_3j), where H ∈ ℝ^3×50 = [h _ij ]. (A red circle: RELN; A red circle and blue diamonds: genes associated with the Reelin signaling pathway; Black dots: other genes)

Figure 2

Visualization of genes associated with the Alzheimer's disease pathway. Three-dimensional representations of 50 genes, which were obtained from NMF (k = 3) using an initial matrix H built from known genes associated with the Alzheimer's disease pathway. The j-th gene is located at (h_1j, h_2j, h_3j), where H ∈ ℝ^3×50 = [h _ij ]. (A red circle: APP; A red circle and blue diamonds: genes associated with the Alzheimer's disease pathway; Black dots: other genes)

On the other hand, by using the proposed initialization scheme, the NMF achieved the maximal F-measure value 49 times. Using known biological knowledge, we were able to improve the probability that the NMF converges to a solution which reflects the known knowledge very well. In addition, since the convergence criterion Eq. (6) (below) is sometimes not tight enough for true convergence, guaranteed and faster convergence by the proposed NMF initialization scheme is required.

Tables 4 and 5 show the influence of the reduced dimension k on the LSI/SVD and GR/NMF retrieval performance. Recall, precision, and F-measures were computed when 10 genes were retrieved. Both cases showed that small k was enough to generate high F-measure values. GR/NMF showed comparable performance with LSI/SVD. By using the k-selection scheme in LSI/SVD, k = 3 was chosen for the Reelin signaling pathway and the Alzheimer's disease pathway. As for GR/NMF, we chose k = 3 for both pathways by using the k-selection scheme and the initialization strategy. Tables 4 and 5 show that LSI/SVD and GR/NMF exhibited excellent retrieval for indirect or unrecognized genes with k = 3.

Table 4 Influence of the reduced dimension k on gene retrieval of the Reelin signal pathway. Recall, precision, and F-measure were computed when 10 genes were retrieved.

Table 5 Influence of the reduced dimension k on gene retrieval of the Alzheimer's disease pathway. Recall, precision, and F-measure were computed when 10 genes were retrieved.

Practical applications

The LSI/SVD and proposed GR/NMF can elucidate unrecognized gene-gene interactions (i.e. edges in a gene interaction graph) from some known gene relationships. There are several types of identified gene-gene interactions. Firstly, a gene relationship identified by our methods can be a completely novel direct gene-gene interaction so that it needs to be confirmed by wet-laboratory biochemical experiments. Secondly, it can be an indirect gene-gene interaction implicit in a priori knowledge. For instance, if gene-A activates gene-B and gene-B inhibits gene-C, then gene-A and gene-C have an indirect gene-gene interaction. However, one need to be careful since there is still a possibility that gene-A and gene-C have a direct gene-gene interaction. Thirdly, it can be an explicitly known direct gene-gene interaction that is available in public databases although it was not recognized by users in advance.

Conclusion

In this paper, we have shown the utility of the SVD and the NMF extracting unrecognized gene-gene relationships from the biomedical literature. We have introduced a way to incorporate a priori knowledge into LSI/SVD and NMF in order to retrieve unrecognized documents related with a query document. Specifically, we have established the reduced dimension k estimation schemes for LSI/SVD and GR/NMF, which are generally applicable to information retrieval using the SVD and the NMF when there are some known relationships between a query document and other documents. The proposed GR/NMF takes advantage of a priori knowledge of cluster structure in its initialization step. It could retrieve unrecognized genes by using known genes associated with a biological pathway, which showed comparable performance with LSI/SVD. Extracting potential gene relationships from the biomedical literature may be helpful in building biological hypotheses that can be explored further experimentally.

Methods

Gene-document collection

In [3], a gene-document is generated by concatenating all titles and abstracts of the PubMed IDs cross-referenced in the human, mouse, and rat Entrez Gene IDs for each gene. However, this concatenation may cause skewed encoding in favor of a gene that has more PubMed IDs. In our study, to identify unrecognized gene relationships for n genes, a term-by-gene_document matrix A of size m × n is generated by the following scheme. A total of 50 genes were considered in three broad categories: (1) Alzheimer's disease; (2) cancer; (3) development (see Table 1). These 50 genes are the same genes as those used in [3]. For each Enrez Gene ID, we downloaded up to 10 of the most recent titles and abstracts, which were available as of July, 2006. Table 6 shows Entrez Gene IDs for human, mouse, and rat and the number of PubMed citations for each gene. As an intermediate step, we constructed a term-by-PudMed_document matrix A _p of size 8, 316 × 1, 273 in the form of MATLAB sparse arrays generated by Text to Matrix Generator (TMG) [20], where a PubMed-document is a document generated by concatenation of a title and an abstract for a PubMed ID. Then, from A _p , we built a term-by-gene_document matrix A of size 8, 316 × 50, where 50 gene-documents are centroid vectors for 50 genes.

Table 6 The number of PubMed citations associated with Entrez Gene IDs for each gene

We applied common filtering techniques (e.g. removal of common words, removal of words that are too short or too long, etc.) for the purpose of reducing the size of the term dictionary. Stemming was also applied. The m × n _p term-by-PubMed_document matrix A _p = [] was provided by using a log-entropy weighting scheme [2]. The elements of A _p are often assigned two-part values = l _ij * g _i , where l _ij is the local weight for the i-th term in the j-th PubMed-document, and g _i is the global weight for the i-th term. The local weight l _ij and the global weight g _i can be computed as

   l _ij = log₂(1 + f _ij ),

where f _ij is the frequency of the i-th term in the j-th PubMed-document, p _ij is the probability of the i-th term occurring in the j-th PubMed-document, and n _p is the number of PubMed-documents in the collection.

A gene-document vector a _i (the i-th column of A) can be easily compared with another gene-document vectors a _j (1 ≤ j ≤ n) in the full dimensional space. The similarity scores between two gene-documents (a _i and a _j ) can be computed as

(2)

Gene-document vectors having the higher cosine values are deemed more relevant to each other. In this gene retrieval method, a query gene vector is one of column vectors of A. This method tries to retrieve genes relevant to the given query gene. In order to compare gene retrieval methods quantitatively, we used the following performance measures. We defined the relevant genes, which include the query gene itself as well as genes related with the query gene. The recall and precision are defined as

The weighted harmonic mean of precision and recall, the traditional F-measure is defined as

Reduced rank estimation for gene retrieval via LSI/SVD

LSI is based on the assumption that there is some underlying latent semantic structure in the term-by-gene_document matrix that is corrupted by the wide variety of words used in gene-documents. This is referred to as the problem of polysemy and synonymy. The basic idea is that if two gene-documents represent the same topic, they will share many associating words, and they will have very close semantic structures after dimension reduction via the SVD. In LSI/SVD, if the matrix A has its SVD,

   A = U ΣV,

then its rank k approximation for some k < rank(A),

   A = U _k Σ _k

is considered, where the columns of U _k are the leading k left singular vectors, Σ _k is a k × k diagonal matrix with the k largest singular values in nonincreasing order along its diagonal, and the columns of V _k are the leading k right singular vectors. Then, S _k is the reduced dimensional representation of A, or equivalently, a gene-document vector a ∈ ℝ^m×1can be represented in the k-dimensional space as = a. Then, the similarity scores between two gene-documents ( _i and _j ) in the k-dimensional space can be computed as

Gene-documents having the higher cosine values in the reduced k-dimensional space are deemed more relevant to each other.

Here, we suggest a method to estimate the reduced rank k in LSI/SVD in order to retrieve unrecognized genes related with a query gene. If we can capture known gene-gene relationships in the reduced dimensional space obtained from LSI/SVD, we expect that the low-rank representations of gene-document vectors would be reliable to extract other gene relationships as well. This reduced rank k estimation scheme computes the following recall (), precision (), and -measure only from known genes relevant to the given query gene:

(3)

for various k values. It chooses the smallest k that shows the highest -measure value in order to retrieve unrecognized genes related to the given query gene.

Gene retrieval via NMF (GR/NMF)

In this section, we describe a gene retrieval method based on the NMF (GR/NMF) including the initialization of H and the reduced dimension k selection scheme.

NMF based on alternating non-negativity-constrained least squares (NMF/ANLS)

Given a non-negative matrix A ∈ ℝ^m×n, the NMF based on alternating non-negativity-constrained least squares (NMF/ANLS) starts with the initialization of H ∈ ℝ^k×nwith non-negative values. Then, it iterates the following ANLS until convergence:

(4)

which fixes H and solves the optimization with respect to W, and

(5)

which fixes W and solves the optimization with respect to H. Paatero and Tapper [21] originally proposed using a constrained alternating least squares algorithm to solve Eq. (1). Lin [22] discussed the convergence property of alternating non-negativity-constrained least squares and showed that any limit point of the sequence (W, H) generated by alternating non-negativity-constrained least squares is a stationary point of Eq. (1). After convergence, the columns of the basis matrix W are normalized to unit L₂-norm and the rows of H are adjusted so that the approximation error is not changed. Here, we adopt a fast algorithm for large scale non-negativity-constrained least squares (NLS) problems [23] to solve Eqs. (4–5). Bro and de Jong [24] made a substantial speed improvement to Lawson and Hanson's algorithm [25] for large scale NLS problems. Van Benthem and Keenan [23] devised an algorithm that further improves the performance of NLS. This algorithm deals with the following NLS optimization problem given B ∈ ℝ^m×kand A ∈ ℝ^m×n:

where G ∈ ℝ^k×nis a solution. It is based on the active/passive set method. More detailed explanations of this algorithm can be found in [23].

A method for initialization

Most NMF algorithms require initialization of both W and H, whereas NMF/ANLS described in this paper requires only initialization of H. In our approach, we incorporate a priori knowledge of gene relationships into the initialization of H. A gene-document is represented as a linear combination of basis vectors. For gene clustering by NMF, gene-documents that are dominated by the same basis vector belong to the same cluster. Here, we propose the following NMF initialization strategy. The elements of the first row of the initial matrix H ∈ ℝ^k×nare set to 1 only if the columns corresponding to a set of known genes _g are related with one another, otherwise set the elements to 0. For the other rows of H, the elements are set to 0 only if the columns correspond to _g , otherwise the elements are set to random numbers ∈ (0.25, 0.75). For instance, we know that RELN is related to DAB1, LRP8, VLDLR, and FYN. Thus, the elements of the first row of H have 1 only if the columns correspond to a set of genes _g = {RELN, DAB1, LRP8, VLDLR, FYN}, otherwise set the element to 0. The elements of the other rows of H have 0 only if the corresponding columns are related to _g , otherwise set the elements to random numbers ∈ (0.25, 0.75). Let us assume that the 4th, 5th, 6th, 7th, and 8th columns of H correspond to RELN, DAB1, LRP8, VLDLR, and FYN. Then, when k = 3, we can construct an initial matrix H as

where the values in the location of • are random numbers. The columns of the initial matrix H are normalized to unit L₂-norm. This initial matrix H contains a priori cluster structure.

Gene retrieval

NMF/ANLS is used to obtain the final W and H from the initial matrix H. Convergence is tested at every five iterations. The Frobenius norm of the error, i.e. f = ||A - WH|| _F , is computed at each convergence test. The convergence criterion is

(6)

where f _prev and f _curr are the Frobenius norms in the previous and current convergence tests respectively. The final matrix H ∈ ℝ^k×ncontains the low-rank representation of the term-by-gene_document matrix A. Hence, the similarity scores between two genes (i and j) in the k-dimensional space can be computed as

where h _j is the j-th column of the final matrix H. Genes having the higher cosine values in the k-dimensional space are deemed more relevant to each other.

NMF can generate different final H matrices because of the random numbers in the initial matrix H. Therefore, it is natural to repeat NMF with different initial matrices to obtain final H matrices. GR/NMF selects one of the final H matrices, which generates the highest -measure value using only known genes related with a query gene (see Eq. (3)). If there are several final H matrices that yield the same maximal -measure value, it chooses one producing the highest average of cosine values between the query gene and other known genes related with the query gene.

Reduced dimension estimation for GR/NMF

The determination of the reduced dimension k is also an open problem in the NMF. As with the k-selection scheme for LSI/SVD, we can estimate the reduced dimension k for GR/NMF by making use of known gene-gene relationships. The reduced k-dimensional representations of n gene-documents are obtained from the NMF. Then, the -measure is calculated only from known genes relevant to a query gene, after retrieving genes by cosine similarities in the reduced k-dimensional space. Even with a single k value, the NMF can generate different -measure values owing to the random numbers in the initial matrix H. Thus, to determine an -measure value for each k, the k-selection scheme selects the highest -measure value after computing -measure values with different initial matrices. This determines -measure values for various k values and then chooses the smallest reduced dimension k that produces the highest -measure value.