BACKGROUND
gwas have resulted in the identification of novel genetic loci associated with a variety of diseases and clinical phenotypes.
however, a disease or clinical phenotype is the end point of the behaviour of numerous genes and pathways in addition to environmental influences.
this at least partly explains the general observation that the effect size of genetic association with clinical phenotypes is rather small.
spurred by recent technological developments in the field of metabolomics, interest in genome wide association studies with metabolite levels in blood
<cit>  is gathering momentum.
metabolites are intermediate phenotypes, entities that lie between genes and clinical end points
<cit> .
due to their proximity to an enzyme/gene, metabolites may offer greater effect sizes for gwas than clinical phenotypes
<cit> .
moreover, the pathways in which the metabolite plays a role may provide insight into the underlying biological mechanism responsible for the development of the associated disease.
typically, in metabolomics gwas, hundreds of metabolites are tested for genetic association.
however, association of all snps with all measured metabolites comes with considerable multiple testing problems.
recent publications have also shown that testing ratios of metabolites for genetic association results in much larger effect sizes; however this further exacerbates the multiple testing problem which precludes genuine snp-metabolite pairs from reaching genome-wide significance.
several approaches like gene based tests
<cit>  and pathway analysis
<cit>  have been proposed to overcome this limitation of inadequate statistical power in gwas.
all these approaches have been suggested in the context of gwas with clinical phenotypes but genetic association with metabolites presents its own set of unique opportunities and challenges.
herewith, we explore the utility of background knowledge present in metabolic pathway databases to increase the power in identification of metabolite quantitative trait loci .
our approach involves selective testing of snps near genes in pathways supposedly relevant to the metabolite levels, as a way to reduce the multiple testing burden in gwas.
background knowledge pertaining to a metabolite is retrieved through systematic interrogation of metabolic pathway databases which describe biochemical pathways, reactions, and enzymes relevant to human metabolism.
several pathway databases have been created by groups around the world, while the intent of these efforts remains the elucidation of biological mechanism, the databases however, differ quite significantly in their content, size, user accessibility, download formats and most importantly availability and type of web services for machine-enabled interrogation of the database
<cit> .
in this publication, as a proof of principle, we have chosen to focus on two important metabolic pathway databases, kegg
<cit>  and biocyc
<cit> .
kegg is an integrated database resource of seventeen databases which provide system, genomic and chemical information.
the pathway database consists of both metabolic and non-metabolic pathways and is constructed by a team of curators based on information available in the literature.
biocyc is a collection of pathway/genome databases that describe the genome and metabolic pathways of several organisms.
the database that describes human genomes and pathways, humancyc was interrogated in this study.
in our approach, for every metabolite under consideration, genes acting in the vicinity of the metabolite are determined using knowledge present in databases mentioned above.
we thus generate an integrated set of genes that represent entities with influence over the metabolite.
a workflow management system called taverna
<cit>  was used to generate these gene sets and the snps associated with these genes.
the workflows that were designed for this purpose have been submitted to a workflow repository at http://www.myexperiment.org/packs/ <dig> html <cit> .
a previously published metabolomics dataset by illig et al.
2010
<cit>  was analyzed to evaluate the sensitivity of the method in picking true positives and to identify novel snp-metabolite pairs that had hitherto been obscured in the gwa list given the stringent threshold for significance.
in addition to validating a novel bioinformatics workflow analysis tool, we identified a new locus of human metabolic individuality, aldehyde dehydrogenase family <dig> l <dig> .
this locus was found associated with serine/glycine ratios, a metabolic trait that functionally matches the gene function.
candidate genes identified through the analysis of illig et al.
dataset were taken up for replication in a separate study published by demirkan et al.
<cit> .
we report gpam  and cbs  as novel loci associated with phosphatidylcholine moieties.
RESULTS
our approach can be divided into three stages:  generate a non-redundant gene set for every metabolite considered using knowledge in pathway databases like kegg and biocyc applying interrogation schemes as shown in figure
<dig> and outlined below.
for every gene in the set, generate the set of snps within the gene and 50 kb flanking sequences, and create a snp set for each metabolite  match snps generated for a metabolite with the gwas for the same metabolite and store the matches with the p-values reported for the association .
analysis strategy of databases and interrogation schemes
to retrieve a prioritized list of candidate genes associated with metabolite levels, gene sets were generated for each metabolite through the pathway scheme and the reaction scheme [figure
1a and
1b] for the kegg and biocyc databases .
the pathway scheme generates a list of genes that participate in pathways relevant to the synthesis or degradation of the metabolite.
in the reaction scheme, the metabolite is used as a seed node and shells of reactions around the metabolite are explored.
the list of genes that catalyse the reactions are retrieved and form the gene set for the given metabolite.
for every gene set, a corresponding snp set is generated by retrieving snps within the flanking 50 kb of every gene.
in the final step, the snp set for a metabolite is matched with the gwas dataset for the same metabolite.
at this stage, the sensitivity of the method is evaluated and potential novel discoveries are explored.
results for each of three classes of metabolites  are shown in table
<dig>  for example, for glycine, interrogation of the kegg database identified  <dig> and  <dig> genes using the pathway and reaction schemes respectively, whereas the corresponding numbers of genes were  <dig> and  <dig> for the biocyc database.
the union of all the four interrogation schemes results in a gene set consisting of  <dig> genes relevant to glycine metabolism .
for all the three classes of metabolites,  <dig> unique genes were found,  <dig> are common to kegg and biocyc, the number of genes unique to each of the two databases are  <dig> and  <dig> respectively .
1
2
3
the number of genes for each metabolite and the corresponding database:interrogation scheme is shown.
<dig> the size of the union of the gene set obtained from all the four database:interrogation schemes.
<dig> the size of the corresponding snp set.
<dig> the number of tests is the same as the size of the snp set for the amino acids whereas for aggregated entities like the lipids and carnitine the snp set is multiplied by the number of compounds present in that class.
statistical threshold
the number of unique snps generated for each of the metabolites is shown in table
<dig>  for aggregated metabolites like phosphatidylcholines, sphingomyelins and carnitines the size of the unique snp set is multiplied by the number of metabolites that fall within each class to yield the total number of tests.
for example, the size of the unique snp set for carnitine is  <dig> ; this is multiplied by the number of carnitines which is  <dig>  to yield a total number of  <dig>  tests for these compounds, as shown in the last column of table
<dig>  the sum of all snps derived from our set of metabolites is  <dig> , <dig>  the multiple testing threshold for metabolite concentrations using a bonferroni correction at a nominal p-value of  <dig>  is  <dig> e- <dig> .
in contrast, the p-value threshold for significant association of snps with the same metabolite concentrations in the illig et al.
study would be  <dig> e- <dig> .
this represents a reduction of the multiple testing burden by about two orders of magnitude, regardless of the dependency between the snps or metabolites.
it has been demonstrated that gwas of metabolite ratios offer robust statistical associations and point to biological mechanisms related to the interconversion of metabolite pairs.
to investigate the association of snps with metabolite ratios, we generated the union of snp sets for all combinations of metabolites .
in the case of aggregated metabolites like the lipids and carnitines, the union of the snp set is multiplied by the number of compounds that fall within each class.
for example, the union of the snp set for arginine and carnitine is  <dig> , this is multiplied by  <dig> to yield the total number of  <dig>  tests for this group of ratios.
the number of tests for ratios of compounds within classes such as phosphatidylcholines is equal to the size of the unique snp set multiplied by the number of combinations, n*/ <dig>  which in this case would be 92*91/2 =  <dig>  in choosing combinations of ratios, we have assumed that the association p-value for a linear regression model using a metabolite ratio of a/b is equivalent to that computed using it’s reciprocal, b/a.
the evidence for lack of independence of a ratio and its reciprocal is provided by the illig et al.
study where a comparison of associations computed using untransformed and log-scaled ratios did not detect significant differences.
this implies that we may consider the p-values computed using a/b and b/a to be approximately equal.
the sum of the number of tests for all ratios is  <dig> , <dig> as shown in additional file
1: table s <dig>  the multiple testing threshold for the ratios using bonferroni correction at nominal p-value of  <dig>  is  <dig> e- <dig>  this represents a multiple threshold reduction by two orders of magnitude over the genome-wide threshold estimated by illig et al.
which is  <dig> e- <dig>
proof of principle: sensitivity
the sensitivity of the method was evaluated based on its ability to identify the top hits in the previously published illig et al.
genome-wide association study.
the overall sensitivity of the method as well as the interrogation specific breakdown is shown in table
<dig>  for example, for the biocyc pathway scheme the size of the unique gene set generated for all the metabolites is shown to be  <dig>  the number of genes that are among the  <dig> top hits in the illig et al.
study for this database:interrogation scheme is  <dig> which results in a sensitivity measure of  <dig> .
a metabolite specific breakdown of each of these schemes and the genes with a p-value cut-off of 1e- <dig> is shown in additional file
1: table s <dig>
1
2
3
snapshot of the matches between our method and the association data from the illig et al.
<dig> study for each of the database:interrogation scheme.
1corresponds to the unique set of genes generated for all the metabolites for the given database:interrogation scheme.
2corresponds to the top hits in the illig et al.
publication that were present in the gene set for the given database:interrogation scheme.
3sensitivity is a measure of the actual positives that have been captured by our method and is equal to the ratio of the number of top hits identified by the method over the total number of top hits in the illig et al.
publication which is  <dig>
overall, combining the results from the four database:interrogation schemes helped identify  <dig> of the  <dig> top associations  published by illig et al.
novel discovery in the illig et al dataset
analysis of the first stage or the “discovery stage” dataset of  <dig> samples from the illig et al.
dataset yielded several associations with p-values indicative for association, but that did not meet the significance threshold applied by illig et al.
associations with p-value less than 1e- <dig> were evaluated in the combined “replication stage” dataset with  <dig> samples.
analysis of snps in the aldh1l <dig>  gene locus lowered the p-value of association with serine/glycine ratio from  <dig> e- <dig> in the discovery dataset to  <dig> e- <dig> in the combined dataset.
this is well below our threshold of  <dig> e- <dig>  but above the threshold to be applied when considering all associations between snps and metabolite ratios.
furthermore, the original publication did not select this association for replication because of the threshold set in the first stage of the analysis.
this is an example of the method pointing to potential true positives in a genome-wide scan and the association of aldh1l <dig> with the trait is being reported as a novel discovery.
statistical threshold in the replication study
the analysis of the illig et al.
dataset identified several biologically relevant candidate genes with p-values less than 1e- <dig>  a list of  <dig> of these genes associated with phosphatidylcholines and sphingomyelins were investigated in an independent study in the gwas dataset of phospholipids published by demirkan et al.
the number of matches between the two datasets was:  <dig> phosphatidylcholines and  <dig> sphingomyelins.
demirkan et al.
also performed gwas for within class molar proportions for these moieties.
we took these into consideration in addition to the gwas of absolute concentrations.
therefore, the total number of metabolites and proportions investigated in the demirkan et al.
gwas dataset was  <dig>  a principal component analysis based on the method proposed by li et al.
<cit>  was performed on this set of metabolites resulting in  <dig> effectively independent variables.
as we considered  <dig> independent snps in the candidate loci for these metabolites, the statistical threshold, applying bonferroni correction at a nominal p-value of  <dig> , for the replication study was  <dig> e- <dig> .
novel discoveries in the replication study
table
<dig> shows the top hits in the meta-analysis of candidate genes identified in the illig et al.
dataset for replication.
the meta-analysis was performed using stouffer’s z-score based method of combining p-values
<cit> .
since the snps in the loci replicated in the demirkan et al.
dataset had relatively low r <dig> values with the snps reported in the illig et al.
dataset, we could not perform a traditional meta-analysis where strict linkage disequilibrium criteria are applied.
therefore, we combined the lowest p-value per gene and sought additional supporting evidence for potential allelic heterogeneity .
as mentioned earlier, the p-value threshold for the replication study is set at  <dig> e- <dig>  snps in the vicinity of the genes cbs, gpam, adcy <dig>  cnr <dig>  hsd17b <dig>  mboat <dig>  pecr, plcb <dig> and tecr pass this threshold.
1
2
3
cbs
*
gpam
*
top hits from the meta-analysis of candidate genes identified in the illig et al.
study and replicated in the demirkan et al.
dataset.
<dig> ,3p-value of association of the snp with the trait in the illig et al., demirkan et al.
and combined p-value respectively.
*
indicates genes for which further evidence was found.
discussion
genome wide association studies with metabolites as phenotypes have identified several loci that explain human metabolic individuality.
however, the large metabolite panel being tested results in a severe multiple testing burden that precludes genuine snp-metabolite pairs from consideration when they fail to reach the stringent threshold for statistical significance.
our method aims to address this problem by selectively testing genes that operate in reactions and pathways relevant to the metabolite.
the goal is to reduce the severity of the multiple testing burden and identify potential true positives in the list of genome-wide associations.
taverna, a workflow management system was used to generate the snp-metabolite pairs.
we have deposited the workflows at a repository called myexperiment.org, making it easier for the scientific community to interpret, repeat and reproduce the result.
the sensitivity of the method, defined as retrieval of previously identified associations, is high, as evident from the proof of principle study carried out on the genome scan published by illig et al.
replication studies on some of the promising snp-metabolite pairs identified by the method pointed to a novel and statistically significant association at the aldh1l <dig> locus with serine/glycine ratios.
additional replication studies of phosphatidylcholines and sphingomyelins uncovered significant gene-wise associations with cbs, gpam, adcy <dig>  cnr <dig>  hsd17b <dig>  mboat <dig>  pecr, plcb <dig> and tecr.
databases, interrogation schemes and software tool
the pathway databases have technical and conceptual differences
<cit>  that mandate interrogation of multiple databases and integration of the results.
interpretation of these results requires a close coordination between biologists and computer scientists.
workflow management systems in general and taverna [additional file
1: s2] in particular is an example of a software tool that is intuitive enough for the biologist, while at the same time offering the flexibility for exploring the algorithmic aspects for the computer scientist
<cit> .
in using taverna as a software tool and depositing the workflows in the repository myexperiment.org, we have attempted to make the method and the rationale transparent to users, thus facilitating its retrieval, reuse and reproduction by other independent scientists
<cit> .
sensitivity of the method
as a sensitivity measure of our method, we evaluated its ability to pick the top hits in the illig et al.
publication
<cit> .
some 60% of the top associations were identified successfully.
a similar analysis of gwas dataset published by suhre et al.
<cit>  yielded a sensitivity of 54%  .
however,  <dig> of the “misses” in the suhre et al.
dataset were peptide fragments that do not have an entry in the pathway databases, which is a prerequisite for our method to work.
we interpret the high sensitivity of our method in three ways; first it reinforces the rationale that gwas with metabolomic phenotypes provides a functional approach to the study of human genetic variation
<cit> .
in other words, the known function of the associated gene and the biochemical characteristics of the affected metabolite support each other in ways that lends itself to a narrative on the underlying biological mechanism.
second, while the pathway databases have a long way to go in achieving a comprehensive annotation and delineation of biological processes, they, however, are a good resource of information in so far as the top hits in a gwas with metabolomic phenotypes are concerned.
only two out of the  <dig> top hits in the study by illig et al.
were genes with unknown functions , and two others were hitherto uncharacterized solute transporters .
third, a good sensitivity measure is a validation of our method and reflects its comprehensive data collection ability through integration of disparate data sources and utilization of appropriate interrogation strategies.
novel discoveries
our analysis of the gwas dataset of the illig et al.
publication based on the first step of the “discovery design” yielded several interesting associations that had not been reported among the top hits in the publication.
we selected a few of the promising associations for replication in the combined dataset of  <dig> subjects.
one of the genes, aldehyde dehydrogenase family 1 l <dig>  was found associated with the ratio of serine/glycine with a p-value of  <dig> e- <dig> in the combined set of  <dig> subjects.
aldh1l <dig> also known as 10-formyltetrahydrofolate dehydrogenase  catalyzes the nadp+ dependent oxidation of 10-formyltetrahydrofolate to co <dig> and tetrahydrofolate
<cit>  as shown in figure
<dig>  it plays an important role in folate metabolism
<cit> .
among other functions, aldh1l <dig> has been known to deplete cellular 10-formyltetrahydrofolate pool resulting in a loss of de novo purine biosynthesis
<cit> , maintain cellular folate concentrations by regulating the availability of thf
<cit> , but most importantly, it has been shown to compete with the enzyme serine hydroxymethyl transferase  for the polyglutamyltetrahydrofolates
<cit> .
the latter enzyme catalyzes the conversion of serine to glycine as shown in figure
<dig>  it has also been shown that glycine to serine inter-conversion by shmt accounts for approximately 41% of whole body glycine flux inclusive of both mitochondrial and cytoplasmic processes
<cit> .
to further investigate the potential of our approach to uncover novel genetic associations, we extended the analysis to an additional independent gwas dataset
<cit> .
candidate genes identified in the illig et al.
dataset in association with phosphatidylcholines and sphingomyelins were considered for replication in the dataset provided by demirkan et al.
<cit> .
we discuss here two novel findings for which additional evidence was obtained.
snps near glycerol- <dig> phosphate acyltransferase  are associated with pc ae c34: <dig> moieties in the illig et al.
and demirkan et al.
datasets with p-values of  <dig> e- <dig> and  <dig> e- <dig>  respectively, with a meta-analysis p-value of  <dig> e- <dig>  gpam encodes a mitochondrial protein that esterifies the acyl group from acyl-coa to the sn- <dig> position of glycerol-3-phosphate.
it is a rate-limiting enzyme that catalyzes the initial step in the biosynthesis of triacylglycerols and phospholipids
<cit> .
a recent study showed that in breast cancer, gpam expression is strongly correlated with survival rates, clinico-pathological features as well as metabolomic and lipidomic profiles
<cit> .
interestingly, the study identified the metabolite pc c34: <dig> as the most significantly altered metabolite with respect to gpam expression in breast cancer patients.
this suggests that, for this particular example, genetic control is primarily at the level of gene expression, with secondary effects on enzyme levels and metabolic conversion rates.
the example also highlights the potential influence of genetic variation of metabolic pathways on disease.
a large number of genes identified by our method in the context of phospholipids participate in fatty acid metabolism and are therefore likely to affect the levels of groups of phosphatidylcholines and sphingomyelins.
for example, gpam esterifies the acyl group from acyl-acp to the sn- <dig> position of glycerol-3-phosphate, and is therefore relevant to both acyl-acyl and acyl-alkyl moieties.
the lowest p-value of association, at this locus, with a phosphatidylcholine moiety in the illig et al.
study is with pc ae c36: <dig>  while in the demirkan et al.
study it is pc aa c36: <dig>  since both associations make biological sense, future work should incorporate joint modelling of suitable phospholipid moieties to help identify loci that are biologically relevant but fail to reach the statistical threshold in gwas analysis.
we have reported such best case associations for phosphatidylcholines in additional file
1: table s <dig>
snps near cystathionine beta-synthase  are associated with pc ae c40: <dig> moieties in the illig et al.
and demirkan et al.
datasets with p-values of  <dig> e- <dig> and  <dig> e- <dig>  respectively, with a meta-analysis p-value of  <dig> e- <dig>  mutations in cbs cause hyperhomocysteinemia
<cit> , which is marked by elevated levels of homocysteine.
several studies have associated altered phosphatidylcholine biosynthesis with hyperhomocysteinemia/cbs deficiency
<cit> .
in one of the studies
<cit> , phosphatidylcholine levels and the activity of the enzyme lecithin-cholesterol acyltransferase  were significantly lower in cbs deficient mice than in wild type mice.
while there is considerable literature evidence for the role of cbs in phosphatidylcholine metabolism, the stringent p-value threshold obscures this association in the list of gwas results.
the low r <dig> values for significant snps in gpam, cbs and other loci between the illig et al.
and demirkan et al.
datasets could be explained by allelic heterogeneity.
the latter is a phenomenon where multiple alleles from one gene influence a trait.
however, in some cases it may be that the two apparently independent snps are tagging a third snp
<cit> .
this may be the case for the two snps  near cbs which have an r <dig> of  <dig>  and are associated with c40: <dig> phosphatidylcholines in both the datasets.
however, both snps are in ld with cis-eqtls in the region .
this is suggestive of the snps exerting their effect through the expression levels of the cbs enzyme, as was suggested for gpam.
apparent allelic heterogeneity may preclude identification in a standard meta-analysis, but would justify further investigation of independent or dependent signals at loci showing this phenomenon.
challenges and future direction
in general, our effort was directed at exploring the utility of machine-enabled interrogation of metabolic pathway databases in prioritizing snp-metabolite associations in a gwas dataset.
while the method’s sensitivity and ability to make novel discovery are encouraging, considerable progress needs to be made in metabolite disambiguation to achieve a relevant and comprehensive gene set for a given metabolite.
this problem is particularly acute for phospholipids like phosphatidylcholines and sphingomyelins and various forms of the fatty acid transporters of l-carnitine.
for example, the metabolomics technology used in the illig et al.
study differentiated more than  <dig> forms of phosphatidylcholines based on alkyl or acyl bonds and single or double bonds on the side chains.
however, the pathway databases do not yet contain information for the complex structures.
this forces users to analyze these metabolites at a higher aggregation level.
another issue that requires attention is the bias introduced in selecting genes for inclusion in the gene set.
we have formulated simple rules for interrogation [additional file
1: s1] that facilitates unbiased generation of gene sets for any given metabolite.
another challenge arises due to the high correlation between metabolites, particularly the phospholipids like phosphatidylcholines and sphingomyelins.
these moieties are associated with loci relevant to fatty acid metabolism.
while the variation at these loci effects the levels of fatty acids and thereby the phospholipid pool, to a large extent, these loci are not specific for any particular phospholipid moiety.
as a result, several loci exhibit a pleiotropic effect for biologically related metabolic phenotypes in general and phospholipids in particular [shown in additional file
1: table s7] we have demonstrated that background knowledge and evidence-based approach is ideally suited to identify such candidate genes, however future work should focus on statistical methodologies with sufficient power to detect such pleiotropic loci in gwas of intermediate phenotypes.
in summary, future work includes integration of more pathway databases, metabolite disambiguation, consideration of allelic heterogeneity and multivariate statistical techniques that take into account the high degree of correlation between the metabolites.
CONCLUSIONS
a measurement of metabolites as intermediate phenotypes is a potentially powerful approach to uncover the influence of genetic variation on disease susceptibility and progression.
however, we still face many hurdles in the interpretation of gwas data.
in this study, we investigated the utility of background knowledge present in pathway databases in extending our understanding of the genetic basis of metabolic variation.
we developed a bioinformatics method that prioritizes snp-metabolite associations in a gwas based on metabolic pathway information present in the kegg and biocyc databases.
the validity of the method is demonstrated by re-analysing published gwas datasets and identifying previously known associations.
we report a new locus of human metabolic individuality, aldh1l <dig>  associated with serine/glycine ratios.
replication studies in an independent gwas of phospholipids identified gpam  and cbs  as novel loci, and this was further supported by additional literature evidence.
the utility of a workflow management system in facilitating novel biological discoveries and as a tool for efficient sharing of computational protocols is demonstrated.

 gwas have resulted in the identification of novel genetic loci associated with a variety of diseases and clinical phenotypes. however, a disease or clinical phenotype is the end point of the behaviour of numerous genes and pathways in addition to environmental influences. this at least partly explains the general observation that the effect size of genetic association with clinical phenotypes is rather small. spurred by recent technological developments in the field of metabolomics, interest in genome wide association studies with metabolite levels in blood is gathering momentum. metabolites are intermediate phenotypes, entities that lie between genes and clinical end points. due to their proximity to an enzyme/gene, metabolites may offer greater effect sizes for gwas than clinical phenotypes. moreover, the pathways in which the metabolite plays a role may provide insight into the underlying biological mechanism responsible for the development of the associated disease. typically, in metabolomics gwas, hundreds of metabolites are tested for genetic association. however, association of all snps with all measured metabolites comes with considerable multiple testing problems. recent publications have also shown that testing ratios of metabolites for genetic association results in much larger effect sizes; however this further exacerbates the multiple testing problem which precludes genuine snp-metabolite pairs from reaching genome-wide significance. several approaches like gene based tests and pathway analysis have been proposed to overcome this limitation of inadequate statistical power in gwas. all these approaches have been suggested in the context of gwas with clinical phenotypes but genetic association with metabolites presents its own set of unique opportunities and challenges. herewith, we explore the utility of background knowledge present in metabolic pathway databases to increase the power in identification of metabolite quantitative trait loci. our approach involves selective testing of snps near genes in pathways supposedly relevant to the metabolite levels, as a way to reduce the multiple testing burden in gwas. background knowledge pertaining to a metabolite is retrieved through systematic interrogation of metabolic pathway databases which describe biochemical pathways, reactions, and enzymes relevant to human metabolism. several pathway databases have been created by groups around the world, while the intent of these efforts remains the elucidation of biological mechanism, the databases however, differ quite significantly in their content, size, user accessibility, download formats and most importantly availability and type of web services for machine-enabled interrogation of the database. in this publication, as a proof of principle, we have chosen to focus on two important metabolic pathway databases, kegg and biocyc. kegg is an integrated database resource of seventeen databases which provide system, genomic and chemical information. the pathway database consists of both metabolic and non-metabolic pathways and is constructed by a team of curators based on information available in the literature. biocyc is a collection of pathway/genome databases that describe the genome and metabolic pathways of several organisms. the database that describes human genomes and pathways, humancyc was interrogated in this study. in our approach, for every metabolite under consideration, genes acting in the vicinity of the metabolite are determined using knowledge present in databases mentioned above. we thus generate an integrated set of genes that represent entities with influence over the metabolite. a workflow management system called taverna was used to generate these gene sets and the snps associated with these genes. the workflows that were designed for this purpose have been submitted to a workflow repository at http://www.myexperiment.org/packs/ html. a previously published metabolomics dataset by illig et al. 2010 was analyzed to evaluate the sensitivity of the method in picking true positives and to identify novel snp-metabolite pairs that had hitherto been obscured in the gwa list given the stringent threshold for significance. in addition to validating a novel bioinformatics workflow analysis tool, we identified a new locus of human metabolic individuality, aldehyde dehydrogenase family l. this locus was found associated with serine/glycine ratios, a metabolic trait that functionally matches the gene function. candidate genes identified through the analysis of illig et al. dataset were taken up for replication in a separate study published by demirkan et al.. we report gpam and cbs as novel loci associated with phosphatidylcholine moieties.  our approach can be divided into three stages: generate a non-redundant gene set for every metabolite considered using knowledge in pathway databases like kegg and biocyc applying interrogation schemes as shown in figure  and outlined below. for every gene in the set, generate the set of snps within the gene and 50 kb flanking sequences, and create a snp set for each metabolite match snps generated for a metabolite with the gwas for the same metabolite and store the matches with the p-values reported for the association. analysis strategy of databases and interrogation schemes to retrieve a prioritized list of candidate genes associated with metabolite levels, gene sets were generated for each metabolite through the pathway scheme and the reaction scheme [figure  1a and 1b] for the kegg and biocyc databases. the pathway scheme generates a list of genes that participate in pathways relevant to the synthesis or degradation of the metabolite. in the reaction scheme, the metabolite is used as a seed node and shells of reactions around the metabolite are explored. the list of genes that catalyse the reactions are retrieved and form the gene set for the given metabolite. for every gene set, a corresponding snp set is generated by retrieving snps within the flanking 50 kb of every gene. in the final step, the snp set for a metabolite is matched with the gwas dataset for the same metabolite. at this stage, the sensitivity of the method is evaluated and potential novel discoveries are explored. results for each of three classes of metabolites are shown in table  for example, for glycine, interrogation of the kegg database identified and genes using the pathway and reaction schemes respectively, whereas the corresponding numbers of genes were and for the biocyc database. the union of all the four interrogation schemes results in a gene set consisting of genes relevant to glycine metabolism. for all the three classes of metabolites, unique genes were found, are common to kegg and biocyc, the number of genes unique to each of the two databases are and respectively. 1 2 3 the number of genes for each metabolite and the corresponding database:interrogation scheme is shown. the size of the union of the gene set obtained from all the four database:interrogation schemes. the size of the corresponding snp set. the number of tests is the same as the size of the snp set for the amino acids whereas for aggregated entities like the lipids and carnitine the snp set is multiplied by the number of compounds present in that class. statistical threshold the number of unique snps generated for each of the metabolites is shown in table  for aggregated metabolites like phosphatidylcholines, sphingomyelins and carnitines the size of the unique snp set is multiplied by the number of metabolites that fall within each class to yield the total number of tests. for example, the size of the unique snp set for carnitine is; this is multiplied by the number of carnitines which is to yield a total number of tests for these compounds, as shown in the last column of table  the sum of all snps derived from our set of metabolites is, the multiple testing threshold for metabolite concentrations using a bonferroni correction at a nominal p-value of is e-. in contrast, the p-value threshold for significant association of snps with the same metabolite concentrations in the illig et al. study would be e-. this represents a reduction of the multiple testing burden by about two orders of magnitude, regardless of the dependency between the snps or metabolites. it has been demonstrated that gwas of metabolite ratios offer robust statistical associations and point to biological mechanisms related to the interconversion of metabolite pairs. to investigate the association of snps with metabolite ratios, we generated the union of snp sets for all combinations of metabolites. in the case of aggregated metabolites like the lipids and carnitines, the union of the snp set is multiplied by the number of compounds that fall within each class. for example, the union of the snp set for arginine and carnitine is, this is multiplied by to yield the total number of tests for this group of ratios. the number of tests for ratios of compounds within classes such as phosphatidylcholines is equal to the size of the unique snp set multiplied by the number of combinations, n*/ which in this case would be 92*91/2 =  in choosing combinations of ratios, we have assumed that the association p-value for a linear regression model using a metabolite ratio of a/b is equivalent to that computed using it’s reciprocal, b/a. the evidence for lack of independence of a ratio and its reciprocal is provided by the illig et al. study where a comparison of associations computed using untransformed and log-scaled ratios did not detect significant differences. this implies that we may consider the p-values computed using a/b and b/a to be approximately equal. the sum of the number of tests for all ratios is, as shown in additional file 1: table s the multiple testing threshold for the ratios using bonferroni correction at nominal p-value of is e- this represents a multiple threshold reduction by two orders of magnitude over the genome-wide threshold estimated by illig et al. which is e- proof of principle: sensitivity the sensitivity of the method was evaluated based on its ability to identify the top hits in the previously published illig et al. genome-wide association study. the overall sensitivity of the method as well as the interrogation specific breakdown is shown in table  for example, for the biocyc pathway scheme the size of the unique gene set generated for all the metabolites is shown to be the number of genes that are among the top hits in the illig et al. study for this database:interrogation scheme is which results in a sensitivity measure of. a metabolite specific breakdown of each of these schemes and the genes with a p-value cut-off of 1e- is shown in additional file 1: table s 1 2 3 snapshot of the matches between our method and the association data from the illig et al. study for each of the database:interrogation scheme. 1corresponds to the unique set of genes generated for all the metabolites for the given database:interrogation scheme. 2corresponds to the top hits in the illig et al. publication that were present in the gene set for the given database:interrogation scheme. 3sensitivity is a measure of the actual positives that have been captured by our method and is equal to the ratio of the number of top hits identified by the method over the total number of top hits in the illig et al. publication which is overall, combining the results from the four database:interrogation schemes helped identify of the top associations published by illig et al. novel discovery in the illig et al dataset analysis of the first stage or the “discovery stage” dataset of samples from the illig et al. dataset yielded several associations with p-values indicative for association, but that did not meet the significance threshold applied by illig et al. associations with p-value less than 1e- were evaluated in the combined “replication stage” dataset with samples. analysis of snps in the aldh1l gene locus lowered the p-value of association with serine/glycine ratio from e- in the discovery dataset to e- in the combined dataset. this is well below our threshold of e- but above the threshold to be applied when considering all associations between snps and metabolite ratios. furthermore, the original publication did not select this association for replication because of the threshold set in the first stage of the analysis. this is an example of the method pointing to potential true positives in a genome-wide scan and the association of aldh1l with the trait is being reported as a novel discovery. statistical threshold in the replication study the analysis of the illig et al. dataset identified several biologically relevant candidate genes with p-values less than 1e- a list of of these genes associated with phosphatidylcholines and sphingomyelins were investigated in an independent study in the gwas dataset of phospholipids published by demirkan et al. the number of matches between the two datasets was: phosphatidylcholines and sphingomyelins. demirkan et al. also performed gwas for within class molar proportions for these moieties. we took these into consideration in addition to the gwas of absolute concentrations. therefore, the total number of metabolites and proportions investigated in the demirkan et al. gwas dataset was a principal component analysis based on the method proposed by li et al. was performed on this set of metabolites resulting in effectively independent variables. as we considered independent snps in the candidate loci for these metabolites, the statistical threshold, applying bonferroni correction at a nominal p-value of, for the replication study was e-. novel discoveries in the replication study table  shows the top hits in the meta-analysis of candidate genes identified in the illig et al. dataset for replication. the meta-analysis was performed using stouffer’s z-score based method of combining p-values. since the snps in the loci replicated in the demirkan et al. dataset had relatively low r values with the snps reported in the illig et al. dataset, we could not perform a traditional meta-analysis where strict linkage disequilibrium criteria are applied. therefore, we combined the lowest p-value per gene and sought additional supporting evidence for potential allelic heterogeneity. as mentioned earlier, the p-value threshold for the replication study is set at e- snps in the vicinity of the genes cbs, gpam, adcy cnr hsd17b mboat pecr, plcb and tecr pass this threshold. 1 2 3 cbs * gpam * top hits from the meta-analysis of candidate genes identified in the illig et al. study and replicated in the demirkan et al. dataset.,3p-value of association of the snp with the trait in the illig et al., demirkan et al. and combined p-value respectively. * indicates genes for which further evidence was found. discussion genome wide association studies with metabolites as phenotypes have identified several loci that explain human metabolic individuality. however, the large metabolite panel being tested results in a severe multiple testing burden that precludes genuine snp-metabolite pairs from consideration when they fail to reach the stringent threshold for statistical significance. our method aims to address this problem by selectively testing genes that operate in reactions and pathways relevant to the metabolite. the goal is to reduce the severity of the multiple testing burden and identify potential true positives in the list of genome-wide associations. taverna, a workflow management system was used to generate the snp-metabolite pairs. we have deposited the workflows at a repository called myexperiment.org, making it easier for the scientific community to interpret, repeat and reproduce the result. the sensitivity of the method, defined as retrieval of previously identified associations, is high, as evident from the proof of principle study carried out on the genome scan published by illig et al. replication studies on some of the promising snp-metabolite pairs identified by the method pointed to a novel and statistically significant association at the aldh1l locus with serine/glycine ratios. additional replication studies of phosphatidylcholines and sphingomyelins uncovered significant gene-wise associations with cbs, gpam, adcy cnr hsd17b mboat pecr, plcb and tecr. databases, interrogation schemes and software tool the pathway databases have technical and conceptual differences that mandate interrogation of multiple databases and integration of the results. interpretation of these results requires a close coordination between biologists and computer scientists. workflow management systems in general and taverna [additional file 1: s2] in particular is an example of a software tool that is intuitive enough for the biologist, while at the same time offering the flexibility for exploring the algorithmic aspects for the computer scientist. in using taverna as a software tool and depositing the workflows in the repository myexperiment.org, we have attempted to make the method and the rationale transparent to users, thus facilitating its retrieval, reuse and reproduction by other independent scientists. sensitivity of the method as a sensitivity measure of our method, we evaluated its ability to pick the top hits in the illig et al. publication. some 60% of the top associations were identified successfully. a similar analysis of gwas dataset published by suhre et al. yielded a sensitivity of 54%. however, of the “misses” in the suhre et al. dataset were peptide fragments that do not have an entry in the pathway databases, which is a prerequisite for our method to work. we interpret the high sensitivity of our method in three ways; first it reinforces the rationale that gwas with metabolomic phenotypes provides a functional approach to the study of human genetic variation. in other words, the known function of the associated gene and the biochemical characteristics of the affected metabolite support each other in ways that lends itself to a narrative on the underlying biological mechanism. second, while the pathway databases have a long way to go in achieving a comprehensive annotation and delineation of biological processes, they, however, are a good resource of information in so far as the top hits in a gwas with metabolomic phenotypes are concerned. only two out of the top hits in the study by illig et al. were genes with unknown functions, and two others were hitherto uncharacterized solute transporters. third, a good sensitivity measure is a validation of our method and reflects its comprehensive data collection ability through integration of disparate data sources and utilization of appropriate interrogation strategies. novel discoveries our analysis of the gwas dataset of the illig et al. publication based on the first step of the “discovery design” yielded several interesting associations that had not been reported among the top hits in the publication. we selected a few of the promising associations for replication in the combined dataset of subjects. one of the genes, aldehyde dehydrogenase family 1 l was found associated with the ratio of serine/glycine with a p-value of e- in the combined set of subjects. aldh1l also known as 10-formyltetrahydrofolate dehydrogenase catalyzes the nadp+ dependent oxidation of 10-formyltetrahydrofolate to co and tetrahydrofolate as shown in figure  it plays an important role in folate metabolism. among other functions, aldh1l has been known to deplete cellular 10-formyltetrahydrofolate pool resulting in a loss of de novo purine biosynthesis, maintain cellular folate concentrations by regulating the availability of thf, but most importantly, it has been shown to compete with the enzyme serine hydroxymethyl transferase for the polyglutamyltetrahydrofolates. the latter enzyme catalyzes the conversion of serine to glycine as shown in figure  it has also been shown that glycine to serine inter-conversion by shmt accounts for approximately 41% of whole body glycine flux inclusive of both mitochondrial and cytoplasmic processes. to further investigate the potential of our approach to uncover novel genetic associations, we extended the analysis to an additional independent gwas dataset. candidate genes identified in the illig et al. dataset in association with phosphatidylcholines and sphingomyelins were considered for replication in the dataset provided by demirkan et al.. we discuss here two novel findings for which additional evidence was obtained. snps near glycerol- phosphate acyltransferase are associated with pc ae c34: moieties in the illig et al. and demirkan et al. datasets with p-values of e- and e- respectively, with a meta-analysis p-value of e- gpam encodes a mitochondrial protein that esterifies the acyl group from acyl-coa to the sn- position of glycerol-3-phosphate. it is a rate-limiting enzyme that catalyzes the initial step in the biosynthesis of triacylglycerols and phospholipids. a recent study showed that in breast cancer, gpam expression is strongly correlated with survival rates, clinico-pathological features as well as metabolomic and lipidomic profiles. interestingly, the study identified the metabolite pc c34: as the most significantly altered metabolite with respect to gpam expression in breast cancer patients. this suggests that, for this particular example, genetic control is primarily at the level of gene expression, with secondary effects on enzyme levels and metabolic conversion rates. the example also highlights the potential influence of genetic variation of metabolic pathways on disease. a large number of genes identified by our method in the context of phospholipids participate in fatty acid metabolism and are therefore likely to affect the levels of groups of phosphatidylcholines and sphingomyelins. for example, gpam esterifies the acyl group from acyl-acp to the sn- position of glycerol-3-phosphate, and is therefore relevant to both acyl-acyl and acyl-alkyl moieties. the lowest p-value of association, at this locus, with a phosphatidylcholine moiety in the illig et al. study is with pc ae c36: while in the demirkan et al. study it is pc aa c36: since both associations make biological sense, future work should incorporate joint modelling of suitable phospholipid moieties to help identify loci that are biologically relevant but fail to reach the statistical threshold in gwas analysis. we have reported such best case associations for phosphatidylcholines in additional file 1: table s snps near cystathionine beta-synthase are associated with pc ae c40: moieties in the illig et al. and demirkan et al. datasets with p-values of e- and e- respectively, with a meta-analysis p-value of e- mutations in cbs cause hyperhomocysteinemia, which is marked by elevated levels of homocysteine. several studies have associated altered phosphatidylcholine biosynthesis with hyperhomocysteinemia/cbs deficiency. in one of the studies, phosphatidylcholine levels and the activity of the enzyme lecithin-cholesterol acyltransferase were significantly lower in cbs deficient mice than in wild type mice. while there is considerable literature evidence for the role of cbs in phosphatidylcholine metabolism, the stringent p-value threshold obscures this association in the list of gwas results. the low r values for significant snps in gpam, cbs and other loci between the illig et al. and demirkan et al. datasets could be explained by allelic heterogeneity. the latter is a phenomenon where multiple alleles from one gene influence a trait. however, in some cases it may be that the two apparently independent snps are tagging a third snp. this may be the case for the two snps near cbs which have an r of and are associated with c40: phosphatidylcholines in both the datasets. however, both snps are in ld with cis-eqtls in the region. this is suggestive of the snps exerting their effect through the expression levels of the cbs enzyme, as was suggested for gpam. apparent allelic heterogeneity may preclude identification in a standard meta-analysis, but would justify further investigation of independent or dependent signals at loci showing this phenomenon. challenges and future direction in general, our effort was directed at exploring the utility of machine-enabled interrogation of metabolic pathway databases in prioritizing snp-metabolite associations in a gwas dataset. while the method’s sensitivity and ability to make novel discovery are encouraging, considerable progress needs to be made in metabolite disambiguation to achieve a relevant and comprehensive gene set for a given metabolite. this problem is particularly acute for phospholipids like phosphatidylcholines and sphingomyelins and various forms of the fatty acid transporters of l-carnitine. for example, the metabolomics technology used in the illig et al. study differentiated more than forms of phosphatidylcholines based on alkyl or acyl bonds and single or double bonds on the side chains. however, the pathway databases do not yet contain information for the complex structures. this forces users to analyze these metabolites at a higher aggregation level. another issue that requires attention is the bias introduced in selecting genes for inclusion in the gene set. we have formulated simple rules for interrogation [additional file 1: s1] that facilitates unbiased generation of gene sets for any given metabolite. another challenge arises due to the high correlation between metabolites, particularly the phospholipids like phosphatidylcholines and sphingomyelins. these moieties are associated with loci relevant to fatty acid metabolism. while the variation at these loci effects the levels of fatty acids and thereby the phospholipid pool, to a large extent, these loci are not specific for any particular phospholipid moiety. as a result, several loci exhibit a pleiotropic effect for biologically related metabolic phenotypes in general and phospholipids in particular [shown in additional file 1: table s7] we have demonstrated that background knowledge and evidence-based approach is ideally suited to identify such candidate genes, however future work should focus on statistical methodologies with sufficient power to detect such pleiotropic loci in gwas of intermediate phenotypes. in summary, future work includes integration of more pathway databases, metabolite disambiguation, consideration of allelic heterogeneity and multivariate statistical techniques that take into account the high degree of correlation between the metabolites.  a measurement of metabolites as intermediate phenotypes is a potentially powerful approach to uncover the influence of genetic variation on disease susceptibility and progression. however, we still face many hurdles in the interpretation of gwas data. in this study, we investigated the utility of background knowledge present in pathway databases in extending our understanding of the genetic basis of metabolic variation. we developed a bioinformatics method that prioritizes snp-metabolite associations in a gwas based on metabolic pathway information present in the kegg and biocyc databases. the validity of the method is demonstrated by re-analysing published gwas datasets and identifying previously known associations. we report a new locus of human metabolic individuality, aldh1l associated with serine/glycine ratios. replication studies in an independent gwas of phospholipids identified gpam and cbs as novel loci, and this was further supported by additional literature evidence. the utility of a workflow management system in facilitating novel biological discoveries and as a tool for efficient sharing of computational protocols is demonstrated. 
this paper explores the opportunities and challenges in the analysis of gwas of metabolomic phenotypes and provides novel insights into the genetic basis of metabolic variation through the re-analysis of published gwas datasets. we demonstrate the utility of automated exploitation of background knowledge present in pathway databases for the analysis of gwas datasets of metabolomic phenotypes. using a pathway-based workflow, confirmed previously identified hits and identified a new locus of human metabolic individuality, associating aldehyde dehydrogenase family l with serine/glycine ratios in blood. in contrast, snps that associate with metabolite levels generally explain a higher percentage of the genetic variation and demonstrate larger effect sizes. re-analysis of the published gwas dataset from illig et al. replication in an independent gwas dataset of phospholipids identified two novel loci supported by additional literature evidence: gpam and cbs. we report novel loci and potential biochemical mechanisms that contribute to our understanding of the genetic basis of metabolic variation and its relationship to disease development and progression. genome-wide association studies have identified many common single nucleotide polymorphisms that associate with clinical phenotypes, but these snps usually explain just a small part of the heritability and have relatively modest effect sizes. still, the discovery of snps associated with metabolite levels is challenging since testing all metabolites measured in typical metabolomics studies with all snps comes with a severe multiple testing penalty.
 genome-wide association studies have identified many common single nucleotide polymorphisms that associate with clinical phenotypes, but these snps usually explain just a small part of the heritability and have relatively modest effect sizes. in contrast, snps that associate with metabolite levels generally explain a higher percentage of the genetic variation and demonstrate larger effect sizes. still, the discovery of snps associated with metabolite levels is challenging since testing all metabolites measured in typical metabolomics studies with all snps comes with a severe multiple testing penalty. we have developed an automated workflow approach that utilizes prior knowledge of biochemical pathways present in databases like kegg and biocyc to generate a smaller snp set relevant to the metabolite. this paper explores the opportunities and challenges in the analysis of gwas of metabolomic phenotypes and provides novel insights into the genetic basis of metabolic variation through the re-analysis of published gwas datasets.  re-analysis of the published gwas dataset from illig et al. using a pathway-based workflow, confirmed previously identified hits and identified a new locus of human metabolic individuality, associating aldehyde dehydrogenase family l with serine/glycine ratios in blood. replication in an independent gwas dataset of phospholipids identified two novel loci supported by additional literature evidence: gpam and cbs. in addition, the workflow approach provided novel insight into the affected pathways and relevance of some of these gene-metabolite pairs in disease development and progression.  we demonstrate the utility of automated exploitation of background knowledge present in pathway databases for the analysis of gwas datasets of metabolomic phenotypes. we report novel loci and potential biochemical mechanisms that contribute to our understanding of the genetic basis of metabolic variation and its relationship to disease development and progression. genome-wide associationmetabolitegenotype-phenotype prioritizationbioinformaticspathway databases 