For every mouse gene, PubMed citations were obtained from the gene2pubmed repository available at NCBI. These citations are assigned either by professional staff at the National Library of Medicine or by the scientific research community via Gene Reference into Function (Gene RIF) portal. Since these citations are manually curated, we expect to have a very high precision for tagging correct citations to genes. We further filtered the non-specific citations by removing PMIDs that referred to more than 10 genes as these citations usually described high-throughput experiments mentioning a large number of genes with little to no significant functional information. After filtering, 21,022 mouse genes with at least one assigned citation remained in the collection. Among these, 1,111 genes were identified as TFs in the AnimalTFDB transcription factor database (Zhang et al., 2012). Out of the remaining 19,911 non TF genes, 7,695 genes were found to share at least one citation with at least one TF out of 994 TFs. A total of 45,229 gene–TF pairs with at least one shared citation were identified. For each such gene–TF pair, an abstract document was constructed by concatenation of titles and abstracts for each shared citation. A total of 26,781 unique citations were utilized in creating the abstract documents.

Text to Matrix Generator parser (Zeimpekis and Gallopoulos, 2006) was used to parse terms from the collection of gene–TF documents. All punctuation (excluding hyphens and underscores) and capitalization were ignored. In addition, articles and other common, non-distinguishing words were discarded using a stop list. Terms less than three characters in length were filtered out. A total of 106,895 terms remained after all the filtering. A 3-mode term × gene × TF sparse tensor was created where the entries of the tensor were frequencies of terms in abstracts shared between genes and TFs. Tensor construction from gene–TF documents is depicted in Figure S1 in Supplementary Material using a toy example with a small number of genes, TFs, and terms. The sparse tensor had 5,451,735 non-zero elements out of possible 817,621,682,850 (106,895 × 7,695 × 994) elements resulting in density of 6.66779 × 10⁻⁰⁶. In order to discount the effect of high frequency common terms in favor of more specific terms that might be better delineators between gene–TF combinations, each tensor entry f_ijk was scaled and transformed into l_ijk: (1)lijk=log2(1+fijk), where f_ijk is the frequency of the ith term in the document corresponding to the jth gene and the kth TF.

Given a 3-mode data tensor X of size m × n × p, with m mode-1 entities (terms), n mode-2 entities (genes), and p mode-3 entities (TFs), and a desired approximation rank k, the PARAFAC (Harshman, 1970) or CANDECOMP (Carroll and Chang, 1970) model approximates X as a sum of k rank-1 sub-tensors, each formed by the scaled outer product of a set of three vectors of lengths m, n, and p. The set of k vectors for each mode are usually grouped together in factor matrices A, B, and C of sizes m × k, n × k, and p × k, respectively. The columns of the factor matrices are normalized to unit length and the accumulated weight stored in a scaling vector λ. In addition, a constraint is imposed on the solution such that λ₁ ≥ λ₂ ≥ … ≥ λ_k. The tensor is expressed as: (2)X=∑i=1k λ i(ai∘bi∘ci), where a_i, b_i, and c_i represent the ith columns of the factor matrices A, B, and C, respectively; and ◦ denotes the outer product.

A common approach to fitting the PARAFAC model to data is an alternating least squares (ALS) algorithm (Tomasi and Bro, 2006), where one cycles over all factor matrices and performs a least-squares update for one factor matrix while holding all the others constant. We implemented a variant of the PARAFAC model called non-negative tensor factorization (NTF) that constrains the factor matrices to be non-negative. The goal of NTF is to find the best fitting non-negative matrices A∈ℝ+m×k, B∈ℝ+n×k, and C∈ℝ+p×k in the PARAFAC model that fit the data in X, corresponding to the following minimization problem: (3)minA,B,C∥X−∑i=1k ai∘bi∘ci∥F where F is the Frobenius norm. The norm of a tensor is similar to that of a matrix: (4)∥X∥2≡∑i,j,k (xijk)2.

A tensor can be matricized or flattened by rearranging the elements in a matrix. X^{(m × np)} represents 1-mode matricization of X, which is a matrix of size m × np where the index n runs the fastest over the columns and p the slowest. The matricized tensor can be expressed as: (5)Xm×np≈A(C⊙B)′, where ⊙ represents the Khatri–Rao product (Smilde et al., 2004) and ′ denotes the matrix transpose operator. The Kronecker product of matrices A and B is given by: (6)A⊗B=[a11Ba12B…a1nBa21Ba22B…a2nB⋮⋮⋱⋮am1Bam2B…amnB] and the Khatri–Rao product is defined as column wise Kronecker product: (7)A⊙B=[a11b1a12b2…a1nbna21b1a22b2…a2nbn⋮⋮⋱⋮am1b1am2b2…amnbn].

The 2-mode and 3-mode matricizations can be similarly expressed as: (8)Xn×mp≈B(C⊙A)′ (9)Xp×mn≈C(B⊙A)′.

To compute a k-factorization with NTF, we first initialized A, B, and C using the absolute values of k leading left singular vectors of the 1-mode, 2-mode, and 3-mode matricizations, respectively. This initialization was performed to make the subsequent NTF computation deterministic as well as to provide a potentially suboptimal starting point that may require less iterations to converge (Boutsidis and Gallopoulos, 2008). Each of the 3 matrices was treated as a non-negative matrix factorization (NMF) sub-problem: (10)minA∈ℝ+m×k∥X(m×np)−A(C⊙B)′∥F (11)minB∈ℝ+n×k∥X(n×mp)−B(C⊙A)′∥F (12)minC∈ℝ+p×k∥X(p×mn)−C(B⊙A)′∥F and solved in succession using the multiplicative update rule (Welling and Weber, 2001) modified to incorporate ∈ for stability. For example, to solve for A the update rule is: (13)Aiρ←Aiρ(X(m×np)Z)iρ(AZ′Z)iρ+∈, Z=(C⊙B).

The ∈ is a small number 10⁻⁹ added to the denominator in order to add stability to the calculation and guard against introducing a negative number from numerical underflow. The approximation rank k > 0 of NTF corresponds to the number of 3-way associations (sub-tensors) whose additive contributions approximate the information in the data tensor. Each triad {a_i, b_i, c_i}, for i = 1, …, k, defines scores for a set of terms, genes, and TFs for a particular 3-way association in the corpus. The scaling factor λ_i (after normalization) indicates the weight of the association for triad i. NTF was applied to the term × gene × TF log weighted frequency tensor using functions from MATAB tensor toolbox (Bader and Kolda, 2012). The procedure is described graphically in Figure 1.

A k-factorization delivers k sub-tensors. A sub-tensor i can be reconstructed via the outer products of the ith columns of the three non-negative factor matrices corresponding to terms, genes and TFs, respectively. For each such triad, the entities corresponding to the high magnitude scores in each column contribute more to the information content in the sub-tensor than the low magnitude ones and are, therefore, deemed more significant. We construed the triplet of sets of such significant genes, TFs, and terms as an annotated transcriptional module (ATM). Each ATM contained genes and TFs representing dominant elements of a putative transcriptional network, along with terms describing the functional interaction between these genes and TFs.

In order for the ATMs to be consistent with observed biological networks in terms of the number of entities, we established upper bounds for the number of genes, TFs, and terms in the ATMs. We observed the distribution of numbers of genes and TFs in 249 KEGG pathways and found a maximum of 184 genes and 21 TFs in a single pathway (after excluding the outliers). For the terms, we chose the upper bound of 300 as this is usually the maximum number of words allowed for publication abstracts and may be adequate to describe a transcriptional network.

For a given column and an upper bound n, we first obtained a truncated list D containing n highest scoring entities. We then calculated the number of significant (high magnitude) entities required to approximate the information content in the list as described in Alter et al. (2000). Briefly, contributions of each score g_i were computed as pi=giS, where S=∑j=1n gj (sum of all scores). Subsequently, we calculated the normalized entropy of the list as: (14)E=−∑i=1n pilog(pi)log(n), which is the fraction of information content in the list relative to a completely random list of the same size. The number of significant entities s was calculated as: (15)s=max(1,|pi>1n⋅E|).

Functional enrichment analysis was performed using two different human curated datasets. For the first set, we downloaded the list of 249 mouse related manually curated pathways and their associated genes present in the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2004). For the second set, we downloaded the list of 12,418 mouse related Gene Ontology (GO) (Ashburner et al., 2000) categories and their associated genes. The set of genes and TFs belonging to each ATM was evaluated for functional enrichment in the aforementioned KEGG pathways and GO categories using a hypergeometric test as described in Tavazoie et al. (1999). As a control, we compared the enrichment frequencies of ATMs to a set of 1,000 randomly generated gene–TF sets containing 8 genes and 2 TFs from the tensor. The numbers 8 and 2 depict the median number of genes and TFs across all 2,861 ATMs. To create each random gene–TF set, 8 genes out of 7,695 genes constituting the tensor were randomly picked. Similarly 2 TFs were randomly picked from 994 TFs in the tensor. This process was repeated 1,000 times.

The area under the curve (AUC) was used as a measure of quality of the terms associated with ATMs. The AUC will have the value of 1 for perfect ranking (all relevant terms at the top), 0.5 for randomly generated ranking, and 0 for the worst possible ranking (all relevant terms at the bottom) (Hanley and McNeal, 1982). The terms in the descriptions of KEGG and GO categories served as the gold standards.

Precision was used as a measure of accuracy of the gene–TF associations in an ATM. It was calculated as the ratio of the number of common entities (genes and TFs) between an ATM and the gold standard set, to the number of entities in the ATM. Mathematically, the formula was defined as: (16)|{genes and TFs in ATM}∩{genes and TFs in gold standard set}||{genes and TFs in ATM}|.

In addition, a set of ATM precision values was tested for significance by computing a right-tailed two-sample Welch’s t-test (Press, 1992) between the precisions for the ATMs being evaluated, and the precisions for 200 × n randomly generated gene–TF sets from the tensor. For each of n ATMs being evaluated, 200 gene–TF sets were randomly generated containing the same distribution of genes and TFs as the ATM.

Redundancy between any two sets of entities of the ATMs was computed as the Jaccard coefficient between the sets, which is defined as the ratio of the number of elements in the intersection and the number of elements in the union. The Jaccard coefficient will have a value of 0 for disjoint sets and a value of 1 for duplicate sets.

A k-factorization expresses the three-way gene regulatory interaction information in the tensor as additive superposition of latent information content in k sub-tensors. An approximation with k = 1 is expected to summarize the information content as the most dominant (or representative) transcriptional network gleaned from the literature, while k = 2 expresses the information content as superposition of two most dominant transcriptional networks. As k increases, the factorizations are expected to reveal more specific networks.

In order to quantitatively examine the information content across different k, we calculated a diversity coefficient d for each type of entity in the ATMs obtained at various k factorizations. For each entity type, we defined the diversity coefficient as the ratio of the number of unique entities in the union of all k ≥ 2 ATMs, and the total number of entities in the tensor. Figure 2 shows the diversity coefficients for genes, TFs, and terms. As expected, we found that the diversity of TFs, genes, and terms in the ATMs increased with increasing k. In addition, we computed the pairwise redundancy using the mean Jaccard coefficient between any two sets of entities in the ATMs for k ≥ 2 (Figure S3 in Supplementary Material). We found that while the overall diversity increases with increasing k, the pairwise redundancy between the entities remained constant, and that for any given k-factorization, the ATMs were disjoint. Taken together, these results indicate that at higher k, we were able to extract unique regulatory modules with more specific functional annotations.

To comprehensively evaluate the functional relevance of ATMs generated by our method, we examined if the genes and TFs in the ATM were significantly (p-value ≤ 0.05, hypergeometric test as described in Materials and Methods) enriched in KEGG or GO categories. We found that more than 94% of the 2,861 ATMs were enriched in at least one KEGG or GO category. This result indicates that tensor factorization across all ks produced biologically relevant ATMs. The median number of enriched categories in all ATMs at each k value was greater than chance (Figure 3). Also, the number of enriched KEGG and GO categories per ATM was higher at lower k. This result indicates that at lower k, ATM genes, and TFs have a broad range of functions, consistent with our earlier observation that the diversity and specificity of ATMs increased with higher k (Figure 2; Figure S3 in Supplementary Material). Consistent with these observations, we found that the diversity of enriched KEGG and GO categories in ATMs increased at higher k (Figure S4 in Supplementary Material), while the pairwise redundancy between categories remained low as k increased (Figure S5 in Supplementary Material). Taken together, these results suggest that the gene–TF ATMs are more functionally specific at higher k.

Interestingly, we found a higher diversity of KEGG pathways than GO categories at all k values (Figure S4 in Supplementary Material). This result suggests that KEGG pathways can be more specific than GO categories. Some KEGG and GO categories were very frequent among all ATMs. For instance, among the top ten overrepresented KEGG pathways, five were related to cancer, four related to infection, and one related to cytokine signaling (Table S1 in Supplementary Material). On the other hand, among the top ten overrepresented GO categories, four were related to development, three were related to regulation of gene expression, two related to cell proliferation, and one related to protein phosphorylation (Table S2 in Supplementary Material).

Finally, we evaluated the accuracy of the terms identified by NTF for each ATM by comparing them with KEGG and GO annotations. For each ATM enriched in at least one KEGG pathway or GO category, we created a gold standard using all of the terms in the titles and descriptions of the enriched pathways and categories. We compared the top ATM terms against these gold standards by examining the area under the curve (AUC) of receiver operating characteristics (ROC) curves. We found that more than 80% of the 2,861 ATMs produced AUCs above 0.6. There was no change in median AUC across k, albeit the range of AUCs increased with increasing k (Figure 4). The deviation of the median AUC from the lowest point increased faster with increasing k than the deviation from the highest point. This indicates that at higher k, some NTF-derived functional annotations are much less consistent with the descriptors of KEGG and GO categories.

A unique advantage of our approach is that it enables literature based discovery, i.e., new potential gene–TF interactions may be deduced from implied associations in the literature in the absence of direct (experimental) evidence. To evaluate the performance of our method in predicting potentially new gene–TF relationships, we used two different types of experimental data as gold standards. Chromatin immunoprecipitation sequencing (ChIP-Seq) is a genome wide technology that identifies direct binding sites for specific transcription factors on gene promoters. On the other hand, microarray expression analysis of tissues with a targeted deletion of specific transcription factor identifies downstream genes whose expression levels are either directly or indirectly altered by the transcription factor.

We used a set of previously published ChIP-Seq data for four transcription factors (Cebpα, E2f4, Foxa1, and Foxa2) in either liver or 3T3-L1 cells (MacIsaac et al., 2010) as gold standard validation sets. To evaluate our method, we focused only on the ATMs which contained the gold-standard TF. The number of ATMs for the four gold-standard TFs ranged from 11 to 54 and the number of genes and TFs in all the ATMs for the gold-standard TFs ranged from 34 to 112 and 8 to 45, respectively (Table 1). To calculate the average precision for each gold-standard TF, we compared the union of genes and TFs in all ATMs against the ChIP-seq validation set. The number of genes in the four validation sets ranged from 938 to 8,595 genes. As shown in Table 1, the average precision of our method ranged between 11% (Foxa1, liver) and 91% (E2f4, 3T3-L1). Except for Foxa2, all gene–TF associations extracted by our method were significantly (p ≤ 0.05, Welch’s t-test) higher than chance. The average precisions for a given TF in two different tissues varied, indicating that binding sites are tissue dependent. These results indicate that NTF is accurate for identifying TF target genes from the biomedical literature. Importantly, we found that 14–39% of the target genes predicted by our method were based on implied associations in the literature. An explicit association refers to the cases when a gene and a TF share an abstract, whereas an implicit association is inferred based on shared terms among genes or TFs in the same ATM. This result suggests that our method is useful for making literature based discoveries.

Our literature mining method identifies gene–TF associations based on functional information in the biomedical abstracts. Thus, it is possible that some of the gene predictions by our method are not the primary gene targets of the TFs, but rather genes in similar functional pathways whose expression levels are affected indirectly via a downstream TF. To test this possibility, we evaluated our method using four microarray datasets as gold standards. These datasets include differentially expressed genes in different mouse tissues (cerebellum, retina, or choroid plexus) as a consequence of targeted deletions in specific TFs (Atoh1 (Ha et al., 2015), Pax6 (Ha et al., 2015), and Otx2 (Omori et al., 2011; Johansson et al., 2013)). The number of genes and TFs in all the ATMs for the gold-standard TFs ranged from 57 to 152 and 34 to 62, respectively (Table 2). The number of genes in the validation sets ranged from 2,137 (Atoh1, cerebellum) to 11,689 (Otx2, choroid plexus). The average precision for all four validation sets ranged from 0.33 to 0.42 and were all significantly (p ≤ 0.05, Welch’s t-test) higher than chance. Importantly, approximately 27–63% of the predictions from our method were based on implied associations extracted from the literature. Since transcriptional regulation by TFs are tissue specific, we narrowed the ATMs to those which explicitly mentioned the experimental tissue. We calculated average precision for ATMs that contained either one or two keywords associated with the tissue. In all but one case (Otx2, choroid), the average precision values improved when we considered only the tissue-relevant ATMs.

Next, we focused on one dataset (Atoh1) to carefully examine the characteristics of each ATM. As indicated above, the average precision across all 38 Atoh1 ATMs was 0.33. However, the precision of individual Atoh1 ATMs ranged between 0 and 0.85 (Figure 5). In general, higher k produced higher precision, indicating that higher classification specificity is correlated with precision. Interestingly, precision was not correlated with whether Atoh1 was ranked first in the ATM. The average precision was 0.23 for the 17 ATMs in which Atoh1 was top ranked. In contrast, Atoh1 was first ranked in only one of the top three ATMs with the highest precision (Figure 5). Finally, six ATMs had a precision of 0. Upon examination, we found that these ATMs were significantly enirched (p ≤ 0.05) for “auditory receptor cell differentiation” and other GO categories related to ear development. Consistent with these GO classifications, the terms in the ATMs were related to hair cell differentiation and development. Thus, it is likely that the precision of these ATMs is low because the gold standard was derived from microarray experiments using cerebellar tissue rather than auditory tissue.

Finally, as a benchmark, we compared the predictions of our method with those obtained by RegNetwork (Liu et al., 2015) for the same experimental ChIP-Seq and TF knockout microarray datasets described above as gold standards. RegNetwork is a database containing target predictions for TFs, which have been sourced from more than 20 interaction databases. Tables S3 and S4 in Supplementary Material show the precisions obtained by NTF along with those obtained by RegNetwork predictions. With the ChIP-Seq gold standards, we found that the average precision of NTF method was 0.48 compared to 0.38 for RegNetwork. NTF outperformed in 4 out of 6 ChIP-Seq datasets. For the two other datasets, RegNetwork outperforms but only marginally. In contrast, with the TF knockout microarray gold standards, the average precision of NTF (0.37) was better than that of RegNetwork (0.16). NTF precision was higher for 3 out of the 4 microarray datasets. Importantly, RegNetwork predicted targets for Foxa1 and Otx2 did not match any targets in the gold-standard datasets, whereas NTF identified target genes with 0.11 and 0.42 precision, respectively (Tables S3 and S4 in Supplementary Material).

Current methods for functional interpretation of high-throughput genomic data rely on manually curated knowledge bases, such as GO, KEGG, and a variety of other resources. It is generally accepted that the rate of manual curation is not sufficient for the amount of genomic data that is being generated in various species (Baumgartner et al., 2007). Moreover, there is a quantifiable bias in manually curated knowledge bases toward more popular genes and a substantial drift in the annotation of genes over time (Gillis and Pavlidis, 2013). Previous work have focused on using text-mining approaches to enhance curation of GO and other knowledge bases (Chagoyen et al., 2006; Couto et al., 2006; Thomas et al., 2015; Peng et al., 2016).

Here, we examine the performance of our NTF approach in predicting (suggesting) GO classifications. GO predictions were based on two approaches: (1) Guilt-by-Association (GBA) or (2) term mapping. In the GBA approach, new genes or TFs were assigned to a GO category which is significantly enriched for one or more ATMs. These new genes or TFs are part of the ATMs but have not been explicitly assigned to the enriched GO category by the curators. In the term mapping approach, novel and more specific GO categories for a given ATM were predicted based on the number of terms from the ATM that overlapped with GO category descriptions. These new GO categories were not significantly enriched for the ATM under consideration. The results for both approaches were evaluated by manual examination of the biomedical literature.

In order to evaluate the GBA approach, for the 38 ATMs containing Atoh1, we randomly chose 3 GO categories that were found to be significantly enriched in some of those ATMs (Table 3). For example, for category “Axon guidance” that was significantly enriched for 24 ATMs containing a total of 104 genes and TFs, only 16 genes and/or TFs were explicitly assigned to the category by the GO curators. Upon manual examination of sentences in Medline abstracts, we found that out of the 104 genes and TFs, 29 were functionally related to axon guidance. Among the 29, 11 were already assigned to the GO category but the rest (18) were not. The latter are candidates for assignment to the “axon guidance” category, potentially increasing the assignment of new genes and TFs to the category by 1.6-fold. Five of the existing GO assignments could not be validated through manual evaluation. Overall, for three randomly selected significantly enriched GO categories (Table 3), the GBA method increased the assignment of new genes by an average of 1.2-fold. The average precision for the GO assignment by this method was 0.27 (0.14–0.38 range). Notably, on average, 52% (26–100% range) of the existing GO annotations were not validated by our manual analysis, indicating that there is substantial error in GO curation (Table 3).

In order to evaluate the term mapping approach, for the 38 ATMs containing Atoh1, we predicted 7 novel GO categories whose names and descriptions were found to have high overlap with the terms from some of those ATMs (Table 4) but were not significantly enriched for any ATM. For example, for category “Spinal cord oligodendrocyte cell differentiation” whose description had high overlap with the terms of 3 ATMs containing a total of 44 genes and TFs, only 1 gene or TF was explicitly assigned to the category by the GO curators. This low number contributes to the reason the category was not significantly enriched for any of the 3 ATMs. Upon manual examination of sentences in Medline abstracts, we found that out of the 44 genes and TFs, 9 were functionally related to spinal cord oligodendrocyte cell differentiation. Among these, 1 was already assigned to the GO category but the rest (8) were not. The latter genes and TFs are candidates for assignment to the category, potentially increasing the assignment of new genes and TFs to the category by 8-fold, which make the category a plausible candidate for being significantly enriched for the ATMs whose terms have high overlap with the category description. On average, 3.9-fold (1- to 8-fold range) more genes were assigned to novel predicted GO categories, but with an average precision of 0.15 (0.06–0.23 range) (Table 4). Most of these categories were more specific than the currently enriched categories.

In general, the terms associated with ATMs provide a greater level of functional specificity than structured categories in ontologies. For example, Figure 6C shows the genes, TFs, and terms associated with one ATM in which Atoh1 is top ranked. As a comparison, the enriched GO and KEGG categories are also displayed. Whereas the GO categories indicate that this group of genes and TFs are associated with inner ear development and auditory receptor cell differentiation, the terms (shown in italics below) in the ATM suggest that they are involved in differentiation of the sensory hair cells in the cochlea located in the organ of Corti. Atoh1 is a basic helix loop helix (bhlh) domain transcription factor involved in early development. Targeted deletion of Atoh1 gene in mice results in loss of hair cells in the organ of Corti (Chonko et al., 2013). Several gene symbols also appear in the top ranked terms such as Sox2, Notch, Jag1, Prox1, and Math1. Interestingly, Math1 is the alias for Atoh1. Direct evidence suggests that Notch1 and Jag1 signaling pathway is required for activation of Sox2 and Atoh1 expression (Neves et al., 2011). Activation of Atoh1 by Sox2 transcription factor is required for hair cell development in cochlea with respect to both expansion of the progenitor cells in the cochlear epithelium and initation of hair cell differentiation (Kiernan et al., 2005; Kempfle et al., 2016). Conversely, Prox1 transcription factor directly suppresses Atoh1 expression. Sox2 is expressed in type 2 vestibular hair cells and in supporting cochlear and vestibular epithelium (Hume et al., 2007).

In order to facilitate manual examination of individual ATMs, e.g., the aforementioned Atoh1 ATMs, the publicly available web tool called Transcriptional Regulatory Modules Extracted from Literature (TREMEL) was developed. It is accessible at http://binf1.memphis.edu/tremel. TREMEL provides a searchable interface for the genes, TFs, and terms for all 2,861 ATMs (Figure 6).

The user can query the tool with either genes, TFs or terms, or a combination of any of the 3 types of entities in additional query boxes. The entity type can be selected from the drop down list to the right of each search box. Clicking the “+” button to the left of the first search box opens an additional search box. A maximum of 3 search boxes are allowed. The gene and TF queries need to be the official symbols designated by the National Center for Biotechnology Information (NCBI). Only one symbol is allowed per search box. The term query can be any single keyword such as “cancer,” “neuron,” “transcription,” etc. All queries are case insensitive. The output of the tool consists of two panels.

The top panel is comprised of a 3-dimensional interactive plot that shows all ATMs containing the search box entities, as points. The axes of the plot correspond to the NTF approximation rank k, the ATM#, and the rank of the queried entity (first search box only) in the ATMs. An ATM can be selected in the panel by clicking on its corresponding point in the 3-D plot. The color of the selected point changes to red and the colors of the remaining points corresponding to all other ATMs are depicted in terms of similarity to the selected ATM. The most similar ATMs are colored in shades of red while the least similar ones are colored in shades of blue. The similarity between any two ATMs is calculated as the Jaccard coefficient between the sets of genes and TFs in the respective ATMs. Upon initial search completion, one ATM is preselected and colored in red. This initial selection is performed in a manner such that the ATM with the lowest queried entity rank is picked. Ties are resolved in favor of the ATM with the lowest NTF approximation rank.

The bottom panel contains several sub-panels, each corresponding to an ATM point in the first panel. The top sub-panel corresponds to the selected ATM, and the remaining ATMs are ordered according to their similarity to the selected ATM. Each sub-panel displays the ranked genes, TFs and terms of the corresponding ATM, as well as the enriched GO categories and KEGG pathways. Clicking a sub-panel expands it to display its contents, and closes the previously open sub-panel. The contents of only one sub-panel are viewable at a time.