An integrated multi-omics analysis of sleep-disordered breathing traits across multiple blood cell types

Overall study design and purpose

The overall purpose of the study was to investigate the multi-omics signaling mechanisms underlying SDB traits to better understand possible drivers of morbidity in SDB. The study design and purpose of each analysis component is illustrated in Figure 1. Briefly, panel A demonstrates the set of associations investigated: SDB phenotypes lead to transcriptional changes which in turn lead to metabolic changes; panel B describes the analysis steps taken to study the potential chain of associations, and the goal of each of these steps. To optimize the available sample size and leverage the fact that transcription is, to some extent, genetically determined, we utilized two separate cohorts to identify the biological components associated sleep exposure to metabolomic changes. Figure 2 further illustrates potential causal relationships underlying a set of measures, and the assumptions that we used to interrogate some of them. Thus, we first performed transcriptome-wide association studies for SDB phenotypes in MESA. For each transcript associated with a SDB trait (FDR p-value <0.1), we used genetic data to construct a transcript Polygenic Risk Score (tPRS) to serve as a predictor of the transcript. Next, to reduce false positive associations in subsequent analyses, we constructed these tPRS in the WHI and tested their association and generalization with their transcripts in whole blood. We proceeded with tPRS results that generalized (p-value <0.05), and constructed and tested them latter for association with SDB phenotypes in HCHS/SOL. If a tPRS was associated with the SDB phenotype in HCHS/SOL, it was interpreted as evidence of reverse association, i.e., the transcript may contribute to SDB. We then calculated the association of the tPRS with metabolites in HCHS/SOL. Lastly, we used another analytic step to support the existence of an association chain linking a sleep exposure, a transcript, and a metabolite: we required evidence of association between the sleep exposure and the metabolite in HCHS/SOL (i.e. any of the potential pathways in column A of Figure 2). If the tPRS was associated with the sleep exposure in HCHS/SOL, it lent support to association chains where the transcript affects both SDB and metabolite levels.

• Download figure
• Open in new tab

Figure 1: Overall study design of the reported analysis.

Flow charts illustrating the methodology and purpose of the analysis. The chart portrays the three SDB phenotypes evaluated in the analysis, the three blood tissues with transcript expression measurement, and demonstrates each step of the analysis with associated goals, cohorts, and reasons why step was performed. Penal A: conceptual linking between SDB measures, transcript expression, and metabolites. Panel B: analytic steps supporting the study of the conceptual links. PBMC tPRS analysis sample sizes in MESA correspond to data from two visits (some individuals were used twice, appropriately accounted for by mixed models).

• Download figure
• Open in new tab

Figure 2: Potential association “chains” of SDB traits, transcript expression, and metabolites, addressed in this study.

The figure illustrates association chains, or pathways, potentially linking an SDB trait, a transcript, and a metabolite. Here we assume that an association between a sleep trait and a transcript was detected in MESA and is assumed “known” for follow-up analysis in WHI in HCHS/SOL. Column A demonstrates potential forms of causal associations between the sleep trait and the transcript, including (A.4) the settings where an association exists due to a common cause, e.g. BMI. Our metabolomics analysis may only detect transcript-metabolite associations, i.e. any sleep-metabolite link is via transcript levels. Column B demonstrate a potential conclusion from an association between a tPRS, validated in WHI and used as an instrumental variable (IV) of a transcript, and a sleep trait: if an association is detected, it provides evidence that changes in transcript levels are upstream (a cause of) changes in sleep trait levels. Column C demonstrates potential conclusions from analyses linking tPRS and a sleep trait to a metabolite. A tPRS is used to link a transcript to a metabolite, and an association, if exists, is likely causal. A sleep-metabolite association should exist if the sleep-transcript and tPRS-transcript associations hold, and therefore observing such an association validates the existence of any of A pathways. Further association between the tPRS and the sleep trait narrows down the potential association chains to A.2 or A.3.

Participating studies

As described in Figure 1, our analysis included three studies: MESA, WHI, and HCHS/SOL each contributing to different analytical steps. The three studies are described in the Supplementary Information. In brief, MESA, our primary study used for discovery of SDB-transcript associations, is a longitudinal cohort study (27). The 1^st and 5^th MESA exams took place between 2000-2002, and 2010-2012, respectively, and whole blood was drawn from participants in both exams. For about 1,000 participants, blood was used later for RNA extraction in at least one of the exams. In addition, a sleep study ancillary to MESA occurred shortly after MESA exam 5 during 2010-2013. Sleep study participants underwent single night in-home polysomnography, as previously described (28). The number of individuals with each type of data and at each time point (exam 1 and exam 5) varies. Figure S1 in the Supplementary Information visualizes the data flow and overlaps across the various measures used in this study: whole-genome genotyping, RNA-seq, and sleep. All MESA participants provided written informed consent, and the study was approved by the Institutional Review Boards at The Lundquist Institute (formerly Los Angeles BioMedical Research Institute) at Harbor-UCLA Medical Center, University of Washington, Wake Forest School of Medicine, Northwestern University, University of Minnesota, Columbia University, and Johns Hopkins University.

• Download figure
• Open in new tab

Figure S1: MESA data flow across various measures in the 1^st and 5^th exam, and the sleep ancillary study.

The WHI was here used to identify tPRS that could be confidently used as IVs for their traits. It is a prospective national health study focused on identifying optimal strategies for preventing chronic diseases that are the major causes of death and disability in postmenopausal women (29). 11,071 WHI participants have whole-genome sequencing data via TOPMed, and 1,274 of these participants have RNA-seq measured in venous blood via TOPMed. All WHI participants provided informed consent and the study was approved by the IRB of the Fred Hutchinson Cancer Research Center.

The HCHS/SOL was used to establish association chains that include an SDB trait, a transcript, and a metabolite. It is a longitudinal cohort study of U.S. Hispanics/Latinos (30, 31). The HCHS/SOL baseline exam occurred on 2008-2011, where 16,415 participants were enrolled.

HCHS/SOL individuals who consented further participated in an in-home sleep study, using a validated type 3 home sleep apnea test recording airflow (via nasal pressure), oximetry, position, and snoring. Genetic data were measured and imputed to the TOPMed freeze 5b reference panel, as previously described, for individuals who consented at baseline (32, 33). Metabolomic data were also measured for n=∼4,000 individuals selected at random out of those with genetic data (34). Figure S2 in the Supplementary Information provides the data flow in HCHS/SOL, focusing on individuals with genetic data and wide consent for genetic data sharing. The HCHS/SOL was approved by the institutional review boards (IRBs) at each field center, where all participants gave written informed consent, and by the Non-Biomedical IRB at the University of North Carolina at Chapel Hill, to the HCHS/SOL Data Coordinating Center. All IRBs approving the study are: Non-Biomedical IRB at the University of North Carolina at Chapel Hill. Chapel Hill, NC; Einstein IRB at the Albert Einstein College of Medicine of Yeshiva University. Bronx, NY; IRB at Office for the Protection of Research Subjects (OPRS), University of Illinois at Chicago. Chicago, IL; Human Subject Research Office, University of Miami. Miami, FL; Institutional Review Board of San Diego State University, San Diego, CA.

• Download figure
• Open in new tab

Figure S2: HCHS/SOL data flow across genotyping, metabolomics, and sleep data.

RNA sequencing

For both MESA and WHI, RNA-seq was performed via the Trans-Omics in Precision Medicine (TOPMed) program. In MESA, RNA-seq was generated from three blood cell types: peripheral blood mononuclear cells (PBMCs; ∼n=1,200 measured in blood from visits 1 and 5), and specific components: T-cells, and monocytes (referred to as T-cell and Mono, for both: n=416), measured in blood from visit 5. Samples were sequenced at the Broad Institute and at the North West Genomics Center (NWGC). Both centers used harmonized protocols. RNA samples quality was assessed using RNA Integrity Number (RIN, Agilent Bioanalyzer) prior to shipment to sequencing centers. QC was re-performed at sequencing centers by RIN analysis at the NWGC and by RNA Quality Score analysis (RQS, Caliper) at the Broad Institute. A minimum of 250ng RNA sample was required as input for library construction, performed using the Illumina TruSeq^TM Stranded mRNA Sample Preparation Kit. RNA was sequenced as 2x101bp paired-end reads on the Illumina HiSeq 4000 according to the manufacturer’s protocols. Target coverage was of ≥40M reads. Comprehensive information about the RNA-seq pipeline used for TOPMed can be found in https://github.com/broadinstitute/gtex-pipeline/blob/master/TOPMed_RNAseq_pipeline.md under MESA RNA-seq pilot commit 725a2bc. Here we used gene-level expected counts quantified using RSEM v1.3.0 (35). RNA sequencing for WHI (whole blood) was performed at the Broad Institute using the unified TOPMed protocols. More information about RNA-seq in WHI is provided in the Supplemental Information.

Metabolomics data in HCHS/SOL

Metabolomics profiling using fasting blood samples was conducted at Metabolon (Durham, NC) with Discovery HD4 platform in 2017. Serum metabolites were quantified with untargeted, liquid chromatography-mass spectrometry (LC-MS)-based quantification protocol (36, 37). The platform captured a total of 1,136 metabolites, including 782 known and 354 unknown (unidentified) metabolites. Detailed methodologic information is provided elsewhere (34).

Phenotypic measures of sleep-disordered breathing (SDB)

We used three SDB measures, as measured by overnight sleep studies in MESA and HCHS/SOL (methods above): (1) the Apnea-Hypopnea Index (AHI), defined in MESA as the number of apneas (breathing cessation) and hypopneas (at least 30% reduction of breath volume, accompanied by 3% or higher reduction of oxyhemoglobin saturation) per hour of sleep, and in HCHS/SOL, due to differences in the recording montage compared to MESA, as the number of apnea or hypopnea events with 3% desaturation per hour of sleep; (2) minimum oxyhemoglobin saturation during sleep (MinO2), and (3) average oxyhemoglobin saturation during sleep (AvgO2).

Testing the association between SDB and blood cell-specific transcriptome-wide gene expression

We used the Olivia R package (38) to perform association analyses of gene expression in PBMCs, monocytes, and T-cells with each of the three SDB measures, separately and in a joint analysis in MESA. SDB phenotypes were treated as the exposures. We followed the recommended Olivia pipeline. Briefly, we performed median normalization, and then filtered lowly expressed gene transcripts defined by removing transcripts with proportion of zero higher than 0.5, median value lower than 1, maximum expression range value lower than 5, and maximum expression value lower than 10. Transcript counts were log transformed after counts of zero were replaced with half the minimum of the observed transcript count in the sample. The analyses were adjusted for age, sex, study center, race/ethnic group, and batch variables: plates, shipment batch, and study site. Because BMI is a strong risk factor for SDB and is assumed to be part of the causal chain, we conducted additional analyses adjusting for BMI (BMIadj). We computed empirical p-values to account for the highly skewed distribution of SDB phenotypes, which may lead to false negative associations if ignored. Finally, we accounted for multiple testing by applying False-Discovery Rate (FDR) correction to each of the association analyses using the Benjamini-Hochberg (BH) procedure (39). We carried forward transcript associations with FDR p-value<0.1 for additional analyses and visualized their association with SDB phenotypes via a hierarchically-clustered heatmap.

Transcript polygenic risk scores (tPRS) construction and validation

To develop tPRS, we first performed a genome-wide association study (GWAS) for each SDB-associated transcript using the MESA TOPMed WGS dataset; each GWAS adjusted for age (years), sex, study site, self-reported race/ethnic background, and 11 principal components, and analyses were restricted to genetic variants with a minor allele frequency of at least 0.05 (due to low sample size). For each GWAS, we used the fully-adjusted two-stage procedure for rank-normalizing residuals in association analyses (40) to identify genetic variants associated with transcript expression. For PBMCs, we used transcript measures from the two MESA visits with RNA-seq data to increase power. To do this, we removed related individuals, and used a random effect model that accounted for individuals. Summary statistics from the GWAS for each transcript were used to develop PRS weights for the corresponding transcript. Next, we constructed tPRS in MESA. We applied clump and threshold implemented in PRSice2 v2.3.1.e (41) using clumping parameters R²=0.1, distance of 250Kb, and three p-value thresholds (5x10^{- 8},10^-7, 10^-6). For each transcript, we constructed the three tPRS in WHI. A tPRS with the smallest p-value in association with the transcript in WHI, and also having p-value<0.05/3=0.017, was selected and considered validated. We also computed FDR-adjusted p-values based on all constructed tPRS (3 candidate tPRS per gene across all genes). To test the association of the tPRS with transcript in WHI, we used logistic mixed models, executed with the GENESIS R package (42) version 2.16.1. Each tPRS served as the exposure, and transcripts served as the outcome, here too using the two-stage procedure for rank-normalization (40). Relatedness was modeled via a sparse kinship matrix among TOPMed WHI individuals. We selected transcripts with p-value <0.017 for follow-up analysis.

We validated that our approach to construct tPRS is robust. We compared a few polygenic prediction models developed using bulk RNAseq in monocytes. First, the prediction model developed using prediXcan based on the MESA dataset (43, 44), with weights provided in the predictDB database (http://www.predictdb.org/). Second, our approach above using genome-wide SNPs (including trans-eQTLs), and third, a similar clump and threshold approach as above limited to cis-eQTLs defined as SNPs within 1Mbp of the start and end position of the transcript (the definition used by prediXcan). We focused on monocytes for this comparison because prediXcan models were only published based on monocytes.

Using tPRS to identify reverse association between gene expression and SDB traits

We constructed generalized tPRS in HCHS/SOL. We used HCHS/SOL genotypes imputed to the TOPMed freeze5b reference panel. Prior to tPRS construction, we filtered SNPs with imputation quality <0.8, minor allele frequency <5%, missingness rate >0.01. As illustrated in column B of Figure 2, We identified potential reverse causation, where gene expression alters SDB, by using the tPRS constructed in HCHS/SOL as instrumental variables (IVs) and testing their association with their respective SDB phenotypes in HCHS/SOL. We used logistic mixed models, executed with the GENESIS R package (42) version 2.16.1. Each tPRS served as the exposure, and the relevant SDB phenotype served as the outcome. To account for skewness of the SDB phenotypes, we used the two-stage procedure for rank-normalization (40). Relatedness was modeled via a sparse kinship matrix, household sharing, and block unit sharing among HCHS/SOL individuals. Association analyses were adjusted for age, sex, study site, Hispanic/Latino background, the first 5 PCs of the genetic data, and log of the sampling weights used to sample HCHS/SOL individuals into the study. Because the tPRS represent a genetic proxy for gene expression, if a tPRS was found to be associated with a SDB phenotype (p-value<0.05), it provided evidence that the transcript contributed to the SDB phenotype, rather than vice versa. However, as illustrated in diagrams A.2 and A.3 in Figure 2 for sleep-transcript association, bidirectional associations are also plausible.

Associations between tPRS and metabolites

Treating tPRS as genetic IVs for gene expression, we estimated associations between tPRS and all identified (named) metabolites with < 25% missing values in HCHS/SOL. We used robust survey models implemented in the R survey package version 4.0 (45), accounting for HCHS/SOL study design (probability sampling and clustering) and providing associations generalizable to the HCHS/SOL target population. For each metabolite, we first imputed observations with missing values of that metabolite with its minimum value observed in the sample, under the assumption that missing values are due to concentrations being below the detection limit, and then rank-normalized it across the sample. We used the same covariates as before: age, sex, study site, Hispanic background, and the first 5 PCs of the genetic data.

Furthermore, we adjusted for BMI depending on the original association of the SDB phenotype and the transcript (BMI unadjusted or BMI adjusted). For each transcript, we corrected metabolite associations to account for FDR using the Benjamini-Hochberg (BH) procedure (39). Associations were considered significant if the FDR p-value was <0.05.

Association analyses of SDB traits with selected metabolites to verify a complete association chain

To further validate a complete association “chain” as detailed in Figure 2, we performed association analyses between the SDB phenotypes and metabolites identified in the tPRS analysis. Associations between SDB phenotypes and metabolites used a survey sampling approach to account for HCHS/SOL sampling design and obtained estimates generalizable to the HCHS/SOL target population. Thus, we used the survey R package (46) with each individual weighted by their sampling weights, and clustering accounted for when computing robust standard errors. Analyses were adjusted for age, sex, study site, Hispanic/Latino background. and BMI depending on the original detected SDB-transcript association (BMI unadjusted or BMI adjusted). If an SDB phenotype was associated with the metabolite (p-value<0.05), we interpreted this as validation of a SDB association with this metabolite via the transcript-level chain.

Sample characteristics

Characteristics of the MESA population that participated in the TOPMed omics study, the sleep study, and the smaller T-cells and monocytes analyses are provided in Table S1; characteristics of the HCHS/SOL participants with genetic and metabolite data are provided in Table S2.

MESA individuals are a multi-ethnic sample, 69 years old on average during MESA exam 5, and 52% female. HCHS/SOL individuals are from diverse Hispanic/Latino backgrounds with a mean age of 46 years during the baseline exam, and 59% female. SDB phenotypes were more severe in MESA, with average AHI=18.6, MinO2=83, and AvgO2=94.1, in contrast to HCHS/SOL with average AHI=6.4, MinO2=87.1, and AvgO2=96.4, consistent with the older age of the MESA sample. Characteristics of the WHI participants with RNA-seq data used to validate the transcript PRS are provided in Table S3. WHI individuals are from a multi-race and ethnic sample and are 80 years old on average at the Long-Life Study exam when RNA was extracted, and are all females.

SDB phenotypes for oxyhemoglobin saturation and AHI are linked to tissue-specific changes in the transcriptome

In MESA, we identified 96 and 24 differentially expressed transcripts (Table S4, S5 and Table S6, S7) with FDR p-value < 0.1 in unadjusted and adjusted BMI analyses, respectively, in the different cell types. Table 1 reports the top differentially expressed transcripts (FDR p-value <0.05). Three transcripts, AJUBA (Ajuba LIM Protein), ZNF665 (Zinc Finger Protein 665), and TMC3-AS1 (TMC3 Antisense RNA 1, a long non-coding RNA), are significantly associated with AvgO2 and AHI in both analyses, in the direction of reduced expression with worse SDB measures (higher AHI, lower AvgO2).

To visualize gene expression and compare across SDB traits and cell types, Spearman correlation of log-fold change in expression of all SDB-associated transcripts (n=96 transcripts FDR p<0.1) was illustrated with a heatmap for BMI unadjusted (Figure S3) and BMI adjusted (Figure S4) analyses, clustered using hierarchical clustering based on the correlation between the log-fold estimates. These results illustrate concordant and discordant patterns of differential gene expression by cell type (PBMCs, monocytes, and T-cells) and SDB trait (AvgO2, MinO2, and AHI). There are a few striking differences in gene expression, particularly the increased expression of FAM106A, DNAJA3, BCDIN3D-AS1, TGFBRAP1, BEND5, TTC24, TMC3-AS1,

• Download figure
• Open in new tab

Figure S3: Heatmap of estimated log-fold gene expression change with SDB phenotypes across tissues without BMI adjustment in MESA

• Download figure
• Open in new tab

Figure S4: Heatmap of estimated log-fold gene expression change with SDB phenotypes across tissues from BMI-adjusted analysis in MESA

LINC00235, SEC14L2, ARHGEF9, TSHZ1, and LA16c-312E8.4 in T-cells compared to monocytes and PBMCs. To further investigate the overall patterns in gene expression in relation to tissue type and SDB traits, a heatmap of the Spearman correlation of the log-fold expression estimates of SDB phenotypes was plotted Figure 3. Within cell types, the SDB traits AHI and MinO2 had the highest correlation for gene expression (Spearman R² between 0.91 to 0.97), whereas AvgO2 had lower correlations with AHI and MinO2, especially in monocytes. In addition, a heatmap of estimated log-fold gene expression change (FDR p-value < 0.1) with SDB phenotypes across tissues adjusted for BMI is shown in Figure S5. The correlations between the SDB effect estimates for gene expression across cell types are different, and generally higher, from the phenotypic correlations between the SDB phenotypes, which are at the range of 0.53 to 0.73 Spearman R² (Figure S6). When computing correlations over all genes, estimated associations between AvgO2 with gene expression had almost no correlation with the other phenotypes (Figure S7). The same patterns (although slightly attenuated) are observed in BMI adjusted analyses (Figure S8).

• Download figure
• Open in new tab

Figure S5: Spearman correlations between estimated log-fold changes in gene expression across SDB phenotypes and tissues in BMI adjusted analysis in MESA (top genes)

Heatmap illustrating the Spearman correlations of log-fold change of transcript expression by tissue type (monocytes, T-cells, PBMCs) and SDB phenotype (AvgO2, MinO2, AHI) in MESA. Correlations were computed over genes with FDR p<0.1. Color legend portrays Spearman R² (no/weak correlation = light yellow; complete/strong correlation = green). Estimated AHI effect sizes are flipped prior to computation of correlations so that they match the direction of MinO2 and AvgO2.

• Download figure
• Open in new tab

Figure S6. Correlation between the SDB phenotypes

Heatmap illustrating the Spearman correlations SDB phenotype (AvgO2, MinO2, AHI). Color legend portrays Spearman R² (no/weak correlation = light yellow; complete/strong correlation = green).

• Download figure
• Open in new tab

Figure S7: Spearman correlations between estimated log-fold changes in gene expression across SDB phenotypes and tissues in analysis without BMI adjustment (all genes).

Heatmap illustrating the Spearman correlations of log-fold change of transcript expression by tissue type (monocytes, T-cells, PBMCs) and SDB phenotype (AvgO2, MinO2, AHI). Color legend portrays Spearman R² (no/weak correlation = light yellow; complete/strong correlation = green).

• Download figure
• Open in new tab

Figure S8: Spearman correlations between estimated log-fold changes in gene expression across SDB phenotypes and tissues in analysis BMI adjustment (all genes).

Heatmap illustrating the Spearman correlations of log-fold change of transcript expression by tissue type (monocytes, T-cells, PBMCs) and SDB phenotype (AvgO2, MinO2, AHI). Color legend portrays Spearman R² (no/weak correlation = light yellow; complete/strong correlation = green).

• Download figure
• Open in new tab

Figure 3: Spearman correlations between estimated log-fold changes in gene expression across SDB phenotypes and blood cell types without BMI adjustment in MESA

Heatmap illustrating the Spearman correlations of log-fold change of transcript expression by tissue type (monocytes, T-cells, PBMCs) and SDB phenotype (AvgO2, MinO2, AHI) in MESA. Correlations were computed over genes with FDR p<0.1. Color legend portrays Spearman R² (no/weak correlation = light yellow; complete/strong correlation = green). Estimated AHI effect sizes were flipped prior to computation of correlations so that they match the direction of MinO2 and AvgO2.

Construction and validation of transcript PRS

We constructed tPRS for gene expression in monocytes using a few methods, focusing on transcripts that were associated with SDB exposures in our analysis. The performance of constructed tPRS was evaluated against whole-blood gene expression levels in n=1,269 WHI participants. Figure S9 visualizes the results, demonstrating that tPRS constructed using the clump and threshold for genome-wide SNPs, including trans-eQTLs and tPRS focusing on cis-eQTLs have similar results, and the same generalization rate as that of the prediXcan-based tPRS. However, prediXcan tPRS had opposite direction of association with one of the transcripts in WHI, and, both cis-eQTLs based tPRS (prediXcan and clump and threshold) were not available for some transcripts due to lack of transcript-associated SNPs near the coding region. Thus, we moved forward with the genome-wide approach. Of the 96 tPRS (BMI unadjusted analysis) and 24 tPRS (BMI-adjusted analysis) tested, 26 and 9 tPRS were associated (p<0.017) with gene expression (Tables S8-S9) in whole blood and considered validated as IVs.

• Download figure
• Open in new tab

Figure S9. Comparison of the associations between monocyte-based tPRSs and whole-blood gene expression in WHI

For each transcript associated with an SDB phenotype, the figure provides the estimated association effect, 95% confidence interval, and p-value, of tPRSs constructed in different approaches with whole-blood transcript expression in WHI. tPRS1 and tPRS3 were constructed using the clump and threshold approach implemented in PRSice2 using summary statistics from GWAS of transcript expression in monocytes in MESA, and with clumping guided by LD in MESA (the same individuals used for GWAS). tPRS1 allows for genome-wide SNPs, and tPRS3 focused on cis-eQTLs. tPRS2 is the prediXcan model (also using cis-eQLTs only). For tPRS1 and 3 we considered three p-value threshold (5x10^-8, 10^-7, and 10^-6), and the tPRS with smallest p-value is displayed. PRS associations were estimated in models adjusted for sex, age, study site, race/ethnic background, batch effects, and 11 ancestral principal components.

Evidence of causal association between transcripts and SDB phenotypes

We tested the association of the validated tPRS, constructed in a cell-specific manner, with SDB phenotypes in HCHS/SOL (Table S10). Of the 26 tested in BMI unadjusted analysis, 3 tPRS showed evidence of reverse association with SDB phenotypes (p-value<0.05), supporting a causal relationship between expression of these transcripts and SDB traits. Among them, the strongest association was of the tPRS for P2RX4 (Purinergic Receptor P2X 4) in PBMC, in its association with AvgO2, one standard deviation (SD) increase in the PRS was associated with increased 1.9% AvgO2. Additionally, tPRS for SEC14L2 (SEC14 Like Lipid Binding 2) was negatively associated with AHI in T-cells and tPRS for TUBB6 (Tubulin Beta 6 Class V) was positively associated with MinO2 in monocytes. After BMI adjustment, only P2RX4 remained positively associated with AvgO2 in PBMCs (p-value <0.05), as shown in Table S11.

Evidence of causal association between transcripts and metabolites

We tested the relation between each validated tPRS and metabolites in HCHS/SOL. The tPRS for P2RX4 and CTD-2366F13.1 (also known as MOCS2-DT, MOCS2 Divergent Transcript) were associated with a total of 6 and 7 metabolites in unadjusted BMI and adjusted BMI analyses (FDR p-value <0.05, Table S12 and Tables S13), respectively; the association “chains” are visualized in Figure 4. Of 7 metabolites, 3 of them (butyrylcarnitine, linoleoyl-arachidonoyl-glycerol (18:2/20:4), and palmitoleoyl-linoleoyl-glycerol (16:1/18:2)) were also associated with AvgO2 (Table S15). However, the AvgO2-metabolite associations did not remain after BMI adjustment, suggesting that BMI, rather than SDB, may be driving these associations (Table S16). Of the transcripts, P2RX4 had evidence of a complete chain of association with SDB and metabolites (p-value <0.05) in the BMI unadjusted analysis.

• Download figure
• Open in new tab

Figure 4. Identified association chains between AvgO2, transcripts, and metabolites

Diagram illustrating the “association chain” relationships between AvgO2, tPRS, and metabolites in BMI unadjusted and BMI adjusted analyses. Dark blue color indicates association between tPRS and metabolites in both BMI unadjusted and adjusted analyses; The light blue color indicates the association between tPRS and metabolite only in BMI adjusted analysis.