Genome analysis and pleiotropy assessment using causal networks with loss of function mutation and metabolomics
BMC Genomics volume 20, Article number: 395 (2019)
Many genome-wide association studies have detected genomic regions associated with traits, yet understanding the functional causes of association often remains elusive. Utilizing systems approaches and focusing on intermediate molecular phenotypes might facilitate biologic understanding.
The availability of exome sequencing of two populations of African-Americans and European-Americans from the Atherosclerosis Risk in Communities study allowed us to investigate the effects of annotated loss-of-function (LoF) mutations on 122 serum metabolites. To assess the findings, we built metabolomic causal networks for each population separately and utilized structural equation modeling. We then validated our findings with a set of independent samples. By use of methods based on concepts of Mendelian randomization of genetic variants, we showed that some of the affected metabolites are risk predictors in the causal pathway of disease. For example, LoF mutations in the gene KIAA1755 were identified to elevate the levels of eicosapentaenoate (p-value = 5E-14), an essential fatty acid clinically identified to increase essential hypertension. We showed that this gene is in the pathway to triglycerides, where both triglycerides and essential hypertension are risk factors of metabolomic disorder and heart attack. We also identified that the gene CLDN17, harboring loss-of-function mutations, had pleiotropic actions on metabolites from amino acid and lipid pathways.
Using systems biology approaches for the analysis of metabolomics and genetic data, we integrated several biological processes, which lead to findings that may functionally connect genetic variants with complex diseases.
Lack of knowledge of underlying biological processes in genome wide association studies and disease endpoints has led to a focus on intermediate phenotypes, such as metabolites. Metabolomic profiles are integrated readouts of many biological processes and can functionally connect genetic variants to disease risk factors and complex disease endpoints [1,2,3]. Furthermore, metabolomics can be used to screen for early disease-related changes  and assess effects of external stimuli on living cells . Considering metabolites as an intermediate molecular mediator between genes and clinical endpoints offers potential to illuminate mechanisms underlying a specific single nucleotide polymorphism (SNP)/gene, as well strengthen the association of a gene with the intermediate trait [6, 7]. Thus, metabolite profiles are ideal intermediate biochemical phenotypes for genome wide association studies (GWAS) [6, 8,9,10] and whole genome sequencing studies [11, 12].
Examining metabolomic relationships facilitates understanding of functional links between genetic variants and disease endpoints, which goes beyond traditional association analyses that examine one variable at a time. Analyses that consider the joint effects of multiple traits can have greater statistical power than single-trait analyses of genomic association studies [13,14,15]. In addition, assigning probabilistic models to expression levels of a group of genes (usually within one pathway) enables development of more accurate classifiers and clusters [16,17,18,19]. For high-dimensional metabolomic data, systems approaches based on Mendelian randomization principles can provide an even more comprehensive analysis. To advance a deeper understanding of how genes and metabolites might be interconnected, Yazdani et al. (2016) introduced and applied systems biology approaches that can reveal the underlying relationships and provide insights in metabolomics system [20, 21]. We capitalize on this prior work and introduce an approach to improve genome analysis and assess pleiotropic gene actions.
We investigated the effects of genetic variants on the human metabolome in the populations of white and non-white from the Atherosclerosis Risk in Communities (ARIC) study . In particular, we focus on loss-of-function (LoF) variants that are predicted to result in a non-viable transcript or a greatly truncated protein product . To analyze the LoF-metabolite relationships, we utilized two approaches: 1) a single variant test, and 2) a convex-concave rare variant selection (CCRS) approach [11, 24]. This latter approach selects genetic variants through a penalization linear model, which assumes sparsity in genomic associations while considering the local linkage disequilibrium. We then applied a systems approach called Genome Directed Acyclic Graph (G-DAG)  to model the metabolomic relationships in causal networks. This G-DAG method assured that the necessary assumptions were met in order to then use structural equation modeling to assess genetic findings. This integrative approach facilitates a mechanistic understanding of how genes and metabolites are related by incorporating and modeling relationships of a large number of metabolites and genetic variants, improves genome-metabolite pathway identification, and identifies genes with multiple functions (i.e. pleiotropic actions). Some of the metabolomic causal pathways identified using the G-DAG algorithm matched the availableknowledge [20, 21] and some of the novel findings of the G-DAG algorithm are validated clinically .
To investigate the association of LoF genetic variants with metabolites, we used linear regression to first adjust all metabolites for the covariates age, gender, body mass index, phase (two different time points that the metabolites were measured), and ten principal components (to adjust for population stratification). We then identified mutations in the coding sequence (individually or in aggregate within a gene) with significant effect on metabolites. To achieve this, we used two methods to select strong associations. One approach is a penalized model based on gene-level analyses (simultaneously evaluating all LoF variants in a gene), and is called the CCRS method (Convex-Concave Rare variant Selection). The CCRS method uses a tuning parameter to select LoF-metabolite associations. The set of tuning parameters included 1.5, 0.3 and 0.01, and the tuning parameter that resulted in the minimum BIC was chosen in the final model. The findings are provided in Tables S4 and S5 for EA and AA respectively. The second approach was based on single variant tests, and the level of statistical significance for the test was based on Bonferroni correction for 122 metabolites, 451 and 372 LoF mutations in EA (European-American/white) and AA (African-American/non-white) populations respectively. The findings are provided in Additional file 1: Tables S6 and S7 for EA and AA respectively. We then focused on the variants that were commonly selected by both approaches to reduce the false discovery rate.
The EA and AA participants differed in several ways. The size of covariate effects differed between the EA and AA subjects, the EA and AA individuals were from different geographical regions, and different diets and environment could impact biochemical pathways. Furthermore, ancestry can affect metabolites. In an attempt to reduce confounders, including environmental and regional dietary variations that may impact metabolites, we focused on AA participants from Jackson, Mississippi. We identified the underlying metabolomic relationships separately for the EA and AA groups and we did not expect that the two causal networks EA and AA be identical. However, the two causal networks were very similar with respect to the metabolites with essential roles in the metabolomic system. For example, the relationship between leucine, valine, and isoleucine and the high impact of these three metabolites in the metabolomic system were the same in both causal networks. Among fatty acids, however, the relationships between AA and EA were not very similar. Because differences between EA and AA participants might make it impossible to compare the results between these two populations, we avoided any comparison below.
We integrated the results of LoF-metabolite investigation with the identified metabolomic causal network in each population. To identify underlying relationships among metabolites, we utilized the G-DAG algorithm and identified a metabolomic causal network over 122 metabolites distributed across multiple functional classes in each population. Because the principles of the G-DAG algorithm are based on Mendelian randomization, where genetic variants act as instrumental variables, hence anchoring the direction of causation, the results from this algorithm provide directions of effects in a network. This provided a valid way to subsequently use structural equation modeling, see the method section and . Figure 1 illustrates the steps of the analysis. The findings are presented in Tables 1 and 2 for EA and AA, respectively and some of them are validated using an independent sample.
Genome analysis of AA population
Some of the results from our genome analysis of the AA population are illustrated and described below.
In the genome analysis, GPR97 showed significant effects on the metabolites glycocholenate-sulfate, oleate, and eicoseneate from the lipid pathway; therefore, this gene was hypothesized to have pleiotropic action. Looking at the relationships among the metabolites (Fig. 2), we observed the metabolites oleate and eicoseneate to have direct relationshisps. Using structural equation models to estimate the effect of GPR97 on eicoseneat, we found this effect to not be statistically significant. Therefore, we concluded that GPR97 does not have a pleiotropic action on both metabolites oleate and eicoseneate, but rather GPR97 directly influences only the oleate metabolite.
From the genome analysis, we found that BNIPL has significant effects on the two metabolites octanoylcarnitine and decanoylcarnitine from the lipid pathway (sub-pathway carnitine metabolism). Therefore, we hypothesized that BNIPL has a pleiotropic effect. From the metabolomic causal network, we found that the two metabolites are directly associated (Fig. 3). To assess the hypothesized pleiotropic effect of BNIPL, we modeled the relationship between the two metabolites using structural equation models. The results from this analysis showed that effect of BNIPL on octanoylcarnitine did not remain statistically significant. Therefore, the pleiotropic hypothesis was rejected.
Note that Figs. 2 and 3 both are showing the same causal network. In Fig. 2, we focused on the metabolites affected by GPR97 and in Fig. 3, we focused on the metabolites affected by BNIPL. In each Figure, we highlighted the relationship between the metabolites of interest and brought those metabolites to the border of the causal network to show the relationships.
The impact of LoF mutations on AA and EA metabolomics
Further analyses of the metabolomic causal network identified five modules (densely connected metabolites). We noticed that the identified modules each referred to a set of metabolites with similar function (e.g. metabolites in a pathway) that work together to achieve a coordinated biologic outcome. Below we summarize the effect of LoF mutations on the identified modules.
One of the modules includes hippurate, p-cresol sulfate, methylcatecholsulfate, catecholsulfate, and phenylacetylglutamine, with urea and glutarylcarnitine as neighbors. Investigation of these metabolites indicates that they mostly result from gut metabolism of a diet rich in protein and polyphenols. Many of these metabolites represent major gut contributions to uremic solutes (e.g. methylcatechol sulfate) and uremic toxicity (e.g. p-cresol sulfate) [21, 28]. In the AA causal network, the LoF-metabolite findings revealed that these metabolites were highly influenced by genes harboring LoF mutations: phenylacetylglutamine is influenced by DCLK3 (p = 4e-09), ZSWIM1 (p = 4e-15), and TMPRSS3 (p = 2e-16); urea is influenced by LTK (p = 6e-09), AVEN (p = 2e-08) and CLSPN (p = 1e-13), and glutarylcarnitine is influenced by GUCA1C (p = 4e-11) (see Fig. 4). In contrast, these metabolites were not influenced by any LoF loci in the EA causal network.
In the AA causal network, hormone-related metabolites were significantly influenced by LoF mutations (Fig. 5). These hormone-related metabolites did not influence other metabolites in the metabolomic causal network. Therefore, the LoF loci associated with hormone-related metabolites did not influence other metabolites in the metabolomic causal network.
In the EA causal network, the hormone-related metabolites glycocholenate sulfate, pregn steroid monosulfate, and androsten-3-beta-17-beta-diol-disulfate 1 were influenced by LTK (p = 6e-09), DSE (p = 1e-07), and PLAC4 (p = 4e-09,) respectively.
For the diet-related metabolites heptanoate, adipate, azelate, sebacate, dodecanedioate, and glutarate, the gene CAPN9 (p = 1e-7) influenced heptanoate in the AA population. In the EA causal network, the gene, MYO1A, influenced (p = 4e-09) heptanoate and the gene PDE4DIP influenced dodecanedioate (p = 1e-14). Among 16 long chain fatty acids in the analysis, oleate and nonadecanoate, were affected by LoF loci in the AA causal network. Oleate was affected by two genes, GPR97 (p = 4e-07) and MAP10 (p = 4e-9), and the nonadecanoate was affected by CLEC4C (p = 4e-10). In the EA causal network, long chain fatty acid myristoleate was affected by two genes, ZNF211 (p = 1e-9) and EGFL8 (p = 8e-13).
In the module including amino acid related metabolites in the AA causal network, two peptides from the gamma-glutamyl sub-pathway were influenced by LoF loci: leucine by OR11G2 (p = 5e-10) and CYP2A6 (p = 1e-11), and glutamate by COASY (p = 4e-09) and STPG1 (p = 1e-12). In the EA population, no LoF mutation had any impact on metabolites from the amino acid module.
Findings in EA and AA populations based on metabolite pathways
In AA, 34 metabolites were significantly affected by LoF loci; in EA, 19 metabolites. A summary of the pathway-based metabolites is provided in Table 3.
Figure 6 shows the minor allele frequency (MAF%) for LoF loci and the selected LoF loci in each of the populations. This figure does not show any differences in MAF% between AA and EA populations.
To validate the findings, we used a validation set from the same AA population including 672 samples with 154 LoF carriers with MAF > 7. In the replication analysis, 8 of the gene-metabolite findings were replicated (p < 2.6E-06), depicted with an asterisk in Table 2.
To make fully informed inference about the likely function of genomic variants, systems biology approaches that help to understand the pathway of the genetic effects on disease through intermediate molecular levels offer improvements over standard single-variable analyses commonly used in GWAS. We investigated the effects of LoF genomic variants on the human metabolome in the two ethnic groups of EA and AA participants from the ARIC study. To investigate LoF-metabolite relationships, we selected genetic variants that were found to be statistically associated with metabolites by both a penalized model (CCRS) and single variant tests. We identified metabolomic causal networks using the G-DAG algorithm that is based on a systems biology method. This allowed us to integrate genomic and metabolomic relationships, built on directed acyclic graphs that portray directions, thus meeting the assumptions of the structural equation modeling that we used to assess the LoF-metabolite findings. Because differences between EA and AA participants, such as the size of covariate effects, different geographical regions and as a result different diet, could impact biochemical pathways, we avoided any comparisons of EA and AA metabolomic causal networks in this study.
Although many of our findings were for rare genetic variants, which made it difficult to corroborate our results from published reports, we were able to find published reports that supported some of our results. We identified some LoF mutations with high impact on fatty acid metabolism that may play a role in the pathogenesis of metabolic syndrome and cardiovascular disease. For instance, KIAA1755 has been found to be associated with heart rate from exome chip meta analyses . We found a strong relationship between KIAA1755 and eicosapentaenoate (p < 5e-14). The metabolite eicosapentaenoate is related to essential hypertension, which is the most common type of hypertension with no known cause [30, 31]. From the identified metabolomic causal network based on Mendelian randomization principles, we found the metabolite eicosapentaenoate was among 4 metabolites with high impact on arachidonic acid . It has also been validated clinically that arachidonic acid has the greatest positive impact on triglyceride levels, a risk factor of cardiovascular disease . To assess the effect of triglycerides on essential hypertensive patients, in a separate previous clinical study, 900 patients were examined and a link was found between increased plasma triglyceride levels with more fatal events in essential hypertensive patients . These relationships are depicted in Fig. 7 and may provide some insights into the disease process. The gene KIAA1755 is a gene/protein of unknown function and findings here may reveal new avenues into the gene function and the understanding of the disease etiology.
We found a significant influence of PDE4DIP on the metabolite dodecanedioate from a lipid pathway. Dodecanedioate supplementation in non-insulin-dependent diabetic patients (but not healthy controls) reduced glucose levels without altering insulin levels . Another interesting finding was that the gene MYO1A was found to have a strong impact on the metabolite heptanoate from a lipid pathway. This metabolite has been used as a treatment for Long-chain L-3 hydroxyacyl-CoA dehydrogenase deficiency, a condition in which the body is unable to break down certain fats .
In Table 2, we showed significant associations between CLDN17 with increased levels of isoleucine (pvalue = 2e-12), and APOA1BP (pvalue = 9e-13) with increased levels of leucine, in the AA population. The relationships were through heterozygous LoF variants, suggesting that haploinsufficiency of these genes results in an increase of the corresponding metabolites. Figure 8 illustrates that these two metabolites were highly connected in the causal network, meaning that their effects can spread across the metabolomic system by influencing multiple metabolites directly and indirectly. Note that all participants in AA population were from Jackson, Mississippi, in an attempt to control for environmental confounders, such as regional dietary variations that can influence the metabolome.
Leucine and isoleucine are branched-chain essential amino acids. Studies have shown that higher levels of leucine or isoleucine are strongly related to insulin resistance, obesity, and higher risk of type 2 diabetes [35, 36]. We found a strong association between the gene CLDN17 (Claudin-17) and isoleucine. This gene is a tight junction protein that facilitates anion-selective paracellular transportation and is predominantly expressed in proximal nephrons . There is evidence suggesting the electrogenic transport of amino acids including leucine and isoleucine by proximal tubules . APOA1BP is an apolipoprotein A-I (ApoA-I) binding protein, which is essential in maintaining cholesterol homeostasis and plays a role in cardiovascular disease . Low levels of ApoA-I and HDL cholesterol are linked to risk of cardiovascular events , whereas increases in several amino acids, including leucine, are linked to low levels of HDL cholesterol and insulin resistance in patients with renal dysfunction .
In the AA population, we found the metabolite gammaglutamylleucine from the peptide pathway was influenced by a LoF mutation in CYP2A6 (p-value = 1e-11). The other related findings were relationships between CYP2A13 and the metabolites pyroglutamine (p-value = 5.13e-07) and hydroxyphenyllactate (pvalue = 1.81e-08). The LoF mutation in CYP2A13 revealed a pleiotropic action on two metabolites, pyroglutamine and hydroxyphenyllactate. CYP2A6 and CYP2A13 are both members of cytochrome P450 superfamily of enzymes, which are involved in many catalytic reactions and drug metabolism including nicotine metabolism. Genetic variants in these genes present a lower risk of developing tobacco-related lung cancer in multiethnic populations including African-American smokers [42,43,44]. An increased level of gammaglutamylleucine has been identified as an indicator of anti-obesogenic metabolism and the strongest risk-decreasing predictor of death . Additionally, the inhibition of CYP2A6- and CYP2A13-mediated metabolism of nicotine can increase the activity of gammaglutamylcysteine syntheses that results in the overproduction of pyroglutamate . It has been proposed that pyroglutamate has a xenobiotic detoxifying property . The second metabolite associated with CYP2A13 was hydroxyphenyllactate, which is a tyrosine metabolite, and the L-form of it has been reported to be highly elevated in the urine of patients with pheylketonuria and tyrosinemia. Taken together, using our approach, we identified two separate metabolites gammaglutamylleucine and pyroglutamine that had significant associations with LoF variants of two genes, CYP2A6 and CYP2A13, respectively from the same family of genes. Thus, these metabolites could be functionally related given their protective functions. Furthermore, our analysis method was able to uncover the potentially pleiotropic effects of the gene CYP2A13 on two different metabolites, pyroglutamine and hydroxyphenyllactate, that may act independently of each other. Nonetheless, further study is required to understand the biological functions of these genes on related metabolites and diseases. Identifying pleiotropic genes reveals important biological relationships among molecular components or clinical phenotypes and leads to understanding of complex biological mechanisms, disease pathogenesis and underlying co-morbidities [46, 47]. This understanding may improve and speed drug development and furthermore, predict possible side effects [48, 49].
Although systems approaches to molecular pathway identification will not replace mechanistic experiments, they are complementary and hypothesis-generating, especially in the era of large scale data to narrow the search space, identify targets for intervention, and define pathways that spread the effect of any intervention in the system. The integrated approach introduced here facilitates a mechanistic understanding through incorporating and modeling relationships of a large number of metabolites in genome analysis, improves genome pathway identification, and identifies genes with pleiotropic actions. The focus of the current study was on LoF variants. For future studies, we aim to extend the work to additional previously-identified associated genetic variants  and utilizing new clustering approaches [50, 51].
Metabolomic data function as an intermediate at the molecular level to illuminate mechanisms underlying a specific genetic variant, to identify biological pathways linking the genome to disease, and to discover valuable clinical biomarkers. By integrating results of the genome analysis with metabolomic relationships, and as a result, improving the focus on biological systems, including the identification of pleiotropic gene actions, we improved genome analysis and studied how naturally occurring genetic variants can affect metabolites.
Metabolomic and genomic data
The data was from the ARIC Study, which enrolled 15,792 middle-aged individuals (45 to 64 years) at baseline. For more details of the cohort characteristics see https://www2.cscc.unc.edu/aric/system/files/ CohortCharacteristics.pdf and Additional file 1: Table S1. Serum metabolomic and genomic data were available on 1376 individuals from the AA and EA populations in the study. The dbGAP accession number for ARIC study data is phs000668.v3.p1. A total of 602 metabolites were detected and semi-quantified by Metabolon Inc. (Durham, North Carolina). After excluding metabolites with at least 50% missing values across all samples, metabolites with unknown chemical structures, and metabolites (or any transformation of them) that did not follow a normal distribution, we focused on 122 named metabolites, Additional file 1: Table S8. Details of metabolite assessments are provided in the next section and Supplementary, section Metabolite measurement. Replication samples were available; and we utilized them to identify and impute the missing values with an optimal approach . Since covariates could have a profound effect on metabolites, we adjusted metabolites for age, gender, and body mass index, phase (two different time points that the metabolites were measured), and ten principal components (to adjust for population stratification) using a linear regression. In addition, we carried out the analysis for AA ad EA populations separately. LoF variants were defined as sequence changes caused by single nucleotide variants or small insertions and deletions. The exomes were sequenced on the individuals and 372 and 451 LoF variants were identified for AA and EA populations respectively . Variants were annotated using ANNOVAR  and dbNSFP v2.0  according to the reference genome GRCh37 and National Center for Biotechnology Information Reference Sequence.
After exome capture with VCRome 2.1 (NimbleGen, Inc., Madison, WI), sequencing was carried out using Illumina HiSeq instruments. Using Burrows-Wheeler aligner , sequences were aligned to Genome Reference Consortium Human Build 37. Allele calling and variant-call file construction were carried out with the Atlas2 suite (Atlas-SNP and Atlas-Indel). Single-nucleotide variants (SNVs) were removed using the following criteria: SNVs with a SNP posterior probability less than 0.95, a total depth of coverage less than 6×, an allelic fraction of < 0.1, fewer than three variant reads, 99% reads in a single direction and homozygous reference alleles with < 6× coverage. More stringent filtering was applied on low-quality single-nucleotide substitutions with total depth less than 10; low-quality indels with the following differences: allelic fraction < 0.2 for heterozygous variants and < 0.8 for homozygous variants, minimum total depth less than 60, and variant reads smaller than 30.
All LoF varients in this study were confirmed to have a premature stop codon in a coding exon, disruption of an essential splice site, or an indel predicted to disrupt the downstream reading frame .
Serum metabolomic and genomic data were available on an additional subset of individuals from the same population, 672 African-Americans, which were used to validate the findings in AA population.
Metabolic profiling was carried out on fasting serum samples from the baseline examination stored at − 80 °C. Metabolites were measured using untargeted gas chromatography-mass spectrometry and liquid chromatography-mass spectrometry (GC-MS and LC-MS-based quantification protocols by Metabolon, Inc., Durham, NC). Because the laboratory work for this study is expansive and spanned many days, a data normalization step was performed to correct for variation resulting from differences in instrument tuning from day to day. Metabolites were identified by comparison to library entries of purified standards or recurrent unknown entities. Several types of internal controls were analyzed in concert with the experimental samples. Additional file 1: Tables S2 and S3 describe these quality assurance and quality control samples and standards. Further QC consisted of four major components: the LIMS, the data extraction and peak-identification software, data processing tools for QC and compound identification. The proposed methods, although robust, assume normality of the metabolomics and risk factor data. Log transformation may not be the best transformation for all metabolites. Different transformations were assessed, such as square, root square, log etc., and different metabolites were transformed to a normal distribution using different transformations.
To select the genetic variants associated with metabolites, we applied.
1- Convex-Concave Rare variant Selection (CCRS). The CCRS approach  is a penalized model with some constraints on the design matrix and coefficient vector to provide parsimony of selected covariates, similar to lasso penalization. Using this approach, we had a multivariable regression model which included all LOF variants in the model at a time. Therefore, in contrast to single variant tests, the CCRS approach is a multiple regression method that simultaneously selects the most promising genome variants and estimates their effects on metabolites of interest. Applying the CCRS not only increases the power of identification due to avoiding multiple comparison adjustments, it also takes into account local linkage disequilibrium and prevents overestimation.
2- In addition to the CCRS, to identify mutations in the coding sequence with significant effect on metabolites, we sought association with the metabolites using a linear regression model for each variant at a time. For rare variants, we aggregated their effects by summing the number of LoF variants in each gene and then tested the association of this sum with metabolites by linear regression.
3- We then focused on the variants that were commonly selected by the two steps above to reduce the false discovery rate.
Metabolomic causal networks
Metabolomic causal networks represent metabolites as nodes connected by directed edges indicating the relationships among metabolites. A missing link between two metabolites in the causal network means no relationship, and a link between two metabolites represents the relationship after conditioning on the effect of other metabolites in the analysis (conditional analysis). In a causal network, directions represent cause and effect relationships in observational studies identified based on Mendelian principles/instrumental variables, using variation in the system that is free of confounding [55,56,57]. The G-DAG (genome directed acyclic graph) algorithm [21, 25] with clinically validated novel findings  is based on the principle that the genome inherited variation is a causal factor of metabolomic changes and not the other way around. Some of applications of the G-DAG algorithm are [3, 20, 21, 25, 58]. The G-DAG algorithm was employed to identify metabolomic causal networks, see below for the steps of the G-DAG algorithm. To enforce the assumptions of instrumental variable application and identify robust directions among metabolites, the G-DAG algorithm has the following features :
Extracting information from multiple genetic variants using principal component analysis to create strong instrumental variables;
Independent instrumental variables due to application of principal component analysis;
Multiple independent instrumental variables used for each metabolite, to make overall instrumental variables even stronger.
For identification of the AA and EA metabolomic causal networks, in total 353 and 412 instrumental variables were utilized, respectively. These instrumental variables that were generated from whole exome sequencing data represent the global genomic effects on metabolites despite an individual genetic variant effect. Additional file 1: Figure S1 represents the outcomes of the G-DAG algorithm for EA and AA populations.
To obtain the best fit for the networks, we employed structural Hamming distance , a well-established assessment for the quality of fit in networks, e.g. see [60, 61]. The distance is a function of sample size, variables, and neighbors; a smaller distance indicates a better fit. We considered a set of tuning parameters including 0.0005, 0,001, 0.005, 0.01, 0.05. For the setting with alpha 0.001 and 0.005 the average of the distance on 45 replications was minimum. Therefore, we carried out the analysis at statistical significance level 0.001.
The G-DAG algorithm
The G-DAG algorithm  uses principal components (generated from genome variation) as strong instrumental variables to identify a stable network over the metabolites. Note that the primary aim was identification of a causal network over the set of variables of interest (here metabolites) and the genome information is used as a tool to aid in identifying directionality among the metabolites.
First, we start by reducing the number of SNPs by considering the fact that some SNPs are nearly perfectly correlated (> 0.80) with others, so that one SNP can thereby serve as a proxy for many others in the analysis. To determine a proxy, we use hierarchical clustering and a measure of linkage disequilibrium .
Second, the genomic information was summarized using principal component analysis.
Third, the set of principal components responsible for more than 90% of the variation was selected.
Fourth, the tuning parameter in the G-DAG algorithm was set using the structural Hamming distance, which is a function of variables, sample size and neighboring.
Fifth, in a constraint-based algorithm, the G-DAG identified a causal network over the principal components as instrumental variables and metabolites.
The pleiotropic effect of an inherited gene can refer to a single nucleotide polymorphism (SNP), an entire gene, a large segment of the genome containing multiple genes , or regulatory motifs across the genome . We point out the definition of pleiotropy considered in this study as genes with more than one functionality, the ability of a gene to cause distinct phenotypic traits . The word “distinct” should not be interpreted as independent since some dependency can arise from shared environmental influences and direct physiologic and biochemical relationships . Alternatively, in association studies, where one trait at a time is analyzed, a gene associated with more than one trait cannot be identified with pleiotropic action since this relationship might be due to association between the two traits and not different functionality of the gene. We assessed the pleiotropic action through identification of the genes with direct effect on more than one metabolite.
When a study has a large number of traits, such as metabolites, illumination of underlying relationships provides an opportunity to improve the power of genome analysis, and consequently identification of pleiotropic gene actions. Using causal networks established in principles of Mendelian randomization and application of structural equation modeling, we identified pleiotropic actions by removing findings that are due to metabolomic relationships and not the gene’s direct effect on metabolites.
Structural equation modeling
Here, structural equations were used to model genetic variants with pleiotropic effects based on the relationships of genetic variants and metabolites. In the general form, with M representing metabolites and G representing genetic variants, the structural equations are
where m = 1,…, h; j = h + 1,…, l; i = l + 1; and λij ≠ 0 is equivalent with (j ⟶ i) and γim ≠ 0 is equivalent with (m ⟶i). The notation AM(KR) is called casual parameter . It stands for Assignment Mechanism given the Knowledge about Response . To assure that the assumptions of structural equation modeling were met, we first identified the metabolomic causal networks. The causal networks are illustration of Assignment Mechanisms behind observations [67, 68] which can be illuminated by gathering any Knowledge about Response variables. The reason for conditioning the equations above on the causal parameter is to show that the equations are written based on an identified causal network that meet the required assumption. Note that in a metabolomic causal network, each metabolite is a response variable and the causal network represents metabolites that were involved in assigning value to each response variable (a metabolite). For applications of structural equation modeling see [27, 69].
- ARIC study:
Atherosclerosis Risk in Communities
Convex-Concave Rare Variant Selection
Genome Directed Acyclic Graph
Loss of Function
Structural Equation Modeling
Lewis GD, Wei R, Liu E, Yang E, Shi X, Martinovic M, et al. Metabolite profiling of blood from individuals undergoing planned myocardial infarction reveals early markers of myocardial injury. J Clin Invest. 2008;118(10):3503–12.
Blasco H, Nadal-Desbarats L, Pradat PF, Gordon PH, Madji Hounoum B, Patin F, et al. Biomarkers in amyotrophic lateral sclerosis: combining metabolomic and clinical parameters to define disease progression. Eur J Neurol. 2016;23(2):346–53.
Yazdani A, Yazdani A, Samiei A, Boerwinkle E. A causal network analysis in an observational study identifies metabolomics pathways influencing plasma triglyceride levels. J Biomed Inform. 2016;63:337–43.
Miller MJ, Kennedy AD, Eckhart AD, Burrage LC, Wulff JE, Miller LAD, et al. Untargeted metabolomic analysis for the clinical screening of inborn errors of metabolism. J Inherit Metab Dis. 2015;38(6):1029–39.
Shayanfar S, Broumand A, Pillai SD. Acid stress induces differential accumulation of metabolites in Escherichia coli O26:H11. J Appl Microbiol. 2018.
Suhre K. Genetics meets metabolomics: from experiment to systems biology. Vol. 9781461416, Genetics Meets Metabolomics: From Experiment to Systems Biology. 2012. 1–318 p.
Shah SH, Newgard CB. Integrated metabolomics and genomics: systems approaches to biomarkers and mechanisms of cardiovascular disease. Circ Cardiovasc Genet. 2015;8(2):410–9.
Shin S-Y, Fauman EB, Petersen A-K, Krumsiek J, Santos R, Huang J, et al. An atlas of genetic influences on human blood metabolites. Nat Genet. 2014;46(6):543–50.
Rhee EP, Ho JE, Chen M-H, Shen D, Cheng S, Larson MG, et al. A genome-wide association study of the human metabolome in a community-based cohort. Cell Metab. 2013;18(1):130–43.
Kettunen J, Tukiainen T, Sarin A-P, Ortega-Alonso A, Tikkanen E, Lyytikäinen L-P, et al. Genome-wide association study identifies multiple loci influencing human serum metabolite levels. Nat Genet. 2012;44(3):269–76.
Yazdani A, Yazdani A, Liu X, Boerwinkle E. Identification of rare variants in metabolites of the carnitine pathway by whole genome sequencing analysis. Genet Epidemiol. 2016;40(6):486–91.
Yousri NA, Fakhro KA, Robay A, Rodriguez-Flores JL, Mohney RP, Zeriri H, et al. Whole-exome sequencing identifies common and rare variant metabolic QTLs in a middle eastern population. Nat Commun. 2018;9(1).
Schaid DJ, Tong X, Larrabee B, Kennedy RB, Poland GA, Sinnwell JP. Statistical methods for testing genetic Pleiotropy. Genetics. 2016;204(2):483–97.
Yazdani A, Yazdani A, Giráldez RM, Aguilar DSL. A multi-trait approach identified genetic variants including a rare mutation in RGS3 with impact on abnormalities of cardiac structure/function. Nat-Sci Rep. 2019.
Marchini J, Howie B, Myers S, McVean G, Donnelly P. A new multipoint method for genome-wide association studies by imputation of genotypes. Nat Genet. 2007;39(7):906–13.
Broumand A, Esfahani MS, Yoon BJ, Dougherty ER. Discrete optimal Bayesian classification with error-conditioned sequential sampling. Pattern Recogn. 2015.
Broumand A, Dadaneh SZ. Sequential Sampling for Optimal Bayesian Classification of Sequencing Count Data. In: 52nd Asilomar Conf Signals. CA, USA: Syst Comput Pacific Grove; 2018. p. 1357–61.
Yazdani H, Ortiz-Arroyo D, Choros K, Kwasnicka H. Applying bounded fuzzy possibilistic method on critical objects. In: CINTI 2016 - 17th IEEE international symposium on computational intelligence and informatics: proceedings. 2017.
Yazdani H. Bounded fuzzy Possibilistic method. arXive. 2019.
Yazdani A, Yazdani A, Boerwinkle EA. Causal network analysis of the fatty acid metabolome in African-Americans reveals a critical role for Palmitoleate and Margarate. Omi A J Integr Biol. 2016;20(8):480–4.
Yazdani A, Yazdani A, Samiei A, Boerwinkle E. Identification, analysis, and interpretation of a human serum metabolomics causal network in an observational study. J Biomed Inform. 2016;63:337–43.
Investigators TARIC. The atherosclerosis risk in communities (ARIC) study: design and objectives. The ARIC investigators. Am J Epidemiol. 1989;129(4):687–702.
Liu X, Jian X, Boerwinkle E. dbNSFP v2.0: a database of human non-synonymous SNVs and their functional predictions and annotations. Hum Mutat. 2013;34(9).
Yazdani A, Yazdani A, Boerwinkle E. Rare variants analysis using penalization methods for whole genome sequence data. BMC Bioinformatics. 2015;16(1):405.
Yazdani A, Yazdani A, Samiei A, Boerwinkle E. Generating a robust statistical causal structure over 13 cardiovascular disease risk factors using genomics data. J Biomed Inform. 2016;60:114–9.
Yazdani et al. Arachidonic acid as a target for treating hypertriglyceridemia reproduced by a causal network analysis and an intervention study. Metabolomics. 2018.
Yazdani A, Yazdani A, Boerwinkle E. Conceptual aspects of causal networks in an applied context. J Data Mining Genomics Proteomics. 2016;07(02):2–4.
Patel KP, Luo FJG, Plummer NS, Hostetter TH, Meyer TW. The production of p-cresol sulfate and indoxyl sulfate in vegetarians versus omnivores. Clin J Am Soc Nephrol. 2012;7(6):982–8.
van den Berg ME, Warren HR, Cabrera CP, Verweij N, Mifsud B, Haessler J, et al. Discovery of novel heart rate-associated loci using the exome Chip. Hum Mol Genet. 2017.
Wang S, Ma A, Song S, Quan Q, Zhao X, Zheng X. Fasting serum free fatty acid composition, waist/hip ratio and insulin activity in essential hypertensive patients. Hypertens Res. 2008;31(4):623–32.
Miyajima T, Tsujino T, Saito K, Yokoyama M. Effects of eicosapentaenoic acid on blood pressure, cell membrane fatty acids, and intracellular sodium concentration in essential hypertension. Hypertens Res. 2001;24:537–42.
Turak O, Afşar B, Ozcan F, Öksüz F, Mendi MA, Yayla Ç, et al. The role of plasma triglyceride/high-density lipoprotein cholesterol ratio to predict new cardiovascular events in essential hypertensive patients. J Clin Hypertens. 2016;18(8):772–7.
Karall D, Mair G, Albrecht U, Niedermayr K, Karall T, Schobersberger W, Scholl-Bürgi S. Sports in LCHAD deficiency: maximal incremental and endurance exercise tests in a 13-year-old patient with long-chain 3-Hydroxy acyl-CoA dehydrogenase deficiency (LCHADD) and Heptanoate treatment. In JIMD Reports. 2014;17:7–12.
Greco AV, Mingrone G, Capristo E, Benedetti G, De Gaetano A, Gasbarrini G. The metabolic effect of dodecanedioic acid infusion in non-insulin- dependent diabetic patients. Nutrition. 1998.
Newgard CB, An J, Bain JR, Muehlbauer MJ, Stevens RD, Lien LF, et al. A branched-chain amino acid-related metabolic signature that differentiates obese and lean humans and contributes to insulin resistance. Cell Metab. 2009 Apr;9(4):311–26.
Lotta LA, Scott RA, Sharp SJ, Burgess S, Luan J, Tillin T, et al. Genetic predisposition to an impaired metabolism of the branched-chain amino acids and risk of type 2 diabetes: a Mendelian randomisation analysis. PLoS Med. 2016;13(11).
Krug SM, Günzel D, Conrad MP, Rosenthal R, Fromm A, Amasheh S, et al. Claudin-17 forms tight junction channels with distinct anion selectivity. Cell Mol Life Sci. 2012 Aug;69(16):2765–78.
Jørgensen KE, Kragh-Hansen U, Sheikh MI. Transport of leucine, isoleucine and valine by luminal membrane vesicles from rabbit proximal tubule. J Physiol. 1990 Mar;422:41–54.
Metrustry SJ, Karhunen V, Edwards MH, Menni C, Geisendorfer T, Huber A, et al. Metabolomic signatures of low birthweight: pathways to insulin resistance and oxidative stress. PLoS One. 2018;13(3):e019.
von Weymarn LB, Chun JA, Hollenberg PF. Effects of benzyl and phenethyl isothiocyanate on P450s 2A6 and 2A13: potential for chemoprevention in smokers. Carcinogenesis. 2006;27(4):782–90.
Yu B, Heiss G, Alexander D, Grams ME, Boerwinkle E. Associations between the serum metabolome and all-cause mortality among African Americans in the atherosclerosis risk in communities (ARIC) study. Am J Epidemiol. 2016 Apr;183(7):650–6.
Wassenaar CA, Ye Y, Cai Q, Aldrich MC, Knight J, Spitz MR, et al. CYP2A6 reduced activity gene variants confer reduction in lung cancer risk in African American smokers--findings from two independent populations. Carcinogenesis. 2015 Jan;36(1):99–103.
Park SL, Murphy SE, Wilkens LR, Stram DO, Hecht SS, Le Marchand L. Association of CYP2A6 activity with lung cancer incidence in smokers: the multiethnic cohort study. PLoS One. 2017;12(5):e0178435.
Yuan J-M, Nelson HH, Carmella SG, Wang R, Kuriger-Laber J, Jin A, et al. CYP2A6 genetic polymorphisms and biomarkers of tobacco smoke constituents in relation to risk of lung cancer in the Singapore Chinese health study. Carcinogenesis. 2017 Apr;38(4):411–8.
Weymarn V, B L, Chun JA, Hollenberg PF. Effects of benzyl and phenethyl isothiocyanate on P450s 2A6 and 2A13: potential for chemoprevention in smokers. Carcinogenesis. 2006 Apr;27(4):782–90.
Zhernakova A, Van Diemen CC, Wijmenga C. Detecting shared pathogenesis from the shared genetics of immune-related diseases. Nat Rev Genet. 2009;10:43–55.
Lee SH, Yang J, Goddard ME, Visscher PM, Wray NR. Estimation of pleiotropy between complex diseases using single-nucleotide polymorphism-derived genomic relationships and restricted maximum likelihood. Bioinformatics. 2012;28(19):2540–2.
van de Peppel J, Holstege FCP. Multifunctional genes. Mol Syst Biol. 2005;1(1):E1–2.
Pleiotropy WGC. Natural selection, and the evolution of senescence. Evolution (N Y). 1957;11(4):398.
Yazdani H, Choroś K. Comparative analysis of accuracy of fuzzy clustering methods applied for image processing. In: Advances in intelligent systems and computing. 2019.
Yazdani H, Ortiz-Arroyo D, Choroś K, Kwasnicka H. On high dimensional searching spaces and learning methods. In 2017.
Yazdani A, Yazdani A. Using statistical techniques and replication samples for missing value imputation with an application on metabolomics. J Biostat Epidemiol. 2018.
Wang K, Li M, Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res. 2010;(16):38.
Li H, Durbin R. Fast and accurate long-read alignment with burrows-wheeler transform. Bioinformatics. 2010;26(5):589–95.
Yazdani A, Boerwinkle E. Causal inference at the population level. Int J Res Med Sci. 2014;2(4):1368.
Sheehan NA, Didelez V, Burton PR, Tobin MD. Mendelian randomisation and causal inference in observational epidemiology. PLoS Med. 2008;5(8):1205–10.
Dawid AP. Fundamentals of statistical causality. RSS/EPSRC Grad Train Progr. 2007;(279):1–94.
Yazdani A, Yazdani A, Lorenzi PL, Samiei A. Integrated systems approach identifies pathways from the genome to triglycerides through a Metabolomic causal network. arXiv. 2018;(02):08.
Tsamardinos I, Brown LE, Aliferis CF. The max-min hill-climbing Bayesian network structure learning algorithm. Mach Learn. 2006;65(1):31–78.
Norouzi M, Fleet DJDDJ, Salakhutdinov R, Blei DM. Hamming distance metric learning. Adv Neural Inf Process Syst. 2012:1–9.
Schuster P, Fontana W, Stadler PF, Hofacker IL. From sequences to shapes and Back: a case study in RNA secondary structures. Proc R Soc B Biol Sci. 1994;255(1344):279–84.
Yazdani A, Dunson DB. A hybrid Bayesian approach for genome-wide association studies on related individuals. Bioinformatics. 2015;31(24):3890–6.
Darabos C, Harmon SH, Moore JH. Using the bipartite human phenotype network to reveal pleiotropy and epistasis beyond the gene. Pac Symp Biocomput. 2014;19:188–99.
Chesler EJ, Lu L, Shou S, Qu Y, Gu J, Wang J, et al. Complex trait analysis of gene expression uncovers polygenic and pleiotropic networks that modulate nervous system function. Nat Genet. 2005;37(3):233–42.
Stearns FW. One hundred years of pleiotropy: a retrospective. Genetics. 2010;186:767–73.
He X, Zhang J. Toward a molecular understanding of pleiotropy. Genetics. 2006;173(4):1885–91.
Rubin DB. Causal inference using potential outcomes. J Am Stat Assoc. 2005.
Yazdani ABE. Causal inference in the age of decision medicine. J Data Mining Genomics Proteomics. 2015.
Pearl J. Causal inference in statistics: An overview. Stat Surv. 2009.
The corresponding author, Dr. Yazdani, appreciates Dr. Boerwinkle for providing the data and introducing the pleiotropic topic to her. Dr. Yazdani also thanks the staff and participants of the ARIC Study for their important contributions. ARIC study is supported by the U.S. Public Health Service, National Institutes of Health, contract grant number GM065450.
The first author, Dr. Yazdani was supported by a training fellowship from the Keck Center for Interdisciplinary Bioscience Training of the Gulf Coast Consortia (Grant No. RP140113). DJ Schaid was supported by the U.S. Public Health Service, National Institutes of Health, contract (Grant No. GM065450). MR Kosorok was supported by National Cancer Institute grant P01 CA142538. The funding was for the design and writing of the manuscript.
Availability of data and materials
The dbGAP accession number for ARIC study data is phs000668.v3.p1.
The authors declare that they have no competing interests.
Ethics approval and consent to participate
Consent for publication
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figure S1. European-American metabolomic causal network using the G-DAG algorithm. Table S1. Characteristics of ARIC Cohort. Metabolite assessment and Tables S2. and S3. for QC. Tables S4. and S5. LoF-metabolite relationships using the CCRS approach for EA and AA populations respectively. Tables S6. and S7. LoF-metabolite relationships using the single variant test for EA and AA populations respectively. Table S8. List of 122 named metabolites measured in ARIC study. (DOCX 974 kb)
About this article
Cite this article
Yazdani, A., Yazdani, A., Elsea, S.H. et al. Genome analysis and pleiotropy assessment using causal networks with loss of function mutation and metabolomics. BMC Genomics 20, 395 (2019). https://doi.org/10.1186/s12864-019-5772-4