- Primary research
- Open Access
Integrative analysis of multi-omics data to detect the underlying molecular mechanisms for obesity in vivo in humans
Human Genomics volume 16, Article number: 15 (2022)
Obesity is a complex, multifactorial condition in which genetic play an important role. Most of the systematic studies currently focuses on individual omics aspect and provide insightful yet limited knowledge about the comprehensive and complex crosstalk between various omics levels.
Subjects and methods
Therefore, we performed a most comprehensive trans-omics study with various omics data from 104 subjects, to identify interactions/networks and particularly causal regulatory relationships within and especially those between omic molecules with the purpose to discover molecular genetic mechanisms underlying obesity etiology in vivo in humans.
By applying differentially analysis, we identified 8 differentially expressed hub genes (DEHGs), 14 differentially methylated regions (DMRs) and 12 differentially accumulated metabolites (DAMs) for obesity individually. By integrating those multi-omics biomarkers using Mendelian Randomization (MR) and network MR analyses, we identified 18 causal pathways with mediation effect. For the 20 biomarkers involved in those 18 pairs, 17 biomarkers were implicated in the pathophysiology of obesity or related diseases.
The integration of trans-omics and MR analyses may provide us a holistic understanding of the underlying functional mechanisms, molecular regulatory information flow and the interactive molecular systems among different omic molecules for obesity risk and other complex diseases/traits.
Obesity is a chronic metabolic disorder mainly characterized by excessive body fat. Body mass index (BMI) is widely used in obesity research and clinical diagnosis to quantify an individual’s tissue mass. Epidemiological studies estimate that the elevated BMI level is a driving force for the similarly rapid increase of cardiovascular diseases, insulin resistance, type 2 diabetes (T2D), and certain types of cancer . Heritability studies have demonstrated a substantial genetic contribution to obesity risk (h2 ~ 40–70%) . Identification of the genetic determinants for BMI could lead to a better understanding of the biological basis of obesity.
Most systematic studies currently are focused on single to two omics measurements such as at DNA, RNA, metabolite levels. Although useful, little knowledge has been obtained about the cross-talks between molecules of various omics levels and the underlying biological networks that drive complex phenotypes. With the advance of emerging high-throughput sequencing technology such as whole genome sequencing (WGS), RNA sequencing (RNA-Seq), reduced-representation bisulfite sequencing (RRBS) and liquid chromatography–mass spectrometry (LC–MS), multi-omics data including genomics, transcriptomics, epigenomics and metabolomics are rapidly generated and accumulated . As a result, more and more researchers are currently working on the integration of comprehensive multi-omics data to discovery new and meaningful biological knowledge [4, 5], but those studies focused on obesity are rare.
Peripheral blood monocytes (PBMs) or whole blood cells are emerging as a potent source of transcriptomic and epigenetic biomarkers of diabetes and obesity and their related metabolic alterations [6,7,8]. Therefore, we will use PBMs as an example cell type for illustration for this integrative multi-omics study for obesity.
In the current study, we intend to perform systematic genetics analysis which integrates genomic, transcriptomic (from PBM), epigenomic (from PBM) and metabolomic (serum) data of BMI—a major index for obesity, to identify potential molecular and genomic factors/mechanisms underlying the pathogenesis of obesity at different omics level and to re-construct functional module networks to discover the potential regulatory patterns for obesity. Additionally, we seek to identify the significant interactions/networks/causal regulatory relationships within and especially those between omics molecules, shedding lights into the in vivo functional mechanisms for obesity etiology in humans.
Subjects and methods
The study was approved by Tulane University (New Orleans, USA) Institutional Review Board and all participants were required to sign informed consent documents before taking part in the study. A total of 104 premenopausal Caucasian females aged 25–40 years were derived from Louisiana Osteoporosis Study (LOS), an ongoing cohort recruitment since 2011 . The inclusion and exclusion criteria were detailed in our previous studies [10, 11]. All the participants completed an interviewer-assisted comprehensive questionnaire to collect their baseline information including demographic characteristics (age, weight, and height) and life factors (smoking, drinking and exercise, etc.). Weight was measured in light indoor clothing using a calibrated balance beam scale, and height was measured using a calibrated stadiometer without shoes. BMI was calculated as weight (kg) divided by height squared (m2). For the ease of differential analysis, the subjects were categorized into normal weight group and overweight/obesity group according to the WHO criteria.
Peripheral blood mononuclear cells (PBMCs) were isolated from fresh blood from each subject using Lymphoprep™ (Axis-Shield, Oslo, Norway), PBMs were then isolated from PBMCs with Dynabeads® Untouched™ Human Monocytes kit (Life Technologies, CA, USA) using a previously established and routinely performed protocol [12, 13]. Then genomic DNA and total RNA were extracted from the freshly isolated PBMs with the AllPrep® DNA/RNA/ miRNA Universal Kit (Qiagen, CA, USA) and kept at -80℃ for further use. After the collection of whole blood sample, the blood was left undisturbed in room temperature for 15–30 min for coagulation. Then we centrifuged the blood at 1000–2000×g for 10 min to remove the clot, the resulting supernatant yielded the needed serum. Following the centrifugation, the serum was immediately transferred into a clean polypropylene tube and stored at − 80 ℃ or lower. To systematically illuminate the underlying functional mechanisms of obesity, WGS, RRBS (PBM), RNA-seq (PBM), and LC–MS (serum) were performed on the DNA, RNA, and metabolites, respectively.
WGS, RNA-seq, DNA methylation and metabolomic analysis
The WGS, RNA-Seq and epigenome-wide DNA methylation profiling (identified by RRBS) were performed by Technology Center for Genomics & Bioinformatics (TCGB) at University of California, Los Angeles (UCLA). Libraries for WGS were prepared with KAPA DNA LTP library preparation kit (KAPA Biosystem) on Biomek FX Laboratory Automation Workstation (Beckman Coulter). Libraries for RNA-Seq were prepared with KAPA RNA Hyper kit with RiboErase (KAPA Biosystem, USA) according to manufacturer’s instructions. The detailed library construction procedures, library concentration and quality measurement, sequencing protocols and epigenomic analysis methods were described in our previous studies [11, 14, 15].
The LC–MS based metabolomics platform developed by Dr. Garrett’s lab in the Southeast Center for Integrated Metabolomics at University of Florida was used to perform the metabolomic analysis of the study. The detailed laboratory protocols and metabolomics analyses were described in our previous studies [11, 14, 16, 17]. The detailed quality control information for the three omics were also described in our previous study .
For RNA-seq data, we first performed data filtering and normalization , then empirical Bayes method  was used to fit linear models to identify the differentially expressed genes (DEGs) between two groups. Genes with adjusted P value ≤ 0.01 and absolute values of log transformed fold changes (|logFC|) ≥ 5 were considered as significant genes. Finally, those significant DEGs were subject to the Multiscale Embedded Gene Co-expression Network Analysis (MEGENA)  to identify functional co-expressed gene modules and DE hub genes (DEHGs) associated with obesity. MEGENA were implemented in four steps: Fast Planar Filtered Network construction (FPFNC), Multi-scale Clustering analysis (MCA), Multiscale Hub Analysis (MHA) and Cluster-Trait Association Analysis (CTA). For detailed information about the steps, please refer to the method of original MEGENA paper. Functional enrichment analysis was performed on those significant DEGs to using GO enrichment analysis (http://geneontology.org/docs/go-enrichment-analysis/).
For the RRBS data, we adopted logistic regression to identify the differentially methylated regions/bases (DMRs) of multiple CpG sites between the two groups (normal weight group vs overweight/obesity group). Logistic P-values were adjusted to FDR-based Q-values using the Sliding Linear Model (SLIM) method . Methylation regions with Q-value < 0.01 and percent methylation difference (PMD) more than 10% were considered as significant DMRs. Functional enrichment analysis was performed on the genes those DMRs annotated to using DAVID Bioinformatics Resources 6.8 (https://david.ncifcrf.gov/), which provides a comprehensive set of functional annotation tools for investigators to understand biological meaning behind large list of genes.
For the LC–MS data, we conducted both partial least squares regression-discriminant analysis (PLS-DA)  and Logistic regression analysis to detect the differentially accumulated metabolites (DAMs) between two groups. Metabolites with variable importance in projection (VIP) more than 1 and logistic P values less than 0.05 were considered as significant DAMs. Functional enrichment analysis was performed on the identified significant DAMs using Metabolites Biological Role (MBROLE) 2.0 (http://csbg.cnb.csic.es/mbrole2) , which has been widely used to perform metabolites functional annotation and metabolite–protein and drug–protein interactions.
Quantitative trait loci (QTL) analysis to generate the datasets for mendelian randomization (MR) analysis
To identify genetic variants underlying the variation of various omics profiles, we performed QTL analysis by “Matrix eQTL” R package  to identify the expression quantitative loci (eQTL), methylation QTL (meQTL) and metabolomic QTL (metaQTL) for multi-omics data individually. By performing the association analysis between genotype data and RNA expression/DNA methylation/metabolite data individually, eQTL/meQTL/metaQTL datasets were generated for further MR analysis. To reduce the computational burden, we only included DEHGs, DMRs with Q-value < 0.01 and PMD larger than 15%, and the significant DAMs for QTL and further MR analysis. Given the moderate sample size (n = 104), we defined the QTLs that achieve P < 1E−5 as significant QTLs.
MR analysis among multi-omics data
To better understand the crosstalk among the multi-omics data, we first performed Spearman correlation analysis among DEHGs, DMRs and DAMs, and heatmap was generated to represent and visualize their correlation patterns.
Then we applied the bi-directional MR approach to DEHGs and DMRs, DMRs and DAMs, DEHGs and DAMs to detect the potential casual pairs among multi-omics data. To detect the putative causal pairs between gene expression and DNA methylation, we analyzed each gene-methylation pair twice, defining the exposure as either gene expression or methylation. Specifically, we first selected eQTLs with P < 1E−5 as IVs, then the effect estimates of these eQTLs on DNA methylation were extract from the meQTL datasets. When target SNPs were not available in the methylation datasets, we used proxy SNPs that were in high LD (r2 > 0.8) with the SNPs of interest.
Standard inverse-variance weighted (IVW), simple median and weighted median  approaches were utilized to assess the effect estimates of gene expression on DNA methylation. We also applied MR-Egger regression  and mendelian randomization pleiotropy residual sum and outlier (MR-PRESSO) analysis  to evaluate the overall horizontal pleiotropy among all the IVs. MR analysis was then performed on DNA methylation to gene expression as in the bi-directional MR analyses. Similar analysis was repeatedly performed on gene expression and metabolites, DNA methylation and metabolites. To prioritize the sets of the results, any test pairs with at least two MR methods showed P < 0.01 was considered as significant potential causal signals, and 0.01 < P < 0.05 was considered as suggestive evidence for potential causal association.
Network MR analysis
We then applied the network MR analysis  to investigate whether there was mediation effect in those causal pathways. Network MR assumes that the causal effect of the exposure (X) on outcome (Y) is partially mediated through mediator (M). Therefore, the total effect of exposure on outcome are composed of direct effect and indirect effect. Genetic association between genetic variables with exposure (IVX), mediator (IVM), and outcome (IVY) could be derived from the linear regression performed previously. The significant difference (P < 0.05) between indirect effect and total effect suggests the existence of mediation effect.
Our omics workflow was demonstrated in Fig. 1. We observed significant difference in BMI and exercise between the two groups (P < 0.05), and the detailed information was shown in Table 1. Therefore, the following analysis for different omics all adjusted for “exercise”. For the RNA-Seq data, by fitting the gene expression data and BMI group into the linear regression model (adjusted for exercise), we identified 214 DEGs (adjP < 0.01, Additional file 2: Table S1). By using MEGENA, three scales groups (S1, S2, and S3) were identified that had similar interaction patterns and shared highly connected hubs across different scales. These genes were clustered into 17 gene modules (Table 2 and Additional file 2: Table S2) and 8 genes were identified as DEHGs. The module subnetwork figures (Additional file 1: Figs. S1–S4) were used to present the DEHGs of the specific module interconnected with obesity related genes. Modules hierarchy plot (Fig. 2) was generated to visualize the module hierarchical structure. These genes were enriched in obesity-related terms such as “Glycolysis”, “T cell activation”, “Blood coagulation” and “Integrin signaling pathway”. The results of GO term enrichment analysis were detailed in Additional file 2: Table S3.
By using the threshold of Q-value < 0.01 and PMD larger than 10%, we identified 95 DMRs (Additional file 2: Table S4), 35 of them were hyper-methylated regions/bases (Additional file 2: Table S5), and 60 were hypo-methylated regions/bases (Additional file 2: Table S6) when compared with the control group. Then those 95 DMRs were annotated to 67 nearest genes according to the distance to transcriptome start site (TSS) (Additional file 2: Table S7). After using a more stringent threshold of Q-value < 0.01 and PMD larger than 15%, there were 14 DMRs remained (Table 3), and they were annotated to 12 nearest genes. These genes were enriched in obesity-related terms such as “Cytoplasm”, and “transcription, DNA-templated”. The results of were detailed in Additional file 2: Table S8.
By performing PLS-DA and logistic regression analysis, we identified 12 DAMs for obesity (Table 4). These metabolites were enriched in obesity-related terms such as “Amino sugar and nucleotide sugar metabolism”, “Amino Sugar Metabolism”, “Insulin signaling pathway” and “Type 2 diabetes mellitus”. The results of MBROLE term enrichment analysis were detailed in Additional file 2: Table S9.
Bi-directional mendelian randomization analysis among multi-omics data
Spearman correlation analysis demonstrates significant correlation patterns among gene expression, DNA methylation and metabolites (Fig. 3). By performing association analysis between gene expression, DNA methylation, metabolites and genotype data separately, we successfully identified 3560 eQTLs (P < 1E−5) for the DEHGs, 734 meQTLs (P < 1E−5) for the DMRs, and 9055 metaQTLs (P < 1E−5) for DAMs.
We identified 7 of the 112 (8 × 14) DEHG-DMR site pairs with predicted causal direction from the bi-directional MR analyses. Within our 7 predicted causal pairs, five predicted DEHG to causally influence DMR (Additional file 2: Table S10 and Additional file 1: Fig. S5) and two predict DMR to causally influence DEHG. As an example, we highlighted one predicted gene → methylation causal pair ANO6 → 6.110721178. Elevated ANO6 expression was significantly associated with increasing 6.110721178 methylation [β = 2.324, 95% CI (0.691, 3.957), PIVW = 0.005], while 6.110721178 methylation was not associated with ANO6 expression (Additional file 2: Table S10).
We found 40 of the 168 (12 × 14) DMR-DAM site pars with predicted causal direction. Within those 40 observed causal pairs, 31 predicted DAMs were significantly associated with DMRs, 8 suggested DAMs have suggestively causal association with DMRs, and one predicted DMR was causally associated with DAM (Additional file 2: Table S11 and Additional file 1: Figs. S6–S11). We exemplified one predicted causal pair from metabolite to methylation here: isobutyrylcarnitine → 6.110721178. Increased isobutyrylcarnitine was significantly causally associated with increasing 6.110721178 methylation [β = 3.25E−06, 95% CI (1.464E−06, 5.032E−06), PIVW = 0.0003], while 6.110721178 methylation did not show any reverse causal association with isobutyrylcarnitine (Additional file 2: Table S11).
We detected 42 of the 96 (8 × 12) DEHG-DAM site pairs with precited causal direction. Within our 42 predicted causal relationships, 12 predicted DEHG to causally influence DAM and 30 predicted DAM to causally influence DEHG (Additional file 2: Table S12 and Additional file 1: Figs. S12–S17). We demonstrated one of the gene and metabolite causal pairs as an illustration of the causal relationships. Increased ANO6 expression was significantly associated decreasing fructose [β = − 13.732, 95% CI (− 23.070, − 4.394), PIVW = 0.004], though there was no causal association from fructose to ANO6 expression (Additional file 2: Table S12).
Network MR analysis result
Our network MR identified 18 causal pairs with mediation effect (Table 5), including eight causal pathways DAMs → DEHGs → DMRs, four DAMs → DMRs → DEHGs causal pathways and six causal pairs pathways DEHGs → DAMs → DMRs. There were 20 biomarkers involved in those pathways. We highlighted one of those causal pathways: isobutyrylcarnitine_ANO6_6.110721178. We first assessed the causal association between metabolite level of isobutyrylcarnitine and expression level of gene ANO6, it turned out higher isobutyrylcarnitine level was causally associated with the decreasing expression of gene ANO6 (P = 0.003), however the reverse did not suggest any causal association (P = 0.447). Then we tested the association status between expression level of gene ANO6 and methylation level of the region 6.110721178, and the results showed that elevated ANO6 expression level was causally associated with increasing methylation of region 6.110721178 (P = 0.005), but the reverse MR showed no association. Additionally, the mediation analysis showed significant difference between indirect and total effect (P = 0.012), which suggests that the causal association between isobutyrylcarnitine and 6.110721178 was partially mediated by gene ANO6.
In the current study, we first identified significant DEGs, DMRs, and DAMs for obesity in single omics individually. Then, by integration of the multi-omics data (DEHGs, DMRs and DAMs) for obesity using MR analysis, we observed significant correlation among gene expression, DNA methylation and metabolites and identified several putative causal pathways among the various omics molecules. Finally, the application of network MR enabled us to detect 18 causal pathways with mediation effect among different omics.
Implications of the identified DEHGs on obesity
For the identified DEHGs for obesity, there were six previously reported genes (UGGT1, ANO6, MPEG1, PTGS1, CLU and IQGAP1) and two novel genes (LUZP6 and PLCB2) for obesity. Genes UGGT1 and ANO6 were previously reported to be associated with BMI in gluteal subcutaneous adipose tissue (GSAT) , UGGT1 was also reported as hip-associated gene in GSAT . Studies in animal model showed significant difference in gene expression of MPEG1 between normal and obese mice . Expression of gene PTGS1 was reported up-regulated in the intestinal mucosa of the obese rats compared to lean rats , and study reported that PTGS1 expression showed a significant positive correlation with BMI  in human subcutaneous tissue. Microarray analysis in female subcutaneous adipocytes found that CLU gene expression was upregulated in obese versus lean patients , and serum levels of gene CLU was elevated during T2D and coronary heart disease . Studies reported that knockdown of IQGAP1 significantly reduced the ability of glucose to stimulate insulin secretion from β-cell . The rest two genes LUZP6 and PLCB2 were not implicated earlier in obesity or related diseases. However, an extracellular transcriptome study in saliva demonstrated that four extracellular RNA biomarkers including LUZP6 had the potential to distinguish high and low insulin resistance participants . Gene PLCB2 was shown to exhibit diagnostic value for hepatocellular carcinoma , PLCB2 also have an important role in the progression of Alzheimer’s disease and enriched in another neurodegenerative disorder Huntington’s disease . Furthermore, studies with established evidence have reported the associations between cognitive dysfunction, insulin resistance, hepatocellular carcinoma and obesity [39, 40].
Implications of the identified DMRs on obesity
For the identified DMRs, we focused on the 12 nearest genes they were annotated to. For those 12 genes, six of them (DDO, SEPT9, TMEM45B, RXRα, ZNF395 and AHRR) were previously reported to be implicated in the obesity or related diseases and the rest six were novel (PACRG, LINC00494, KLHL4, DTX1, VCX3A and VSTM1). Specifically, according to the genotype–phenotype association of dbGap in Harmonize dataset (http://amp.pharm.mssm.edu/Harmonizome/), DDO was one of the genes that were associated with obesity. A GWAS in a cohort of 1,895 females reported that the variation in gene SEPT9 was correlated with BMI . Gene TMEM45B was proved differentially expressed in white adipose tissue between autism mouse model and wild type mouse model . Also, research suggested that the methylation of gene RXRα was a diagnose biomarker for childhood obesity . While the rest genes were not implicated earlier in the obesity or related diseases, previous research suggested their roles in other complex diseases. Genes PACRG and VSTM1 were reported to be involved in the Parkinson  and rheumatoid arthritis , respectively. Studies already showed that overweight and obesity might be potential risk factors for Parkinson disease  and rheumatoid arthritis . SNPs located in intron of DTX1 were implicated in the process of immune response to HBV vaccination. Studies demonstrated LINC00494 (long intergenic non-protein coding RNA 494) was enriched in prognosis-related long non-coding RNAs (lncRNAs) module that were involved in chemokine signaling pathway, acute myeloid leukemia and pathways in cancer . SNP variants in gene KLHL4 were associated with bone density  and HDL cholesterol . VCX3A was reported to be associated with an abnormal neurocognitive phenotype, which plays a role in neurogenesis regulation . Meanwhile, another gene PNPLA4 in the same region was related to human obesity . Although those genes were not directly implicated in obesity, previous studies showed their associations with other complex diseases that may be related to obesity risk [39, 53].
Implications of the identified DAMs on obesity
As for the identified metabolites for obesity, previous studies reported ten of them (Table 5) were related to obesity, two were novel (Indole-3-acetate and N-methyl-D-aspartic acid (NMDA)). Oral treatment by glucosamine (GLC) in high-fat diet-induced obese rats demonstrated GLC’s effect in preventing body weight gains . Studies of the pathways involved in the obesity and metabolic disorders have showed that ursodeoxycholic acid (UDCA) is used for the treatment of diseases related to liver disorders, especially cholestasis, obesity and lipemic frames , and studies in mouse model of diet-induced obesity also illustrated the effective of UDCA in the treatment of obesity by alleviating metabolic dysfunction [56, 57]. Previous studies involving LC–MS revealed that the sphingosine level in the adipose tissue was increased in obese mice compared to lean mice, furthermore, the plasma level of sphingosine was also indicated elevated in obese mice . Established evidence reported that high consumption of beverage rich in fructose was directly associated with the obesity development and its complications [59, 60]. A study performed in Turkish population illustrated that the obesity prevalence in children with phenylketonuria and hyperphenylalaninaemia who received phenylalanine-restricted diet treatment was higher than that in the normal population . Metabolomic profiling in both obese adults and children reported elevated isobutyrylcarnitine level in plasma in obese than that in lean subjects . Metabolite aspartate was also reported as BMI-associated metabolite . N-acetylneuraminate belongs to sialic acids (SAs) and takes up the highest content of them, SAs were widely common in human tissues and fluids, and research showed that increased level of SA was positively associated with coronary artery disease (CAD) . While for metabolite 3-(2-hydroxyphenyl)propanoate, microbiome study in obese children and adults showed that propanoate was one of the metabolites that were associated with obese individuals , meanwhile, randomized control study of human diet revealed that propanoate level in plasma decreased with the decrease of the weight and it could be used as an independent predictor for insulin sensitivity . As for plasmenyl-LysoPE, plasmenyl was highly accumulated in nerve, immune and cardiovascular systems, study reported it has the potential to protect the cells from reactive oxygen damage , and metabolomic profiling in obese males reported LysoPE was one of the DAMs . To our acknowledge, this study for the first time reported the association of indole-3-acetate and NMDA with obesity. Metabolomics analysis demonstrated that indole-3-acetate was reported associated with carotid intima-media thickness, a validated surrogate marker of atherosclerotic vascular disease . Furthermore, indole-3-acetate was reported could attenuate indicators of inflammation in macrophages and cytokine-mediated lipogenesis in hepatocytes . NMDA receptors are responsible for the majority of excitatory synaptic transmission in the central nervous system, which have been implicated as mediators of neuronal damage caused by excess glutamate in multiple neurologic disorders, including stroke, epilepsy, trauma, and neurodegenerative disorders .
Potential causal regulatory relationship between significant omics profiles
In this analysis, by innovatively using the bi-directional MR principle, we identified significant causal pairs between DEHGs and DMRs, DMRs and DAMs, DEHGs and DAMs. For these gene expression-methylation causal pairs, five of them have evidence of gene driving methylation and the rest two the reverse, which demonstrates that we cannot simply assume methylation always drives changes in gene expression in a model, this findings may represent true causal relationships of gene expression on methylation or methylation on gene expression, but are in themselves not proof. Our findings demonstrated the complex relationship among gene expression, methylation and metabolite, highlighting that the different omics data for complex disease do not simply model one driving another. Therefore, it would make more sense to identify the causal pathways among different omics rather than simply focus on the global relationship of different omics, which enable us to integrate data from different omics levels to reveal their interrelation and combined potential functional influence and pathways on the disease processes. The application of QTL analysis and the integration with MR approach enabled us to detect the specific causal pathways among different omics, which will provide us novel insights into etiology and potential mechanisms underlying complex diseases.
In the current study, our Spearman correlation analysis first demonstrated the significant correlation between biomarker pairs of different omics, then the bi-directional MR analysis further assessed their causal association. The results were partially validated by the previous evidence of their associations with obesity or related diseases. For the 20 biomarkers included in the mediation causal pairs, 17 of them were reported related to the development of obesity or related diseases, which were illustrated earlier. The rest three biomarkers, NMDA, indole-3-acetate and 6:163,743,051 (PACRG-AS1), although they were not implicated in the obesity, research showed their significance with other complicated obesity-related diseases. The results demonstrate the feasibility of the application of MR and network MR in multi-omics data integration, which deepen our understanding of the cross-talks between different omics of obesity and provide us novel insights into discovering the genetic flow information in the pathological of obesity.
Strengths and limitations
Our study has several strengths. First, the application of MEGENA  in identifying the gene co-expression network not only helped us to prioritize meaningful and relevant co-expressed gene clusters for obesity and meanwhile reduces the computational burden for further QTL analysis. Secondly, MR analysis has been extensively applied to multiple integration analysis of multi-omics data such as gene expression and phenotypes, methylation and phenotypes, metabolite and phenotypes, gene expression and methylation. However, there was no application in detecting causal relationship among multi-omics data sets from gene expression, methylation, and metabolites simultaneously. Finally, and particularly, the application of network MR enables us to detect the mediation effect among the causal pathways, which provide us novel insights in unraveling the complex network underlying the mechanisms of obesity, and the biomarkers included in those pathways may serve as potential targets for future precision medicine. To our knowledge, this is the first reported study to integrate multi-omics data of obesity from same population using MR and network MR. We successful identified 18 mediation causal pathways among different omics, which demonstrates the feasibility of MR approach and its effectiveness in helping develop mechanistic insight into the etiology of obesity and other complex diseases.
There are several limitations in the current study. First, our sample size is relatively moderate so that our findings may be limited to those molecules and pathways with most significant effects. Second, current study subjects only included Caucasian females, which may not generalizable to male and other ethnicities. Last, further experimental validation is needed to confirm the biological functional of the identified potential causal pairs in this study.
With the increasing availability of multi-omics or multilayer datasets for complex traits or diseases, the integration analysis of those datasets would be more helpful and powerful in solving the underlying mechanisms. By the application of MR in multi-omics datasets, we were able to untangle some of the crosstalks among various omics molecules and the underlying biological networks that drive the obesity and other complex phenotypes.
Availability of data and material
The datasets generated during and/or analyzed during the current study are available from the corresponding author upon reasonable request.
Body mass index
Type 2 diabetes
Genome-wide association studies
Whole genome sequencing
Reduced-representation bisulfite sequencing
Liquid chromatography–mass spectrometry
Peripheral blood monocytes
Louisiana Osteoporosis Study
Peripheral blood mononuclear cells
Differentially expressed genes
Multiscale Embedded Gene Co-expression Network Analysis
Differentially methylated regions/bases
Percent methylation difference
Partial least squares regression-discriminant analysis
Differentially accumulated metabolites
Variable importance in projection
Quantitative trait loc
Mendelian randomization pleiotropy residual sum and outlier
Kopelman P. Health risks associated with overweight and obesity. Obes Rev. 2007;8(Suppl 1):13–7.
Visscher PM, Brown MA, McCarthy MI, Yang J. Five years of GWAS discovery. Am J Hum Genet. 2012;90(1):7–24.
Hasin Y, Seldin M, Lusis A. Multi-omics approaches to disease. Genome Biol. 2017;18(1):83.
Ge S, Wang Y, Song M, Li X, Yu X, Wang H, Wang J, Zeng Q, Wang W. Type 2 diabetes mellitus: integrative analysis of multiomics data for biomarker discovery. OMICS. 2018;22(7):514–23.
Xu C, Zhang JG, Lin D, Zhang L, Shen H, Deng HW. A systemic analysis of transcriptomic and epigenomic data to reveal regulation patterns for complex disease. G3 Bethesda. 2017;7(7):2271–9.
Keustermans GC, Kofink D, Eikendal A, de Jager W, Meerding J, Nuboer R, Waltenberger J, Kraaijeveld AO, Jukema JW, Sels JW, Garssen J, Prakken BJ, Asselbergs FW, Kalkhoven E, Hoefer IE, Pasterkamp G, Schipper HS. Monocyte gene expression in childhood obesity is associated with obesity and complexity of atherosclerosis in adults. Sci Rep. 2017;7(1):16826.
Reynés B, Priego T, Cifre M, Oliver P, Palou A. Peripheral blood cells, a transcriptomic tool in nutrigenomic and obesity studies: current state of the art. Compr Rev Food Sci Food Saf. 2018;17(4):1006–20.
Baccarelli A, Ghosh S. Environmental exposures, epigenetics and cardiovascular disease. Curr Opin Clin Nutr Metab Care. 2012;15(4):323–9.
Jeng C, Zhao LJ, Wu K, Zhou Y, Chen T, Deng HW. Race and socioeconomic effect on sarcopenia and sarcopenic obesity in the Louisiana Osteoporosis Study (LOS). JCSM Clin Rep. 2018;3(2):1–8.
He H, Liu Y, Tian Q, Papasian CJ, Hu T, Deng HW. Relationship of sarcopenia and body composition with osteoporosis. Osteoporos Int. 2016;27(2):473–82.
Qiu C, Yu F, Su K, Zhao Q, Zhang L, Xu C, Hu W, Wang Z, Zhao L, Tian Q, Wang Y, Deng H, Shen H. Multi-omics data integration for identifying osteoporosis biomarkers and their biological interaction and causal mechanisms. Science. 2020;23(2):100847–100847.
Lei SF, Wu S, Li LM, Deng FY, Xiao SM, Jiang C, Chen Y, Jiang H, Yang F, Tan LJ, Sun X, Zhu XZ, Liu MY, Liu YZ, Chen XD, Deng HW. An in vivo genome wide gene expression study of circulating monocytes suggested GBP1, STAT1 and CXCL10 as novel risk genes for the differentiation of peak bone mass. Bone. 2009;44(5):1010–4.
Liu YZ, Dvornyk V, Lu Y, Shen H, Lappe JM, Recker RR, Deng HW. A novel pathophysiological mechanism for osteoporosis suggested by an in vivo gene expression study of circulating monocytes. J Biol Chem. 2005;280(32):29011–6.
Yu F, Qiu C, Xu C, Tian Q, Zhao L-J, Wu L, Deng H-W, Shen H. Mendelian randomization identifies CpG methylation sites with mediation effects for genetic influences on BMD in peripheral blood monocytes. Front Genet. 2020;11:60–60.
Euan J, Rodger PA, Stockwell A, Chatter E.J.J.O. Biomedicine, biotechnology, Technical considerations for reduced representation bisulfite sequencing with multiplexed libraries. 2012.
Zhao Q, Shen H, Su K-J, Zhang J-G, Tian Q, Zhao L-J, Qiu C, Zhang Q, Garrett TJ, Liu J, Deng H-W. Metabolomic profiles associated with bone mineral density in US Caucasian women. Nutr Metab. 2018;15:57–57.
Liu H, Garrett TJ, Su Z, Khoo C, Gu L. UHPLC-Q-Orbitrap-HRMS-based global metabolomics reveal metabolome modifications in plasma of young women after cranberry juice consumption. J Nutr Biochem. 2017;45:67–76.
Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11(3):R25.
Smyth GK. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. 2004;3:1–25.
Song WM, Zhang B. Multiscale embedded gene co-expression network analysis. PLoS Comput Biol. 2015;11(11):e1004574.
Wang HQ, Tuominen LK, Tsai CJ. SLIM: a sliding linear model for estimating the proportion of true null hypotheses in datasets with dependence structures. Bioinformatics. 2011;27(2):225–31.
You YS, Lin CY, Liang HJ, Lee SH, Tsai KS, Chiou JM, Chen YC, Tsao CK, Chen JH. Association between the metabolome and low bone mineral density in Taiwanese women determined by (1)H NMR spectroscopy. J Bone Miner Res. 2014;29(1):212–22.
López-Ibáñez J, Pazos F, Chagoyen M. MBROLE 2.0—functional enrichment of chemical compounds. Nucl Acids Res. 2016;44(W1):W201–4.
Shabalin AA. Matrix eQTL: ultra fast eQTL analysis via large matrix operations. Bioinformatics. 2012;28(10):1353–8.
Bowden J, Smith GD, Haycock PC, Burgess S. Consistent estimation in Mendelian randomization with some invalid instruments using a weighted median estimator. Genet Epidemiol. 2016;40(4):304–14.
Bowden J, Smith GD, Burgess S. Mendelian randomization with invalid instruments: effect estimation and bias detection through Egger regression. Int J Epidemiol. 2015;44(2):512–25.
Verbanck M, Chen C-Y, Neale B, Do R. Detection of widespread horizontal pleiotropy in causal relationships inferred from Mendelian randomization between complex traits and diseases. Nat Genet. 2018;50(5):693–8.
Burgess S, Daniel RM, Butterworth AS, Thompson SG. Network Mendelian randomization: using genetic variants as instrumental variables to investigate mediation in causal pathways. Int J Epidemiol. 2015;44(2):484–95.
Pinnick KE, Nicholson G, Manolopoulos KN, McQuaid SE, Valet P, Frayn KN, Denton N, Min JL, Zondervan KT, Fleckner J, McCarthy MI, Holmes CC, Karpe F. Distinct developmental profile of lower-body adipose tissue defines resistance against obesity-associated metabolic complications. Diabetes. 2014;63(11):3785–97.
Xu P, Werner JU, Milerski S, Hamp CM, Kuzenko T, Jahnert M, Gottmann P, de Roy L, Warnecke D, Abaei A, Palmer A, Huber-Lang M, Durselen L, Rasche V, Schurmann A, Wabitsch M, Knippschild U. Diet-induced obesity affects muscle regeneration after murine blunt muscle Trauma-A broad spectrum analysis. Front Physiol. 2018;9:674.
Plaza-Diaz J, Robles-Sanchez C, Abadia-Molina F, Moron-Calvente V, Saez-Lara MJ, Ruiz-Bravo A, Jimenez-Valera M, Gil A, Gomez-Llorente C, Fontana L. Adamdec1, Ednrb and Ptgs1/Cox1, inflammation genes upregulated in the intestinal mucosa of obese rats, are downregulated by three probiotic strains. Sci Rep. 2017;7(1):1939.
Quinkler M, Bujalska IJ, Tomlinson JW, Smith DM, Stewart PM. Prostaglandin synthesis in adipose tissue from women with simple obesity reveals characteristic differences between omental and subcutaneous fat depots. Exp Clin Endocrinol Diabetes. 2006;114(S1):OR5_27.
Bradley D, Blaszczak A, Yin Z, Liu J, Joseph JJ, Wright V, Anandani K, Needleman B, Noria S, Renton D, Yearsley M, Wong STC, Hsueh WA. Clusterin impairs hepatic insulin sensitivity and adipocyte clusterin associates with cardiometabolic risk. Diabetes Care. 2019;42(3):466–75.
Trougakos IP, Poulakou M, Stathatos M, Chalikia A, Melidonis A, Gonos ES. Serum levels of the senescence biomarker clusterin/apolipoprotein J increase significantly in diabetes type II and during development of coronary heart disease or at myocardial infarction. Exp Gerontol. 2002;37(10–11):1175–87.
Hedman AC, Smith JM, Sacks DB. The biology of IQGAP proteins: beyond the cytoskeleton. EMBO Rep. 2015;16(4):427–46.
Zhang Y, Sun J, Li F, Grogan TR, Vergara JL, Luan Q, Park M-S, Chia D, Elashoff D, Joshipura KJ, Wong DTW. Salivary extracellular RNA biomarkers for insulin resistance detection in hispanics. Diabetes Res Clin Pract. 2017;132:85–94.
Wang X, Huang K, Zeng X, Liu Z, Liao X, Yang C, Yu T, Han C, Zhu G, Qin W, Peng T. Diagnostic and prognostic value of mRNA expression of phospholipase C β family genes in hepatitis B virus-associated hepatocellular carcinoma. Oncol Rep. 2019;41(5):2855–75.
Chang WS, Wang YH, Zhu XT, Wu CJ. Genome-wide profiling of miRNA and mRNA expression in Alzheimer’s disease. Med Sci Monit. 2017;23:2721–31.
Solas M, Milagro FI, Ramírez MJ, Martínez JA. Inflammation and gut-brain axis link obesity to cognitive dysfunction: plausible pharmacological interventions. Curr Opin Pharmacol. 2017;37:87–92.
Yaribeygi H, Maleki M, Sathyapalan T, Jamialahmadi T, Sahebkar A. Obesity and insulin resistance: a review of molecular interactions. Curr Mol Med. 2020.
Croteau-Chonka DC, Marvelle AF, Lange EM, Lee NR, Adair LS, Lange LA, Mohlke KL. Genome-wide association study of anthropometric traits and evidence of interactions with age and study year in Filipino women. Obesity (Silver Spring, Md). 2011;19(5):1019–27.
Liu X, Tamada K, Kishimoto R, Okubo H, Ise S, Ohta H, Ruf S, Nakatani J, Kohno N, Spitz F, Takumi T. Transcriptome profiling of white adipose tissue in a mouse model for 15q duplication syndrome. Genomics Data. 2015;5:394–6.
Godfrey KM, Sheppard A, Gluckman PD, Lillycrop KA, Burdge GC, McLean C, Rodford J, Slater-Jefferies JL, Garratt E, Crozier SR, Emerald BS, Gale CR, Inskip HM, Cooper C, Hanson MA. Epigenetic gene promoter methylation at birth is associated with child’s later adiposity. Diabetes. 2011;60(5):1528–34.
West AB, Lockhart PJ, O’Farell C, Farrer MJ. Identification of a novel gene linked to parkin via a bi-directional promoter. J Mol Biol. 2003;326(1):11–9.
Wang D, Li Y, Liu Y, He Y, Shi G. Expression of VSTM1-v2 is increased in peripheral blood mononuclear cells from patients with rheumatoid arthritis and is correlated with disease activity. PLOS ONE. 2016;11(1):e0146805.
Martin-Jiménez CA, Gaitán-Vaca DM, Echeverria V, González J, Barreto GE. Relationship between obesity, Alzheimer’s disease, and Parkinson’s disease: an astrocentric view. Mol Neurobiol. 2017;54(9):7096–115.
Dar L, Tiosano S, Watad A, Bragazzi NL, Zisman D, Comaneshter D, Cohen A, Amital H. Are obesity and rheumatoid arthritis interrelated? Int J Clin Pract. 2018;72(1):e13045.
Pan JQ, Zhang YQ, Wang JH, Xu P, Wang W. lncRNA co-expression network model for the prognostic analysis of acute myeloid leukemia. Int J Mol Med. 2017;39(3):663–71.
Styrkarsdottir U, Halldorsson BV, Gretarsdottir S, Gudbjartsson DF, Walters GB, Ingvarsson T, Jonsdottir T, Saemundsdottir J, Snorradóttir S, Center JR, Nguyen TV, Alexandersen P, Gulcher JR, Eisman JA, Christiansen C, Sigurdsson G, Kong A, Thorsteinsdottir U, Stefansson K. New sequence variants associated with bone mineral density. Nat Genet. 2009;41(1):15–7.
Kathiresan S, Manning AK, Demissie S, D’Agostino RB, Surti A, Guiducci C, Gianniny L, Burtt NP, Melander O, Orho-Melander M, Arnett DK, Peloso GM, Ordovas JM, Cupples LA. A genome-wide association study for blood lipid phenotypes in the Framingham Heart Study. BMC Med Genet. 2007;8(Suppl 1):S17.
Sajan SA, Jhangiani SN, Muzny DM, Gibbs RA, Lupski JR, Glaze DG, Kaufmann WE, Skinner SA, Annese F, Friez MJ, Lane J, Percy AK, Neul JL. Enrichment of mutations in chromatin regulators in people with Rett syndrome lacking mutations in MECP2. Genet Med. 2017;19(1):13–9.
Preumont A, Rzem R, Vertommen D, Van Schaftingen E. HDHD1, which is often deleted in X-linked ichthyosis, encodes a pseudouridine-5’-phosphatase. Biochem J. 2010;431(2):237–44.
Savvidis C, Tournis S, Dede AD. Obesity and bone metabolism. Hormones. 2018;17(2):205–17.
Huang L, Chen J, Cao P, Pan H, Ding C, Xiao T, Zhang P, Guo J, Su Z. Anti-obese effect of glucosamine and Chitosan oligosaccharide in high-fat diet-induced obese rats. Mar Drugs. 2015;13:2732–56.
Hofmann AF. Pharmacology of ursodeoxycholic acid, an enterohepatic drug. Scand J Gastroenterol Suppl. 1994;204:1–15.
Zhang Y, Zheng X, Huang F, Zhao A, Ge K, Zhao Q, Jia W. Ursodeoxycholic acid alters bile acid and fatty acid profiles in a mouse model of diet-induced obesity. 2019; 10(842).
Mazzella G, Bazzoli F, Festi D, Ronchi M, Aldini R, Roda A, Grigolo B, Simoni P, Villanova N, Roda E. Comparative evaluation of chenodeoxycholic and ursodeoxycholic acids in obese patients. Effects on biliary lipid metabolism during weight maintenance and weight reduction. Gastroenterology. 1991;101(2):490–6.
Samad F, Hester KD, Yang G, Hannun YA, Bielawski J. Altered adipose and plasma sphingolipid metabolism in obesity. A potential mechanism for cardiovascular and metabolic risk. Diabetes. 2006;55(9):2579–87.
Pereira RM, Botezelli JD, da Cruz-Rodrigues KC, Mekary RA, Cintra DE, Pauli JR, da Silva ASR, Ropelle ER, de Moura LP. Fructose consumption in the development of obesity and the effects of different protocols of physical exercise on the hepatic metabolism. Nutrients. 2017;9(4):405.
Rizkalla SW. Health implications of fructose consumption: a review of recent data. Nutr Metab. 2010;7:82–82.
Ozturk Y, Gencpinar P, Erdur B, Tokgöz Y, Isik I, Akin SB. Overweight and obesity in children under phenylalanine restricted diet. Hong Kong J Paediatr. 2018;23:169–72.
Baker PR, Boyle KE, Koves TR, Ilkayeva OR, Muoio DM, Houmard JA, Friedman JE. Metabolomic analysis reveals altered skeletal muscle amino acid and fatty acid handling in obese humans. Obesity (Silver Spring, Md). 2015;23(5):981–8.
Cirulli ET, Guo L, Leon-Swisher C, Shah N, Huang L, Napier LA, Kirkness EF, Spector TD, Caskey CT, Thorens B, Venter JC, Telenti A. Profound perturbation of the Metabolome in obesity is associated with health risk. Cell Metab. 2019;29(2):488–500.
Rastam L, Lindberg G, Folsom AR, Burke GL, Nilsson-Ehle P, Lundblad A. Association between serum sialic acid concentration and carotid atherosclerosis measured by B-mode ultrasound. The ARIC investigators. Atherosclerosis Risk in Communities Study. Int J Epidemiol. 1996;25(5):953–8.
Del Chierico F, Abbatini F, Russo A, Quagliariello A, Reddel S, Capoccia D, Caccamo R, Ginanni Corradini S, Nobili V, De Peppo F, Dallapiccola B, Leonetti F, Silecchia G, Putignani L. Gut microbiota markers in obese adolescent and adult patients: age-dependent differential patterns. Front Microbiol. 2018;9:1210.
Tirosh A, Calay ES, Tuncman G, Claiborn KC, Inouye KE, Eguchi K, Alcala M, Rathaus M, Hollander KS, Ron I, Livne R, Heianza Y, Qi L, Shai I, Garg R, Hotamisligil GS. The short-chain fatty acid propionate increases glucagon and FABP4 production, impairing insulin action in mice and humans. Sci Transl Med. 2019;11(489):eaav0120.
Moser AB, Steinberg SJ, Watkins PA, Moser HW, Ramaswamy K, Siegmund KD, Lee DR, Ely JJ, Ryder OA, Hacia JG. Human and great ape red blood cells differ in plasmalogen levels and composition. Lipids Health Dis. 2011;10:101.
Wang Y, Liu D, Li Y, Guo L, Cui Y, Zhang X, Li E. Metabolomic analysis of serum from obese adults with hyperlipemia by UHPLC-Q-TOF MS/MS. Biomed Chromatogr. 2016;30(1):48–54.
Boyd A, Boccara F, Meynard JL, Ichou F, Bastard JP, Fellahi S, Samri A, Sauce D, Haddour N, Autran B, Cohen A, Girard PM, Capeau J. Serum tryptophan-derived quinolinate and indole-3-acetate are associated with carotid intima-media thickness and its evolution in HIV-infected treated adults. Open Forum Infect Dis. 2019;6(12):516.
Krishnan S, Ding Y, Saedi N, Choi M, Sridharan GV, Sherr DH, Yarmush ML, Alaniz RC, Jayaraman A, Lee K. Gut microbiota-derived tryptophan metabolites modulate inflammatory response in hepatocytes and macrophages. Cell Rep. 2018;23(4):1099–111.
Lapteva L, Nowak M, Yarboro CH, Takada K, Roebuck-Spencer T, Weickert T, Bleiberg J, Rosenstein D, Pao M, Patronas N, Steele S, Manzano M, van der Veen JW, Lipsky PE, Marenco S, Wesley R, Volpe B, Diamond B, Illei GG. Anti-N-methyl-D-aspartate receptor antibodies, cognitive dysfunction, and depression in systemic lupus erythematosus. Arthritis Rheum. 2006;54(8):2505–14.
We appreciate the support from Zhengzhou University in providing necessary support for this collaborative project.
CQS was partially supported by China National Social Science Foundation (Grant No. 20BRK041). HWD was partially supported by Grants from the NIH (R01-AR069055, U19-AG055373, R01-MH104680, R01-AR059781, and P20-GM109036). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Ethics approval and consent to participate
The study was approved by Tulane University (New Orleans, USA) Institutional Review Board and all participants were required to sign informed consent documents before taking part in the study.
Consent for publication
All the participants signed the consent documents for publications.
The authors declare they have no conflict of interest.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Zhang, Q., Meng, XH., Qiu, C. et al. Integrative analysis of multi-omics data to detect the underlying molecular mechanisms for obesity in vivo in humans. Hum Genomics 16, 15 (2022). https://doi.org/10.1186/s40246-022-00388-x
- Multi-omics integration