Finding novel relationships with integrated gene-gene association network analysis of Synechocystis sp. PCC 6803 using species-independent text-mining

Networks constructed using microarray and yeast-2-hybrid data

To create a Synechocystis 6803 co-expression network, 68 data sets from a large scale transcriptomics study (Singh et al., 2010) were used. The transcriptome data was collected and stored as fold change (log2 (treatment/control)) of gene expression values in tab-delimited text files. The data was thereafter subjected to further analyzing after importation into the analyzing and visualizing platforms Cytoscape 2.8.2, 3.0.1 and 3.3.1 (Shannon et al., 2003; Smoot et al., 2011), depending on available plugins. The ExpressionCorrelation plugin (Hui et al., 2008) was employed to generate a co-expression network using the expression values. A similarity matrix was calculated using the Pearson correlation coefficient with a strength threshold of ±0.7. In order to identify an appropriate threshold, the ExpressionCorrelation plug-in histogram function was used to select a reasonably sized network. The obtained co-expression based gene network (1,886 nodes and 10,187 edges) is referred to as CoEx. A second yeast two-hybrid (abbreviated Y2H) protein-protein interaction network was constructed by importing into Cytoscape a list of identified protein-protein interactions from an available data set (Fields & Song, 1989; Sato et al., 2007).

Text-mining data

BioContext (Gerner et al., 2012) and EVEX (Van Landeghem et al., 2013b) are one of only few resources providing gene/protein association from the entire literature databases that have been made publicly available. Both EVEX and BioContext apply their text-mining tools to PubMed abstracts and PubMed Central (PMC) full-text articles and are able to extract a similar number of events (36 vs. 40 million for BioContext and EVEX, respectively).

Besides having more extracted events, EVEX also aggregates events of genes in the same family providing hypothetical network where gene/protein association can be conveniently inferred from closely related sequences. Though both resources are considered equivalent extracting the same information from the same data set, we consider EVEX database a more suitable text-mining resource for our studies, due to its data on gene-family association and more recently updated articles.

The network from the EVEX database is composed of two data sets following the different releases of EVEX namely, EVEX-2011 and EVEX-2013. EVEX-2011 (Van Landeghem et al., 2011) is the first public release of the EVEX text-mining database which covers the literature up until June 2011 (http://www.evexdb.org/ (Van Landeghem et al., 2013a; Van Landeghem et al., 2013b; Van Landeghem et al., 2011)). EVEX-2013 (Van Landeghem et al., 2013b) was released with the extended coverage of articles from June 2011 up to June 2012 and an updated gene family assignment. Both of the EVEX data sets (EVEX-2011 and EVEX-2013) were combined and used in the present study.

EVEX data was generated using natural language processing tools primarily based on machine learning (ML) to automatically extract cellular processes and interactions among genes and their products such as RNAs and proteins (genes for short). The tools perform three significant steps namely “name entity recognition”, “event extraction” and “name entity normalization”, which will be discussed here briefly. Firstly, the tools perform name entity recognition by identifying the gene mentions in the documents. The systems then extract the biological events for each gene mention by identifying words or phrases discussing cellular process such as regulation and phosphorylation and link them to corresponding genes. Finally, to be able to link the genes to information in other databases, genes are normalized by mapping to the Entrez Gene database and respective family identifiers. In case of organism ambiguity, i.e., when the organism is not explicitly stated for a particular mention thus preventing it from being normalized to a single unique identifier, the mention is only mapped to a gene family. Full details of the EVEX text-mining pipeline generating has been described previously (Van Landeghem et al., 2013b).

The text-mining network was constructed by retrieving Synechocystis 6803 genes (nodes) and their associations (edges) from the EVEX database. This was extended by the addition of associations of all Synechocystis 6803 gene homologs originating from the Ensembl resource and the GENIA corpus (Kim, Ohta & Tsujii, 2008).

The EVEX network was further enriched with additional information obtained from the EVEX database. The edge attributes included all the organisms for which the relationship was reported, the taxonomic distance between each organism and Synechocystis 6803, fine-grained details of the relationship such as the types of the regulation (positive, negative and unspecified), speculation, negation and text-mining prediction confidence score. The node attributes also include Synechocystis 6803 gene descriptions, symbols, synonyms, and Entrez Gene identifiers.

Finally, in the NCBI Entrez Gene record, the functions of a well-characterized gene are described by human annotators based on experimental evidence. While generally very useful, oftentimes these descriptions lack specificity, e.g., for genes annotated as “hypothetical”. Further, new sequences with no supporting evidence naturally lack this annotation altogether.

To address this deficiency, we also provide descriptions based on the single most prevalent function among the genes in each family. For a small gene family (i.e., <5 genes), the diverse descriptions can be manually combined and selected to represent the common functions of genes in a given family. However, this process is not suitable for a large family with thousands of genes. Here we used the method called “canonicalization” described in (Van Landeghem et al., 2011) to select the representative description of the family. First, we collected the descriptions of all genes in a family from NCBI Entrez Gene records. We then reduced the orthographic differences by lowering the case and removing all non-alphanumeric characters such as empty space, parentheses and apostrophes. The description of the gene family is the most common canonical form of descriptions shared by most genes in the family.

The three networks, CoEX, Y2H and EVEX, were thereafter integrated using the Cytoscape tool “Advanced network merge”. The merge was carried out based on the Entrez Gene identifiers. For those data sets that did not contain such node identifiers, these were obtained by mapping through CyanoBase gene identifiers. The resulting merged network (IntNet) is provided as File S1 and the attribute annotations are listed in File S2. Although the constructed networks inherently contain complex structures (Kashtan et al., 2004), the present work has only focused on simpler triad motifs.

Case study 1—novel candidates with a potential relationship to SigE

SigE (sll1689) is a sigma factor that has been demonstrated to influence central carbon metabolism with broad impact, as evidenced by a shift in the distribution of central carbon metabolites in response to the deletion of sigE or over-expression of SigE (Kloft, Rasch & Forchhammer, 2005; Sundaram et al., 1998). The first neighbour GBA and script-based clusters are shown in Fig. 4A, including several interesting candidates. Firstly, we noted a link between SigE and slr1055 (ChlH), a light- and Mg²⁺-dependent anti-sigma factor shown previously to have specificity for SigE (Osanai et al., 2009). However, this link was not based on the article that demonstrated this relationship in the first place (Osanai et al., 2009). Instead, SigE connects with ChlH through edges of all three network types, a direct Y2H edge, the lead to identifying the role of slr1055 in the first place (Osanai et al., 2009), and indirect edges via sll0306 (SigB, EVEX) and sll1886 (hypothetical protein, CoEx). The experimentally confirmed relationship between ChlH and SigE therefore verifies the conclusion of the relationship that can be drawn from the present network even in the absence of the direct text-mining link.

Several known proteins with an established role in nitrogen-metabolism (e.g., NtcA, PII (Kloft, Rasch & Forchhammer, 2005)), or the circadian clock (KaiB (Hitomi et al., 2005)) were also found to be connected to SigE, in addition to others without any meaningful annotation. The automated script (Fig. 4B) suggested a central role for sll1886 (annotated as hypothetical protein) with a close connection to SigE. Sll1886 harbors a putative zinc binding domain and shows weak homology to di-haem cytochrome C (Vandenberghe et al., 1998), suggesting the possible involvement of electron-transfer. Interestingly, a manganese transport component (MntB, sll1600) was also part of the script-based cluster which is relevant given that ChlH is Mg²⁺-dependent.

In comparison, we also searched for candidate genes to SigE using STRING-db (Franceschini et al., 2013) with sll1689 as input (Fig. 4C). This produced a network of 11 nodes at the default setting. When the script- and STRING-db based networks were compared, the intersection between the two networks contained only three genes; sll1689, sll1423 (ntcA) and sll0687 (sigI). Interestingly, whilst the STRING network contained an association with glnA, the script-based network contained an association with glnB - both genes have important roles in nitrogen metabolism (Herrero, Muro-Pastor & Flores, 2001). Overall, many of the nodes in the STRING network (Fig. 4C) are related to gene transcription (RNA polymerase related gene products), whilst the script-network (Fig. 4B) is dominated by genes with a known role in nitrogen metabolism, as has also been confirmed experimentally (Muro-Pastor, Herrero & Flores, 2001). The former network has no ‘unknown’ members, whilst at least one completely unknown, yet intricately connected, member (sll1886) is present in the latter network. Notably, sll1886 is co-located on the genome to a “two-component sensor histidine kinase” (sll1888) which also is a member of the same CoEx network as sll1886 and ntcA (sll1423) (Figs. 4A, 4B). This strengthens the argument that sll1886 may play an important role in nitrogen metabolism.

Case study 2—NADP(H)-metabolism

The role of the pentose phosphate pathway (PPP) in cyanobacteria under daylight conditions is not entirely clear given that NADP⁺ is a major electron acceptor of electrons generated by water-splitting photosynthesis. A part of the metabolic flux through the carbon fixing CBB cycle has been measured to pass through the oxidative branch of PPP (oxPPP) under daylight conditions (Young et al., 2011) though the optimal solution for biomass flux in stoichiometric models did not incorporate any oxPPP flux (Knoop et al., 2013). We were curious about the metabolic role that key-enzymes responsible for NADP⁺-reduction in fermentative microorganisms may have in an autotrophic system and how they are regulated. The objective in the following analysis was therefore to use the network analysis in order to identify novel candidate genes.

A first neighbour GBA of IntNet with all pre-defined six NADPH key genes generated a complex network of 72 nodes and 194 edges (Fig. 5A) (File S6), including OpcA, the unique cyanobacterial Zwf activator (Hagen & Meeks, 2001). In contrast, only two of the 6 key genes listed for NADP(H)-metabolism were retained by the script (Fig. 5B, 18 nodes and 50 edges): Zwf (slr1834, catalyzing the first committed step of metabolic flux into PPP) and Icd (slr1289), catalyzing the only NADP⁺-reducing step of the TCA-“cycle”.

Looking closer at the script-based network, Zwf forms a motif with slr0952 (annotated as fructose-1,6-bisphosphatase (FBPase)) and sll0508 (annotated ‘unknown protein’) via three different data-types (Fig. 5C). Sll0508 has low similarity to other proteins and there are no hits from a search with the SIB Motif Scan (incl. Pfam, PROSITE, HAMAP etc.). This slr0952-containing motif is interesting as it suggests a link between oxPPP and gluconeogenesis. In other cyanobacteria, multiple FBPases have been identified and some of the encoding genes are co-located with zwf (Summers et al., 1995).

Another interesting candidate gene, found only in the CoEx network, is Slr1194. This node is annotated as a ’1 protein’ that exhibits a high percentage similarity to a ’Mo-dependent nitrogenase family’ protein in Cyanothece sp. PCC 7424, and links to Zwf via slr1793 (talB) and slr1734. The latter gene is a homolog of OpcA, an allosteric regulator and activation factor of Zwf in other cyanobacteria (Hagen & Meeks, 2001).

Zwf also forms several motifs with rpaB (slr0947) that also include the PPP genes gnd (sll0329) and talB (slr1793) (Fig. 5B). RpaB is a regulator involved in controlling energy transfer between phycobilisomes and PSII or PSI. The relationship between RpaB and genes encoding enzymes in PPP suggests the possibility that also PPP flux may be controlled at least in part by RpaB in response to light quality and/or quantity, or another signal reflecting the internal redox-status.

Case study 3—Probing the role of an incompletely known gene or gene set—PntAB

Synechocystis 6803 harbors two genes (slr1239 (pntA) and slr1434 (pntB)) encoding a putative dimeric NADPH:NADH-transhydrogenase. PntAB has been shown to catalyze the proton gradient dependent transfer of electrons from NADH to NADP(H) in E. coli (Sauer et al., 2004). In Synechocystis 6803, we would expect under optimal photosynthetic conditions that NADP⁺ is efficiently reduced by PetH, the Ferredoxin:NADP-oxidoreductase linked to PSI. PntAB may therefore only be important for the supply of NADP(H) under conditions of limiting light (e.g., during the night) and/or in order to re-oxidize NADH formed by NAD(H)-dependent reactions (Kämäräinen et al., 2017). Hence, although PntAB is well-known in fermentative microorganisms it remains unclear what role it may have in cyanobacteria, thereby falling into the category of incompletely known genes.

No motifs satisfying the criteria of the script-based filter were found including either PntA (slr1239) or PntB (slr1434). Nevertheless, a GBA-cluster was extracted using both genes as key genes (Fig. 5D). Both slr1239 and slr1434 form a co-expression based cluster with an operon (slr0144-slr0152) called Pap (Photosystem II assembly proteins) (Wegener et al., 2008) and the essential ferredoxin PetF (slr0150; Fig. 5E). The connection is quite convincing as PntA shows CoEx edges with slr0144 whilst both PntB and PetF share CoEX edges with several of the other genes in the operon, though not slr0144. The presence or absence of the Pap operon does not influence growth under so far tested conditions, although deletion mutants display a reduced capacity to evolve di-oxygen (Wegener et al., 2008). Why would there be a connection between the Pap operon and PntAB? PntAB has the role in fermentative microorganisms of catalyzing electron-transfer between one major electron acceptor-donor and another, though not ferredoxin. Several genes of the Pap operon are predicted to contain Fe-S clusters, co-factors that typically also are involved in electron transfer, the only common theme so far; this connection deserves further experimental attention to resolve.

Case study 4—iron sulphur cluster metabolism

As mentioned above, iron-sulphur (Fe-S) clusters are inorganic protein co-factors that are typically involved in electron transfer. They are assembled in cyanobacteria using the SUF system, even though genes with homology to members of the ISC system (the dominant system in E. coli) also are present in the Synechocystis 6803 genome (Balasubramanian et al., 2006). It has been established that SufR (sll0088) is an Fe-S containing negative transcriptional regulator of the remaining SUF members (sufA, sufB, sufC, sufD, sufS) (Wang et al., 2004). Interestingly, a first neighbour GBA with all of the above key genes (Fig. 6A, File S6) resulted in a single cluster with two divided parts, an upper part containing all the catalytic SUF members, and a second lower part containing SufR. Even though SufR is clearly the transcriptional regulator of the other SUF members, there is surprisingly no direct connection between SufR and the other SUF members. Instead, SufR forms an intense CoEx cluster with a series of genes annotated mainly as ‘hypothetical’. Three of these are iron-related proteins: PerR (slr1738), sll1202 (homolog to iron transporters) and BfrA (sll1341; bacterioferritin homolog). In contrast, the upper SUF operon cluster contains four genes encoding predicted Fe-S containing proteins: The PSI subunit psaA (slr1834), bioB (slr1364), sll0031 (hypothetical) and spoT (slr1325). A possible reason for the lack of a direct association between SufR and the remaining SUF operon may be that SufR is not the only regulatory factor controlling SUF expression, or that its control is conditional.

Filtration of IntNet using all Fe-S key genes generated two smaller clusters (Fig. 6B). Whilst no obvious insight was obtained from the SufR-containing motif, the second cluster contained three SUF operon members connected both by EVEX and CoEx. Interestingly, all text-mining edges originated from a diverse collection of bacteria that did not include any cyanobacteria.

Case study 5—alkane biosynthesis

The two genes encoding the catalytic enzymes of the alkane biosynthesis pathway (Schirmer et al., 2010), and which is present in most but not all cyanobacteria, forms an extended apparent operon in most species where it is found (Klähn et al., 2014). Since the alkane biosynthesis reaction so far does not work as efficiently as needed for economically sustainable fuel production (Eser et al., 2011; Kallio et al., 2014), we were curious whether missing elements required for effective catalysis could be represented in this apparent operon. In Synechocystis 6803, however, only three of the apparent operon members are co-located on the genome, sll0207-sll0209. For the assembly of key genes, we therefore included homologs in Synechocystis 6803 to the most commonly observed members of the alkane biosynthesis operon in cyanobacteria in general (Table 1), even if they are not co-located on the genome in Synechocystis 6803. In this analysis (Fig. 7, File S6), however, most of the operon members did not form a joint cluster with the exception of slr0426. A possible contributing reason for this outcome is that the biosynthetic system is unique to cyanobacteria (Schirmer et al., 2010) and that it has not yet been studied much. Consequently, it is not well-represented in the EVEX network.

Case study 6—Screening for the role of genes annotated as ‘hypothetical’ or ‘unknown’

We considered the possibility to utilize the script in order to obtain an insight into the possible role of all genes that are annotated as ‘hypothetical’ or ‘unknown’. The rationale was that the local context of genes without an annotation may provide insight into its possible role and that the script would allow the most important local context to be identified. All genes without an annotation were therefore employed, one at a time, as an entry gene for the automated script. The criteria of this script demanded as previously that more than one relationship type was present, plus the additional new demand that at least one of the members of the local context had an existing annotation. Over 5% of hypothetical/ unknown genes (112/1913) satisfied these criteria. The combined network with filtered motifs was composed of 331 nodes (Fig. 8A; File S7). Around 60% of these patterns were derived from a combination of CoEx/Y2H and around 40% from EVEX/(CoEx/Y2H). These 112 putative genes represent a list of potentially interesting genes to be studied further (File S8). Many of the entry genes with highest ranking have a local context with a clear single focus. For example, sll0543 forms a cluster with genes encoding three key members of PSI (psaC, psaB, psaD) (Fig. 8B). In contrast, a similar analysis with STRING places sll0543 in a cluster of eight genes annotated as ‘hypothetical protein’ and one as ‘indole-3-glycerol-phosphate synthase’. In another example, the slr0144-48 Pap operon (see case study 3) is once again identified (Fig. 8C). Interestingly, in this search, the Pap operon genes form a cluster together with two PSI subunits (psaB and psaD): the only earlier study linked the Pap operon to PSII, not PSI (Wegener et al., 2008). Other selected findings include unknown genes slr0723 and sll1774 forming an intricate cluster with two genes encoding proteins with a role in pili biogenesis (slr0161, slr0163) and another gene linked to chemotaxis (slr1043). The ‘unknown protein’ slr1187 forms a cluster with three NADH dehydrogenase subunits (slr1279-81) (Fig. 8D), and the ‘hypothetical protein’ slr2003 forms a cluster with two nitrate/nitrite transport system components (slr1450-51) (Fig. 8E).