Information

Why is there a compound highlighted in red in a KEGG Module search?

Why is there a compound highlighted in red in a KEGG Module search?



We are searching data for your request:

Forums and discussions:
Manuals and reference books:
Data from registers:
Wait the end of the search in all databases.
Upon completion, a link will appear to access the found materials.

Considering the following module:
http://www.kegg.jp/kegg-bin/show_module?M00115+C00003
why is the compound C00003 marked in red?


It is just a way of highlighting your search term. If you search M00115+C03722, you will see C03722 in red.


Biochemistry Practical 1

•Genes are composed of exons and introns. Exons are regions retained in the processed mRNA, and are represented by black blocks in the browser, while introns are the regions that are removed during the process of creating the final mRNA, and are represented by lines connecting the blocks.

•The codon ATG in DNA (AUG in mRNA) specifies the amino acid M (Methionine) and is highlighted in green on the "Base Position" track of the Genome Browser. The first Methionine provides the starting signal for protein synthesis.

•The codons TAA, TAG, and TGA in DNA (UAA, UAG, and UGA in mRNA) encode the stop codon (*) and are highlighted in red on the "Base Position" track of the Genome Browser. The stop codons provide the ending signal for protein synthesis.

•Genes may be read either from left to right (top strand of the DNA), or from right to left (bottom strand of the DNA). Arrows on a gene indicate its directionality.


Background

Sunlight is one of the most important abiotic factors for plant growth and development. It can be converted to chemical energy, which is then used for synthesizing organic compounds via photosynthesis altered sunlight conditions can exert a significant influence on the growth and chemical composition of grape berries [1]. Some canopy management practices such as leaf removal, cluster thinning, grapevine training, and leaf moving are widely used to optimize the canopy microclimate, allow varying sunlight exposure, control berry yield, and improve grape berry and wine quality [2]. Among these viticulture practices, leaf removal in a cluster zone (also called basal leaf removal) has been most commonly conducted, primarily because of its ability to promote sunlight exposure and airflow as well as to reduce foliage cover and disease incidence [3, 4]. It has also been found that artificial defoliation has a positive effect on phenolic and volatile compounds in grapes and wine [5, 6].

Leaf removal is generally performed in cool regions with appropriate sunshine and heat accumulation and rainfall [7]. It is typically conducted to selectively or completely strip the foliage from around the bunch zone, and this practice is traditionally implemented at a certain time after fruit set, usually before véraison [6, 8]. In the face of global warming combined with the sensitivity of grape berry ripening to climate change, viticulture management implemented in sunshine- and heat-appropriate regions should be adjusted to adapt to the warming climate [9]. In some strong sunshine and arid regions such as the wine producing regions in northwestern China, grapevine leaf removal in the green-fruit period occasionally causes grape berry sunburn and even leads to lignified and browned stems, which can cause the grape berries to stop growing due to nutrient deficiency. Moreover, the ripening progression of grape berries in this region is always accelerated due to the dry and hot climate [10, 11]. The shorted ripening duration also results in phenolic compound deficiencies, especially anthocyanins and phenolic co-pigments (e.g. myricetin, quercetin, catechin, epicatechin) that are sensitive to changes in climatic conditions and can compromise the color intensity and stability of wine [12]. Accordingly, it is necessary to adjust the timing of cluster sunlight exposure in dry-hot climate viticulture. Our previous study has shown that leaf removal or leaf moving at véraison, which exposes grape clusters to sunlight until harvest, can markedly improve the accumulation of flavon-3-ols and reduce the concentrations of anthocyanins in grape berries grown on the north foot of the Mt. Tianshan region of Xinjiang in northwestern China [5]. The aim of the present study was to dissect the variation in volatile compound metabolome and transcriptome in these exposed grape berries in this dry-hot climate region.

Grape-derived volatile compounds play the most role in evaluating the quality of grapes and wine. Previous studies have reported the effects of basal leaf removal at pre-véraison on the accumulation of monoterpenes and norisoprenoids that contribute to the Muscat varietal aroma and pleasant odor of grape [8, 13, 14]. Moreover, basal leaf removal causes variation in other volatile compounds such as methoxypyrazine [4, 15], thiol [16], and rotundone [17], which impart the vegetal, citrus, and black pepper aromas in grape berries. Indeed, the timing and intensity of sunlight exposure have distinct influences on the volatile compounds produced in grape berries. As Kwasniewski et al. observed [14], only cluster sunlight exposure starting at 33 days past berry set (PBS) significantly increases the concentration of total 1,1,6-trimethyl-1,2-dihydronaphthalene (TDN) and vitispirane, whereas leaf removal at 68 days PBS reduces β-damascenone generation. Additionally, when all basal leaves are removed to completely expose the grape cluster to sunlight, the berries accumulate more β-damascenone and some bound-form terpenoids [6]. Cluster sunlight exposure by apical defoliation approaches, compared with basal leaf removal, can minimally influence wine volatile compounds but reduce wine alcohol content [3]. A limited number of investigations have dealt with the change in volatile C6/C9 compounds in grape berries exposed to sunlight by leaf removal at the early stage of berry development [6, 18, 19] however, the influence of leaf removal at the véraison or ripening stage has not yet been understood. The C6 aldehydes and alcohols can give rise to the characteristic ‘green’ odor, also called ‘green leaf volatiles’ (GLVs). These compounds are induced by the disruption of plant tissues or after plants suffer biotic or abiotic stresses [20]. C9 aldehydes, especially (E)-2-nonenal and (E,Z)-2,6 nonadienal, contribute to the cucumber flavor in plants [21]. Previous studies have also not dealt with the variation in volatile benzenoid-derived compounds in grape berries caused by leaf removal. Such compounds can confer floral and fruity flavors to grape berries and their corresponding wines [22, 23]. Understanding the variation in grape-derived volatile profile benefits an overall evaluation of how leaf removal in regions with intense sunshine and little rainfall will contribute to grape aroma quality improvement strategies.

Leaf removal may eliminate potential assimilated carbon supplements that the fruit receive from neighboring leaves, whereas leaf moving from around clusters enables vines to not only retain the photosynthetic organs but also increase the cluster sunlight exposure. Leaf removing at véraison could significantly promote the accumulation of total anthocyanins and up-regulate related genes [24], but the influence of this performance on the production of volatile compounds remains unclear. Furthermore, a previous transcriptomic study has only focused on the influence of cluster sunlight exposure at the early growth stage of grape berries (E-L 29) [8], whereas the transcriptomic response in grape berries to leaf removal or leaf moving at véraison or the ripening stage is poorly understood.

In this study, four cluster sunlight exposure strategies including leaf removal at the pepper-corn size stage (LR-PS), leaf removal at véraison (LR-V), half-leaf removal at véraison (HLR-V), and leaf moving at véraison (LM-V). A combined analysis of volatile metabolome and transcriptome data was conducted to elucidate the efficiency of these cluster sunlight exposure manipulations on grape berry volatile compound production, and the underlying mechanisms.


Results

Establishment and Validation of the Lasso Penalized Cox Regression Model

The flowchart shows our analysis procedures for the construction of the SE-associated gene-based prognostic model of osteosarcoma ( Figure 1 ). The 349 SE-associated genes based on the dataset of <"type":"entrez-geo","attrs":<"text":"GSE39058","term_id":"39058">> GSE39058 were selected to perform the Lasso penalized Cox regression analysis. Firstly, to identify the best-fit parameter of the λ parameter, two important optional λ values (logλmin = 𢄡.91 and logλlse = 𢄡.48) were calculated from the vertical lines with a minimizing mean-square error and further used to select two groups of genes ( Figures 2A𠄼 ). Moreover, Lasso models were reconstructed according to the λmin and λlse, and survival probabilities were further estimated based on two gene lists ( Figure 2D ). As shown in Figure 2D , the use of the five gene-based prognostic model based on the λmin (Wilcoxon test, p = 1.4e�) facilitated the obvious distinction of survival probabilities between samples obtained from living and expired patients. However, when a two-gene-based prognostic model was used based on the λlse, there was no statistically significant difference observed (Wilcoxon test, p > 0.05). Similarly, the outcome of ROC curve analysis showed that the area under the curve minimum (AUCmin) was 0.92, implying that the five-gene-based Lasso model performed well in predicting the probability of OS of patients with osteosarcoma ( Figure 2E ).

The procedure workflow used to establish and certify the SE-associated gene-based prognostic model for patients with osteosarcoma. KEGG, Kyoto Encyclopedia of Genes and Genomes Lasso, least absolute shrinkage and selection operator Target-OS, Target-osteosarcoma SE, super-enhancer WGCNA, weighted gene co-expression network analysis.

Establishment of the five-gene prognostic model by Lasso regression analysis based on the 349 SE-associated genes from data downloaded from SEdb. (A,B) Lasso coefficient profiles of the 349 SE-associated genes. (C) Two optimal lambda (λ) values (λmin and λlse) were estimated according to the vertical lines with a minimizing mean-square error. The vertical axis represents the mean-square error, while the horizontal axis represents the value of log(λ). (D) The scatter plot of survival status of patients with osteosarcoma based on the five-gene model (left, λmin, p = 1.4e�) or two-gene model (right, λlse, p > 0.05) using the Wilcoxon test. (E) ROC curves of the two prognostic models based on λmin and λlse. The values of AUC were included in the figure. AUC, area under the curve Lasso, least absolute shrinkage and selection operator ROC, receiver operating characteristic SE, super-enhancer SEdb, super-enhancer database.

Validation of Independent Prognostic Factors by the Cox Regression Model

For the validation of the Lasso penalized Cox regression model, we performed univariate and multivariate Cox regression analyses to determine whether these genes were independent prognostic factors for the OS of patients with osteosarcoma. In the univariate Cox regression analysis, all log-rank p-values of these five genes were π.05 ( Figure 3A ). Following multivariate Cox regression analysis, the global p-value (log-rank test) of the five-gene prognostic model was only 0.000171 ( Figure 3B ). The AIC was 59.74, and the C-index was 0.89, indicating that these five genes may be favorable prognostic factors for the OS of patients with osteosarcoma. Moreover, the outcome of K–M survival analysis was consistent with that of the univariate Cox regression analysis ( Figure 3C ). Furthermore, AMN1 and ZP3 may be protective factors in osteosarcoma (HR: 1.26e� and 1.13e�, respectively), whereas LIMS1, SAMD4A, and SPARC appear to be harmful factors in this setting (HR: 699.2, 167, and 298.7, respectively). Thus, a Lasso penalized Cox regression model including five SE-associated genes may be used to predict the OS of patients with osteosarcoma.

Univariate and multivariate Cox regression analyses and K–M survival analysis of the five-gene prognostic model. (A) Univariate Cox regression analysis of each gene of the model. (B) Multivariate Cox regression analysis of the five-gene model. (C) K–M survival curves showing the difference in OS between the relative high- and low-expression groups for each gene according to the median of expression levels. K–M, Kaplan–Meier AIC, Akaike information criterion OS, overall survival.

Establishment and Validation of the Polygenic Risk Score Model

We integrated the expression data of the five genes and corresponding coefficient derived from the above multivariate regression analysis to establish the risk score model. All patients in the training dataset of <"type":"entrez-geo","attrs":<"text":"GSE39058","term_id":"39058">> GSE39058 (N = 41) were divided into high-risk (risk score > 0) and low-risk groups (risk score < 0) ( Figure 4A ). As shown in Figure 4B , survival was more commonly observed in the low-risk group, whereas death was more frequent in the high-risk group. LIMS1, SAMD4A, and SPARC tended to upregulated, whereas AMN1 and ZP3 were downregulated in patients of the high-risk group ( Figure 4C ). The K–M OS analysis predicted an inferior outcome of OS in the high-risk group (log-rank test, p = 0.0006) ( Figure 4D ). Meanwhile, the AUCs of a time-dependent ROC curve calculated using the five-gene-based risk score model were Ϡ.8 ( Figure 4E ), indicating that the forecast model had high sensitivity and specificity.

The five-gene risk score model for the <"type":"entrez-geo","attrs":<"text":"GSE39058","term_id":"39058">> GSE39058 dataset (N = 41). (A) Scatter plot displaying the risk score of each patient in the <"type":"entrez-geo","attrs":<"text":"GSE39058","term_id":"39058">> GSE39058 dataset. Patients were divided into high- or low-risk groups according to the risk score. The blue plots represent the patients in the low-risk group (risk score ≤ 0), while the red plots represent patients in the high-risk group (risk score > 0). (B) The survival status distribution of patients with osteosarcoma in the high- or low-risk groups. (C) The expression status of the five prognostic genes in 41 patients with osteosarcoma in the <"type":"entrez-geo","attrs":<"text":"GSE39058","term_id":"39058">> GSE39058 dataset. (D) K–M survival curves showing the difference in OS between high- and low-risk patients (41 patients) (log-rank test, p = 0.0006). (E) Time-dependent ROC curve analysis for the prediction of survival using the five-gene prognostic model. The AUCs of 1-, 3-, and 5-year OS are shown in the figure. AUC, area under the curve K–M, Kaplan–Meier OS, overall survival ROC, receiver operating characteristic.

The Target-OS dataset was further used to verify the predictive values of the polygenic risk score. Following the multivariate Cox regression analysis, the global p-value (log-rank test) of the five-gene prognostic model for the Target-OS dataset was π.05 (Supplementary Figure S1A). A total of 88 patients were classified into the low- and high-risk groups using the optimal cutoff value of the risk score (Supplementary Figures S1B𠄽). The K–M curves of two groups based on risk scores were significantly different (log-rank test, p < 0.05), and AUCs of time-dependent ROC curves also showed that the five-gene risk model could be a favorable approach to predicting the OS of patients with osteosarcoma (Supplementary Figures S1E,F).

Construction of the Predictive Nomogram Based on the Risk Score Model for Prognostic Prediction

To visualize the Cox regression model, a nomogram was constructed for predicting the 1-, 3-, and 5-year OS probability for patients with osteosarcoma in the <"type":"entrez-geo","attrs":<"text":"GSE39058","term_id":"39058">> GSE39058 dataset (N = 41). The predictors of the nomogram included age, sex, necrosis, recurrence, and risk group ( Figure 5A ). The calibration plot shown in Figure 5B illustrates the performance of the nomogram. The calibration plots were close to the ideal prediction gray line (45° line), indicating that our nomogram performed well in predicting survival probability ( Figure 5B ). To assess the predictive effect of the risk group as an indicator of OS, we applied the nomogram to a specific patient ( <"type":"entrez-geo","attrs":<"text":"GSM954850","term_id":"954850">> GSM954850) in <"type":"entrez-geo","attrs":<"text":"GSE39058","term_id":"39058">> GSE39058 ( Figures 5C,D ). The patient <"type":"entrez-geo","attrs":<"text":"GSM954850","term_id":"954850">> GSM954850 expired at 750 days. According to the model containing the risk group, the predicted probability of death at 1,095 days was 0.815 this value was markedly higher than that obtained for the model lacking the risk group (0.659). This finding implies that the predictive model containing the risk group as a parameter is more accurate than that lacking the risk group.

Nomogram predicting the probability of 1-, 3-, and 5-year OS in patients with osteosarcoma. (A) Nomogram adding up the points identified on the points scale (the upward line) for each variable. The total points projected on the bottom scales indicate the probability of 1-, 3-, and 5-year OS. (B) Calibration plot for predicting the 1-, 3-, and 5-year OS. The gray line represents the ideal condition. (C,D) The nomograms predicting the probability of 1-, 3-, and 5-year OS for the specific patient <"type":"entrez-geo","attrs":<"text":"GSM954850","term_id":"954850">> GSM954850 based on the model containing or not containing the risk group in the <"type":"entrez-geo","attrs":<"text":"GSE39058","term_id":"39058">> GSE39058 dataset. OS, overall survival DFS, disease-free survival.

Identification of SE-Related Hub Genes in Osteosarcoma Using WGCNA and a Protein–Protein Interaction Network

Weighted gene co-expression network analysis is another bioinformatics method for the analysis of clinical traits linked to the levels of gene expression. Cluster analysis was performed to display the heatmap of clinical traits ( Figure 6A ). As shown in Figure 6B , β = 9 was selected as the best soft-thresholding value via prediction of the scale independence and mean connectivity ( Figure 6B ) to construct the gene co-expression network. Subsequently, eigengenes were divided into different colored modules, producing 17 different modules. A module–trait relationship heatmap was developed according to Pearson’s correlation coefficient ( Figure 6C ). Three modules presented higher correlation with clinical traits: MEgray60 with recurrence (r = 0.37, p = 0.01) MEblue with death (r = 0.33, p = 0.03) and MEpink with survival time (r = 0.47, p = 0.002) ( Figure 6C ). A scatter plot also showed that the MM of these three modules was positively correlated with GS ( Figure 6D ). A Venn plot of SE-associated genes derived from the SE-associated gene matrix was drawn, and the genes within three modules were intersected ( Figure 6E ). There were 13, 70, and 16 interaction genes between module gray60, blue, pink, and SE-associated genes. The interaction genes were enriched in numerous cancer-related pathways, such as proteoglycans in cancer and transcriptional misregulation in cancer ( Figure 6F ). In addition, we constructed a gene𠄽rug interaction network using the five genes and related drugs (Supplementary Figure S2). All five genes interacted with JQ1 (an established SE inhibitor). The network depicted that these five genes were probably regulated by SE to a certain extent. Finally, we displayed the signal tracks for the H3K27ac, H3K4me1, and H3K4me3 ChIP–seq profiles of the five genes ( Figure 7 ). We observed five predicted SEs near these five genes. These data suggest that the expression of the five genes is regulated by SE. In addition, the SE inhibitor JQ1 may regulate the expression pattern in U2OS cells.

Identification of SE-related hub genes in osteosarcoma based on the <"type":"entrez-geo","attrs":<"text":"GSE39058","term_id":"39058">> GSE39058 dataset through WGCNA analysis. (A) Cluster analysis between the gene expression in the <"type":"entrez-geo","attrs":<"text":"GSE39058","term_id":"39058">> GSE39058 dataset and sample outliers or clinical traits. Color intensity was proportional to sample outliers, age, gender, necrosis, recurrence, and survival time. (B) Analysis of the scale independence and mean connectivity (vertical axis) for various soft-thresholding powers (β value of horizontal axis). (C) Heatmap of the correlation between modules and clinical traits of osteosarcoma p-values in the table specify the correlation between modules and clinical traits. Three modules had high correlation with clinical traits: MEgray60 with recurrence (r = 0.37, p = 0.01) MEblue with death (r = 0.33, p = 0.03) and MEpink with survival time (r = 0.47, p = 0.002). (D) Scatterplot showing the MM versus GS of module genes related to death in the blue module, recurrence in the gray60 module, or survival time in the pink module (horizontal axis: MM means module membership, vertical axis: GS means gene significance). (E) Venn plot of SE-related genes and eigengenes in the blue, gray60, and pink modules. There are 70, 13, and 16 correlated genes, respectively. (F) KEGG pathway analysis of the above-correlated genes in the Venn plot. KEGG, Kyoto Encyclopedia of Genes and Genomes SE, super-enhancer WGCNA, weighted gene co-expression network analysis.

Signal tracks for H3K27ac (red), H3K4me1 (purple), and H3K4me3 (green) ChIP–seq profiles of the five-SE-associated hub genes visualized using IGV. The regions of SE are shown in a pink bar upon the signal tracks. ChIP–seq, chromatin immunoprecipitation–sequencing SE, super-enhancer IGV, Integrative Genomics Viewer.


Results

Detection of taxonomic biomarkers

To identify core candidate microbiota biomarkers that are present in Crohn’s disease and healthy samples, a cladogram was constructed to demonstrate relative abundance of bacteria. Using LEfSe tool, we identified 40 differential abundant microbial taxonomic features in control samples, stool samples and mucosal tissue samples from CD. Small circle on the cladogram ring represents a taxonomic rank, which has different abundance values among the groups based on the LDA scores. All detected microbial taxonomic features can be presented in cladogram highlighting significant differences across three types of samples (See Fig.  3 (top)). We specifically discuss the results from family and genera biomarkers. The LEfSe analysis found Streptococcaceae, Lactobacillales, and Pseudomonadaceae are differentially abundant in the CDS, whereas Porphyromonadaceae, Shewanellaceae, and Enterobacteriaceae are differentially abundant in CDT. Bacteroidaceae, Lachnospiraceae, Rikenellaceae, and Ruminococcaceae were identified as taxonomic biomarkers for healthy individuals.

Cladograms generated from LEfSe for biomarker detection in taxonomic (top) and metabolic function pathways (bottom)

Detection of metabolic functional biomarkers

In addition to microbial composition, we also compared differentially abundant functional and metabolic characteristics in three microbial samples. Figure  3 (bottom) highlights 135 differentially abundant functional modules detected in the microbial communities corresponding to CDT, CDS and HCS. While various microbial metabolic functions are carried out throughout the human microbiome, specific subsets of this functionality could be enriched in different types of samples. The LEfSe tool highlights these specific metabolic features (KEGG modules) as shown in Fig.  3 (bottom). Modules such as biosynthesis of lysine (M00016), and UMP (M00051) were differentially enriched in healthy control samples. We also found that the glutathione biosynthesis (M00118), metabolism of the sulfur-containing amino acids cysteine (M00338), and methionine biosynthesis (M00017) were significantly enriched in CDT. In addition, several other modules essential for basic life activities of prokaryotic cells, such as central carbohydrate metabolism ( <"type":"entrez-nucleotide","attrs":<"text":"M00002","term_id":"202949","term_text":"M00002">> M00002-M00007) and amino acid metabolism (M00018, M00019, M00020, M00118 and M00338) are highlighted in the cladogram. These results exclusively show that the specific metabolic modules are enriched in distinct biological samples.

Detection of bacterial interactions

We explored the inter-bacterial association networks at the family and genus level in three environments (CDT, CDS, and HCS). Table  1 and Additional fileਁ present the results of the positive and negative associations among bacteria by a Spearman’s correlation approach. Less number of associations were identified in stool samples from Crohn’s disease patients and healthy individuals as opposed to mucosal tissues from Crohn’s disease samples.

Table 1

Inter-bacteria correlations in all sample groups

Taxonomic cladeTaxonomic cladeR 2
CDSf__Bacteroidaceaef__Lachnospiraceae0.94
f__Aerococcaceaef__Fusobacteriaceae0.98
HCSf__Prevotellaceaef__RF160.98
f__Bacillaceaef__Staphylococcaceae0.98
f__Rikenellaceaef__Ruminococcaceae0.97
CDTf__Aeromonadaceaef__Shewanellaceae0.81
f__BA059f__Syntrophobacteraceae0.72
f__Planococcaceaef__Gallionellaceae0.72
f__Porphyromonadaceaef__Pseudomonadaceae0.71
f__Carnobacteriaceaef__Streptococcaceae0.70
f__Moraxellaceaef__Pseudomonadaceae0.68
f__Microbacteriaceaef__Spirochaetaceae0.68
f__BA059f__Gallionellaceae0.68
f__Peptococcaceaef__Alteromonadaceae0.68
f__Peptococcaceaef__Sinobacteraceae0.68
f__Nitrospiraceaef__Syntrophobacteraceae0.68
f__Procabacteriaceaef__Halomonadaceae0.68
f__Veillonellaceaef__Pseudomonadaceae-0.66
f__Porphyromonadaceaef__Shewanellaceae0.66

In mucosal tissues from Crohn’s disease patients, 13 positive with one negative relationships were recognized between bacterial families. A strong positive correlation between Aeromonadaceae and Shewanellaceae was observed in Table  1 and both have common ancestor in their evolutionary lineages. There were also strong positive and negative associations between bacterial genera in CDT and CDS. Like CDT, significantly strong positive interactions were observed in CDS and healthy individual samples. All observed inter-bacterial associations at the genus level have shown high correlation with one another in our result. For example, Prevotellaceae with RF 16, and Bacillaceae with Staphylococcaceae, are highly correlated along with shared evolutionary lineage. Hence, these results demonstrate that there are differences in microbial interactions between CD patients and HCS. Similar differences in bacterial relationship were reflected in other sample groups.

Detection of associations between bacterial taxa and microbial pathway

For all the bacteria with strong associations in the previous results, we identified their highly correlated KEGG orthologues (Tables  2 and ​ and3 3 ).

Table 2

Identifying associations between bacterial families and their microbial pathways with KO density in Crohn’s Disease Stool

Table 3

Identifying associations between bacterial families and their microbial pathways with KO density in Crohn’s Disease Tissue

KEGG ortholog (KO)ModuleDensity
f_Pseudomonadaceaef_Moraxellaceae <"type":"entrez-nucleotide","attrs":<"text":"K00404","term_id":"162411","term_text":"K00404">> K00404, <"type":"entrez-nucleotide","attrs":<"text":"K00405","term_id":"162413","term_text":"K00405">> K00405, <"type":"entrez-nucleotide","attrs":<"text":"K00406","term_id":"162414","term_text":"K00406">> K00406, <"type":"entrez-nucleotide","attrs":<"text":"K00407","term_id":"162416","term_text":"K00407">> K00407M001560.80
f_Pseudomonadaceae <"type":"entrez-nucleotide","attrs":<"text":"K01856","term_id":"174057","term_text":"K01856">> K01856, <"type":"entrez-nucleotide","attrs":<"text":"K03464","term_id":"197513","term_text":"K03464">> K03464, <"type":"entrez-nucleotide","attrs":<"text":"K01055","term_id":"172411","term_text":"K01055">> K01055, <"type":"entrez-nucleotide","attrs":<"text":"K03381","term_id":"210219","term_text":"K03381">> K03381M005680.80
f_Moraxellaceae <"type":"entrez-nucleotide","attrs":<"text":"K01856","term_id":"174057","term_text":"K01856">> K01856, <"type":"entrez-nucleotide","attrs":<"text":"K03464","term_id":"197513","term_text":"K03464">> K03464, <"type":"entrez-nucleotide","attrs":<"text":"K01055","term_id":"172411","term_text":"K01055">> K01055M005680.60
f_Pseudomonadaceaef_Moraxellaceae <"type":"entrez-nucleotide","attrs":<"text":"K00457","term_id":"343194","term_text":"K00457">> K00457, <"type":"entrez-nucleotide","attrs":<"text":"K00451","term_id":"208933","term_text":"K00451">> K00451, <"type":"entrez-nucleotide","attrs":<"text":"K01800","term_id":"161416","term_text":"K01800">> K01800, <"type":"entrez-nucleotide","attrs":<"text":"K01555","term_id":"174453","term_text":"K01555">> K01555M000440.67
f_Pseudomonadaceaef_Moraxellaceae <"type":"entrez-nucleotide","attrs":<"text":"K00166","term_id":"175799","term_text":"K00166">> K00166, <"type":"entrez-nucleotide","attrs":<"text":"K00167","term_id":"174925","term_text":"K00167">> K00167, K09699, <"type":"entrez-nucleotide","attrs":<"text":"K00253","term_id":"175764","term_text":"K00253">> K00253, <"type":"entrez-nucleotide","attrs":<"text":"K00249","term_id":"176448","term_text":"K00249">> K00249, <"type":"entrez-nucleotide","attrs":<"text":"K01968","term_id":"508622","term_text":"K01968">> K01968, <"type":"entrez-nucleotide","attrs":<"text":"K01969","term_id":"171834","term_text":"K01969">> K01969,K1376M000360.62
f_Pseudomonadaceaef_MoraxellaceaeK02274, <"type":"entrez-nucleotide","attrs":<"text":"K02275","term_id":"188926","term_text":"K02275">> K02275, <"type":"entrez-nucleotide","attrs":<"text":"K02276","term_id":"188927","term_text":"K02276">> K02276M001550.60

Several KEGG orthologues related to V/A-type H+/Na+ transporting ATPase subunit A ( <"type":"entrez-nucleotide","attrs":<"text":"K02117","term_id":"146697">> K02117), B ( <"type":"entrez-nucleotide","attrs":<"text":"K02118","term_id":"151839","term_text":"K02118">> K02118), C( <"type":"entrez-nucleotide","attrs":<"text":"K02119","term_id":"153938","term_text":"K02119">> K02119), D( <"type":"entrez-nucleotide","attrs":<"text":"K02120","term_id":"210767","term_text":"K02120">> K02120), E( <"type":"entrez-nucleotide","attrs":<"text":"K02121","term_id":"330029","term_text":"K02121">> K02121), I( <"type":"entrez-nucleotide","attrs":<"text":"K02123","term_id":"335123","term_text":"K02123">> K02123), and K( <"type":"entrez-nucleotide","attrs":<"text":"K02124","term_id":"216170","term_text":"K02124">> K02124) showed positive correlation (Spearman’s correlation Ϡ.6, FDR π.05) with Bacteroidaceae and Lachnospiraceae in CDS (Table  2 ). These strong correlations between the abundances of bacteria taxon and gene abundances (KO) highlight genes relevant to disease phenotype in the bacterial species. A V-type ATPase in prokaryotes (M00159) KEGG module was highly associated (KO density Ϡ.6) with above mentioned KOs in the stool samples from Crohn’s disease patients. In the CDT, Table  3 shows Pseudomonadaceae and Moraxellaceae were found to be positively correlated with several genes (KO). For those significant associations between the taxonomic clades and metagenomic gene familes, 5 strongly associated KEGG modules, viz. Cytochrome c oxidase, cbb3-type (M00156), Catechol ortho-cleavage, catechol ⇒ 3-oxoadipate (M00568), Tyrosine degradation, tyrosine ⇒ homogentisate (M00044), Leucine degradation, leucine ⇒ acetoacetate + acetyl-CoA (M00036) and Cytochrome c oxidase, prokaryotes (M00155), were identified. Additional fileਂ shows the associations of all correlated bacterial genera with their highly correlated KEGG orthologues in CDT. Those three bacteria revealed strong associations with four KEGG modules, viz. Polyamine biosynthesis (M00134), Nucleotide sugar biosynthesis (M00554), PRPP biosynthesis (M00005), and Trans-cinnamate degradation (M00545).

Split graph analysis

The resulting split graph consists of two disjoint sets of nodes, where one set corresponds to correlated microbial communities, and the other set corresponds to their microbial metabolic pathways. We automatically extracted various important components (subgraphs) from the split graphs that model the integrated network in both samples (CDS and CDT). Again, the automatic extraction of such components is implemented by finding high-weighted maximal cliques in the split graph. Due to the independence of the nodes representing the pathways, each clique in the graph contained one node representing a pathway. A high-weighted clique is the graph which will contain a group of bacteria that are highly correlated and a pathway that is highly associated or impacted by such group.

In the CDS split graph, two bacteria at the family level, Bacteroidaceae and Lachnospiraceae, are highly correlated with each other (Spearman’s correlations 0.94, FDR π.05). This clique is associated with V-type ATPase, prokaryotes KEGG module (M00159) (Fig.  4 ). The quantified values for this association were weighted based upon calculating the proportions of KEGG orthologus (KO) for each correlated bacteria. Similarly, in the CDT split graph, a maximal clique of size two was identified with high correlation between Pseudomonadaceae and Moraxellaceae (Spearman’s correlations 0.68, FDR π.05) (Fig.  5 ) and multiple KEGG modules were connected to this clique of bacteria (KO density Ϡ.6). These KEGG modules in CDT are mainly involved in ATP synthesis and amino acid metabolism. Yet, the extent to which bacteria in the clique correspond to distinct functional modules in the split graph has remained largely unclear. We also extracted the split graph where the microbial components were considered at the genus level. At the genus level, we obtained multiple split graphs with different sizes of clique where each pair of bacteria are highly correlated (Spearman’s correlations Ϡ.6, FDR π.05) in CDT (See Additional fileਃ). Figure  6 represent two complete split graphs containing multiple high-weighted maximal cliques. For instance, three bacterial components correlated with PRPP biosynthesis microbial metabolic pathway constitute one of the high-weighted maximal cliques.

Split graph in Crohn’s disease stool samples at the family taxonomic level

Split graph in Crohn’s disease tissue samples at the family taxonomic level

Top two split graphs in Crohn’s disease tissue samples at the genus taxonomic level

We also visualized a heatmap of OTU abundances at the genus level to assess the abundance of bacteria in the samples (Fig.  7 ). The abundance of Blautia and unknown genus from Lachnospiraceae family and that of unknown genus from Ruminococcaceae family were observed to be high in CDT samples. Likewise, the bacterial genera, Veillonella and Bacteroides, were also highly abundant in CDT group. In addition, there were no high-weighted maximal cliques obtained in split graph from CDS samples as none of the KEGG modules were significantly correlated to any of the highly correlated bacteria in CDS (Spearman’s correlation Ϡ.6, FDR π.05).

Heatmap of relative abundance of 23 bacterial genera in Crohn’s disease Stool (CDS), Crohn’s disease tissue (CDT) and control samples (HCS)

Comparison of Two Population Proportions Analysis

In the split graph, the extracted high-weighted maximal cliques of microbial communities with their metabolic pathways were mostly observed in the microbiome profiles obtained mucosal tissue samples of CD patients.

Among all correlated edges in CDT and CDS networks, proportion of the correlated edges with a common family is different between CDT and CDS networks (Table  4 ). In other words, the proportion of correlated edges in CDT with common family is significantly different to the proportion of correlated edges in CDS with common family. Similarly, the proportion of correlated edges with common class (and phylum) is also significantly different between CDT and CDS networks.

Table 4

p-value from proportion test at different level of taxonomy

Previous studies on microbiomes in Crohn’s disease revealed that fecal bacterial ecosystems differ from those in the intestinal mucosal tissue [22, 23]. Studies of microbiome in fecal samples have more challenges in identifying their community associated with respect to disease initiation and progression due to the nature of the environment. Based on these observations, we can infer that microbial dysbiosis is less tended to be shifted toward lumens in a given disease state. In order to gain a better understanding of possible microbial mechanisms, the need to examine tissue biopsies along with stool samples are highlighted.


PathBank LAYOUT AND NAVIGATION

A screenshot montage of the PathBank interface and its search and browsing tools is shown in Figure 1. PathBank has a landing/home page very similar to that of its smaller cousin, SMPDB ( 15). This includes a short textual description of the database, a circulating ‘carousel’ of selected PathBank images and a series of tabs at the top of the page for navigation. An empty search box is located at the top left of each PathBank page. PathBank has five major tabs: (i) Browse (ii) Search (iii) About (iv) Downloads and (v) Contact Us. Under the Browse tab users may browse the website using four different pages: a) Pathways b) Table of Primary Pathways (TOPP) c) Compounds (metabolites or drugs) or d) Proteins. All pathway search and browse results are filterable by species (10 model species available) and by pathway type (Metabolite/Compound and Protein) using two dropdown filters. Additionally, selecting a pathway type will render visible a third dropdown menu to allow filtering of the results by pathway sub-categories (6 Metabolite/Compound pathway sub-categories and 15 Protein pathway sub-categories).

A screenshot montage of different browsing and searching screens taken from PathBank. A more detailed description of the different functions and capabilities of the various browse and search tools in PathBank is given in the text.

A screenshot montage of different browsing and searching screens taken from PathBank. A more detailed description of the different functions and capabilities of the various browse and search tools in PathBank is given in the text.

If users choose to browse by Pathway, a multi-page table is presented that contains four columns: (i) the PathBank identifier (with the thumbnail image) (ii) the pathway details, which includes the pathway name, species, description and download links (iii) the compounds (which are hyperlinked) and (iv) the proteins (which are also hyperlinked). Clicking on the thumbnail image of a given pathway will direct users to the full-size pathway image. Likewise, clicking on the metabolite/compound names (column iii) or protein names (column iv) will open a page with detailed compound descriptions (from the appropriate organism-specific metabolome database such as HMDB or from a specialized organism-specific database) or protein descriptions (from UniProt).

If the TOPP option is selected, a multi-page table is presented showing three columns: (i) the PathBank identifier (ii) the pathway name including the species name (if filtered by all species) and (iii) the link to the pathway view. In the TOPP, similarly named (replicated) pathways are grouped together to facilitate browsing of largely unique pathways. Groupings are represented by a single primary pathway, hence the table name, which has a clickable pathway name displayed in blue text. These hyperlinked pathways may be selected by the user to navigate to a popover view containing all of the related pathways which can be similarly browsed. In the case of lipid pathways, only the generic version of the pathway is listed in the main table the pathways for the specific lipids are browsable in the generic lipid's popover view.

Users can browse PathBank's compounds and proteins in much the same manner. Multi-page tables present compounds and proteins alphabetically in three columns: (i) the compound/protein's identifier along with the corresponding View buttons (ii) the compound/protein's name (with a short description) and (iii) the associated pathways (which are hyperlinked). Clicking on a compound's View button will take a user to the appropriate ‘Compound Card’, corresponding to the database most appropriate for the chosen (or filtered) organism. The default is HMDB if no organism has been selected. Clicking on a protein's corresponding View button will take a user to the appropriate ‘UniProt Card’. Unlike the pathway browse results, compounds and proteins are filterable only by species using the dropdown menu provided. An additional filter in the form of a Search box at the top right of the table allows users to filter the page by compound, protein or pathway name.

PathBank places a strong emphasis on providing users with high quality, artistically pleasing pathway diagrams that are not only correct and informative but also colorful, interactive and richly detailed. All of the pathway diagrams in PathBank use scalable vector graphics (SVG) and a web interface technology inspired by Google Maps. This allows rapid and continuous zooming using a mouse pad or a mouse scroll wheel or through simply clicking on-screen zoom icons. It also allows facile navigation around zoomed-in pathway diagrams through a simple click-and-drag operation or through clicking on-screen up/down or left/right arrows located near the on-screen zoom functions. A full-screen view of each PathBank pathway diagram is also available, which can be toggled off and on by clicking the full-screen icon located between the navigation arrows.

On the right side of each pathway diagram is a pathway display panel with five tabs (Description, Highlight, Analyze, Downloads, Settings/Display). The default view is the Description panel, which describes the pathway in detail and provides one or more literature references. The Highlight panel (viewed by clicking the ‘Highlight’ tab) allows users to select and color different proteins and metabolites for display purposes. The Analyze panel (viewed by clicking the ‘Analyze’ tab) allows users to enter concentration (relative or absolute) data on proteins/transcripts and/or metabolites and to have these colored on the pathway diagram. It is highly recommended that users switch to the black and white version of the pathway when using the highlight and analyze tools in order to clearly see the elements being colored. The download panel (viewable by clicking the ‘Download’ tab) allows users to select the type of file format that they wish to have their pathway image (annotated or untouched) saved as. Options are available for downloading images in BioPAX ( 22) image format (full color SVG, grayscale SVG, simple SVG, large font SVG, and simple large font SVG), BioPAX only format, SBGN ( 23), SBML ( 24), PWML (a custom PathWhiz format), PNG or all of the above. At last, the Display panel (viewed by clicking the ‘Display’ tab) allows users to change the pathway display to suit a particular taste or application. Users may toggle between simple (thick line) and complex (individual lipid) membranes between a dark blue (aqueous) background or a white background between the full color, highly detailed version and the simplified KEGG-like (black and white) representation and finally between full-names of pathway components and their large font abbreviations. There is also an additional fifth option to render the simplified (KEGG-like) pathway representation using the easier-to-read large font abbreviations.

In addition to these browsing functions, PathBank also offers extensive search functions. These include a general text search (available at the top of every PathBank web page) as well as several, more specific, search functions listed under the ‘Search’ tab. The text search supports Boolean logic (AND, OR and NOT operations) along with field-specific searches covering compound names, protein names, compound/protein identifiers, pathway names, and pathway descriptions. Instructions for the text search are accessible using the Search dropdown list. There are four additional searches available: (i) the Path-MAP Advanced Search (ii) the ChemQuery Structure Search (iii) the Molecular Weight Search and (iv) the Sequence Search. The Path-MAP search allows users to enter long lists of compound names, protein names, compound identifiers, protein identifiers or gene identifiers to search for organism-specific pathways enriched with these entities. The result is a list of pathways with enrichment scores calculated using the frequency of matches and the number of pathway entities as scaled using a hypergeometric function. It is also possible to enter similar lists with Path-MAP but with concentrations or relative concentrations, as might be obtained from a typical proteomics, transcriptomics or metabolomics experiment. The result of this type of search is a pathway annotated and colored with the corresponding metabolite/protein concentrations according to a yellow (low)-red (high) concentration gradient. Other searches such as the ChemQuery and Molecular Weight searches are intended for metabolite or compound searching against PathBank's chemical database of 78 271 compounds and are identical to those used by HMDB, SMPDB and many other databases. The Sequence search uses BLAST ( 25) to find sequence matches or sequence similarities to the protein sequences within PathBank's sequence database of 8973 proteins.

PathBank's ‘About’ tab provides additional information about PathBank, its associated release notes, the required citations, database statistics, the PathBank style guide, links to other Pathway Databases, as well as images (via a Pathway legend) for the different components (organelles, organs, tissues) seen in PathBank pathways. PathBank's ‘Download’ tab allows users to download pathways, metabolite names and protein names (in CSV files) as well as all of its pathways in BioPAX, SBGN, SBML and PWML format. A subset of PathBank's pathways (i.e. the primary pathways browsable through the TOPP) are also downloadable in SVG and PNG format. Due to the extremely large file size, users who wish to download all pathway images will need to submit a request using PathBank's Contact form. All of PathBank's sequences (both gene and protein) are available in FASTA text format, all of its chemical structures are available in SDF format and all of its reactions are available in RXN format.


Release Notes

Current release (MBROLE version 2.0) analyzes annotations from the following releases of the databases:

    (Kyoto Encyclopedia of Genes and Genomes) release 54.1. (Human Metabolite Database) version 3.5. (Yeast Metabolite Database) version 1.0. release July 2014. version 17.1 Updated Mar 18, 2014. (Medical Subject Heading) Updated Feb 25, 2014. version 4.0. Release Jan, 2014. (Chemical Entities of Biological Interest) release 116 (June, 2014). (Chemical Entities of Biological Interest) release 47 (April, 2014). version 1.0. (E.coli Metabolome Database) version 1.0.

Background

One important goal of computational systems biology is the development of mathematical models for complex metabolic reaction networks. The type of such models and their predictive capacity depend on the available biochemical knowledge. Generally, one may distinguish two main steps in network modelling: (i) stoichiometric network reconstruction and (ii) network analysis (see figure 1). Necessary prerequisite for any type of mathematical network model is the knowledge of the network stoichiometry, i. e. reactions, membrane transport processes and associated metabolites. For eukaryotic cells, stoichiometric reconstruction of the network includes the assignment of reactions to cellular compartments. Based on the stoichiometric matrix of the network one may perform structural analyses as, for example, identification of elementary flux modes [1], possible routes for a self-consistent expansion of the network starting with some initial seed compounds [2] or calculation of stationary flux distributions based on constraint optimization [3–5]. One central issue in such analyses is to define the possible directionality of a cellular reaction. To decide on this one has to know the (Gibb's) standard free energy change of the reaction and the range within which the ligand concentrations may vary. Compared to structural modelling, a deeper insight into the regulation of reaction networks in response to external variations can be achieved on the basis of kinetic models. For the establishment of kinetic network models rate equations for all reactions and transport processes need to be set up. SABIO-RK and the Chemical Kinetics Database NIST collect kinetic data from the literature. Moreover, application of structural modelling approaches to complex eukaryotic cells requires information on the localization of the reactions in the various intracellular organelles as well as the transport of metabolites among the organelles. The biochemical information required for stoichiometric network reconstruction is spread over various data sources.

Intercorrelation of METANNOGEN with knowledge-bases and network modelling tools. The user manually collects information from biochemical reactions. The network is exported as a SBML file and may be processed with network modelling/analyzing tools.

Comprehensive collections for biochemical reactions are KEGG [6], BIOCYC, [7], REACTOME [8] and UM-BBD [9]. For information on cellular compartmentalization and substrate specificity the enzyme database BRENDA [10] is valuable. Sometimes the cellular compartment of an enzyme reported in databases or publications might not be the site of its action. For example, a mitochondrial enzyme might be firmly attached to the mitochondrial wall and catalyses the biochemical reaction not inside but outside the mitochondrial matrix. Therefore literature search is necessary to obtain detailed knowledge on the biochemical reactions under consideration. There are several approaches to combine information from many sources. The AMAZE LIGHTBENCH combines a variety of information which is accessible with a web browser [11]. Available modelling programs usually allow to edit biochemical networks with a graphical network editor. This approach is sufficient for small networks but for large networks comprising many metabolic reactions an information storage system is required. The recently developed database system META-ALL allows users to enter data on biochemical reactions into a locally running ORACAL database with Web-clients [12]. The resulting model can be exported in SBML format. The program EPE is a visual editor for biological networks including metabolic networks and allows to add annotations.

To facilitate the exploitation of various information sources for the stoichiometric reconstruction of metabolic networks we have developed the interactive computer program METANNOGEN (Figure 1). In contrast to other programs for building user defined reaction network METANNOGEN uses the database KEGG of biochemical reactions as primary information source from which biochemical reactions relevant to the considered network can be selected, edited and stored. The advantages of our approach are that (I) only reaction equations not stored in KEGG need to be entered manually, (II) that the graphical pathway representation of KEGG are used, which look familiar because their layout resembles figures in text books on Biochemistry,(III) that the KEGG database can be updated without affecting the user data and (IV) that the immutable identifiers are provided for compounds and reactions by KEGG.


Drug Repurposing: A Network-based Approach to Amyotrophic Lateral Sclerosis

The continuous adherence to the conventional “one target, one drug” paradigm has failed so far to provide effective therapeutic solutions for heterogeneous and multifactorial diseases as amyotrophic lateral sclerosis (ALS), a rare progressive and chronic, debilitating neurological disease for which no cure is available. The present study is aimed at finding innovative solutions and paradigms for therapy in ALS pathogenesis, by exploiting new insights from Network Medicine and drug repurposing strategies. To identify new drug-ALS disease associations, we exploited SAveRUNNER, a recently developed network-based algorithm for drug repurposing, which quantifies the proximity of disease-associated genes to drug targets in the human interactome. We prioritized 403 SAveRUNNER-predicted drugs according to decreasing values of network similarity with ALS. Among catecholamine, dopamine, serotonin, histamine, and GABA receptor modulators, as well as angiotensin-converting enzymes, cyclooxygenase isozymes, and serotonin transporter inhibitors, we found some interesting no customary ALS drugs, including amoxapine, clomipramine, mianserin, and modafinil. Furthermore, we strengthened the SAveRUNNER predictions by a gene set enrichment analysis that confirmed modafinil as a drug with the highest score among the 121 identified drugs with a score > 0. Our results contribute to gathering further proofs of innovative solutions for therapy in ALS pathogenesis.

This is a preview of subscription content, access via your institution.


Methods

Selection of samples and reads filtering

Illumina sequences were downloaded from Sequence Read Archive (SRA), MG-RAST or JGI Genome portal databases. Quality check and adaptors removal were performed using Trimmomatic (v0.33) and bbduk (version released Nov 2016) (https://jgi.doe.gov/data-and-tools/bbtools/). The composition of the feedstocks used in the different reactors was approximated using substrate information from various sources (Additional file 1). When available, metadata were taken from the publicly accessible description of the respective experiments or full-scale plant operation datasets. Otherwise, reactor feedstock compositions were estimated from the available literature, and were expressed in terms of carbohydrate, protein, lipid and VFA fractions relative to their total solid (TS) content.

Assembly

Reads were assembled using Megahit (v1.1.1) with “−sensitive” mode for samples having less than 40 Gb of sequenced bases and with “–large” for the remaining [50]. Quality of the assemblies was determined using QUAST (v3.1) [51] and the results are reported in Additional file 8.

Binning

Using MetaBAT 2 (v2.12.1) bam files were inspected and each assembly was binned using standard parameters [52]. Minimum size of scaffolds considered for MAGs generation was 1.5 kbp. MAGs were checked for completeness (Cp) and contamination (Ct) using the “Lineage_wf” workflow of CheckM (v1.0.3) [53] and the result obtained for each MAG was determined using the formula: CC3 = Cp − (Ct*3). Removal of contamination from MAGs was performed using RefineM (v0.0.23) [54]. Threshold values used for defining the quality level of MAGs and to assign them to the categories “High Quality” (HQ), “Medium–High Quality (MHQ), “Medium Quality” (MQ) and “Low Quality” (LQ) were defined according to the standards recently described, except for the introduction of the MHQ class (Table 1) [55].

MAGs de-replication

MAGs obtained were de-replicated using Mash (v2.0) [56] on the entire genome sequences with very permissive parameters (0.05 Mash-distance, roughly equivalent to 0.95 ANI and 100/1000 Matching-hashes). Subsequently, a more precise analysis was performed applying the genome-wide Average Nucleotide Identity metric (ANI) using protein-encoding nucleotide sequences only [57]. MAGs were considered as belonging to the same species if they showed ANI value higher than 95% and reaching at least 50% of genome coverage for both strains (on at least one of the two comparisons, “MAG1 vs. MAG2” or “MAG2 vs. MAG1”). Details regarding the assembly and binning procedure are reported in Additional file 2.

Taxonomic assignment

Taxonomic classification was determined for 1635 MAGs obtained after de-replication and belonging at least to the MQ level. This approach was carried out as described previously [4] and more details can be found in the Additional file 2. MAGs were classified by comparison against all taxonomically classified taxa of the NCBI Genome Database (prokaryotic section) using Microbial Genomes Atlas MiGA Online [58].

MAGs coverage calculation and relative abundance

Filtered shotgun reads randomly selected from each sample were aligned back to the entire collection of MAGs. Ordered “bam” files were inspected using CheckM [53] to calculate both the fraction of reads aligned and the relative abundance of each MAG. Analysis was performed using all reads available for each sample and verified using a representative subsample of one million reads per sample. Results obtained using the two datasets of sequences were highly similar (Pearson correlation coefficient was > 0.999 on MAGs representing more than 0.001% of the population). Results obtained using one Mread per sample are reported in Additional file 8. The value (0.001%) was also defined as the arbitrary threshold for considering one MAG as “present in a specific sample”. Coverage values obtained for each MAG were clustered with MeV (v4.9.0) using Pearson correlation and average linkage [59]. The fraction of MAGs shared between different samples was visually represented using CIRCOS (v0.69) [60]. Alpha and beta diversity were determined from the file reporting the number of reads per MAG using Past (v3.21) [61]. The same tool was used for statistical tests and graphical plots.

Gene finding and annotation

Gene annotation was performed using three different procedures: (1) rapid annotation using subsystem technology (RAST annotation server) [62]. These results were reported in a table for comparative purposes (Additional file 14). (2) KEGG annotation and modules completeness were determined using “KEGG Mapping/Reconstructmodule.py” (https://github.com/pseudonymcp/keggmapping). Software assigned to the KEGG modules the results obtained from diamond (v0.9.22.123) alignment only results having max log e-value 1e−5, min bitscore 50, min identity 25 were recovered. Abundance of all the KEGG modules in each experiment was calculated with custom perl scripts (https://sourceforge.net/projects/perl-scripts-kegg/). Cluster analysis on “complete” or “1 bm” KEGG modules identified in HQ and MHQ MAGs was performed using MeV (v4.9.0) [59]. (3) Enzymes involved in carbohydrates utilization were annotated using the carbohydrate-active enzyme database (CAZy) annotation web server dbCAN (dbCAN-fam-HMMs.txt.v4) based on hmmscan. hmmscan-parser.sh was used to filter output file with default parameters (if alignment > 80aa, use E-value < 1e−5, otherwise use E-value < 1e−3 covered fraction of HMM > 0.3) (hmmer.org) [63] (Additional file 12). Abundance of specific functional classes was determined using hypergeometric analysis and p-values corrected using false discovery rate as described previously [64].

MAGs replication rate

Considering the genome size and the total number of reads mapped on each MAG, the coverage of each MAG was determined using Bowtie 2 (v2.2.4). The MAGs having completeness higher than 90%, contamination lower than 5%, a number of scaffolds per Mbp lower than 175 and a coverage value higher than five, were selected in order to determine their index of replication (iRep) applying the iRep software [45]. Pairwise Wilcoxon rank sum test was performed (pairwise.wilcox.test in R software v3.4.4) and p-values were corrected with Bonferroni adjustment. The number of replication origins in archaeal genomes was inspected using Ori-Finder 2 software [65] and those having none or more than one were excluded from further analyses.

Diversity indices, statistics and PCoA

β-diversity (pairwise sample dissimilarity, clustering method UPGMA) was calculated applying the ExpressBetaDiversity (EBD) software (v1.0.7) [66]. Statistical calculations (Mann–Whitney with Bonferroni correction for identification of taxa enriched in different groups and t-test for the comparison of the number of species in reactors fed with different substrate), diversity indexes (including for example Dominance, Simpson, Shannon H, Evenness, Fisher alpha, Berger–Parker, Chao-1) and β-diversity (pairwise sample dissimilarity, Whittaker) calculations were performed using past software (v3.21) [61]. PCoA was performed with past software using Bray–Curtis as distance measure solely acidogenic reactors were excluded from the analysis due to their strongly different microbial composition.


Watch the video: Uso de KEGG (August 2022).