SmCCNet 2.0: A Comprehensive Tool for Multi-omics Network Inference with Shiny Visualization

Summary

Sparse multiple canonical correlation network analysis (SmCCNet) is a machine learning technique for integrating omics data along with a variable of interest (e.g., phenotype of complex disease), and reconstructing multi-omics networks that are specific to this variable. We present the second-generation SmCCNet (SmCCNet 2.0) that adeptly integrates single or multiple omics data types along with a quantitative or binary phenotype of interest. In addition, this new package offers a streamlined setup process that can be configured manually or automatically, ensuring a flexible and user-friendly experience.

Availability

This package is available in both CRAN: https://cran.r-project.org/web/packages/SmCCNet/index.html and Github: https://github.com/KechrisLab/SmCCNet under the MIT license. The network visualization tool is available at https://smccnet.shinyapps.io/smccnetnetwork/.

Cord blood DNA methylation of immune and lipid metabolism genes is associated with maternal triglycerides and child adiposity

OBJECTIVE

Fetal exposures may impact offspring epigenetic signatures and adiposity. The authors hypothesized that maternal metabolic traits associate with cord blood DNA methylation, which, in turn, associates with child adiposity.

METHODS

Fasting serum was obtained in 588 pregnant women (27-34 weeks' gestation), and insulin, glucose, high-density lipoprotein cholesterol, triglycerides, and free fatty acids were measured. Cord blood DNA methylation and child adiposity were measured at birth, 4-6 months, and 4-6 years. The association of maternal metabolic traits with DNA methylation (429,246 CpGs) for differentially methylated probes (DMPs) and regions (DMRs) was tested. The association of the first principal component of each DMR with child adiposity was tested, and mediation analysis was performed.

RESULTS

Maternal triglycerides were associated with the most DMPs and DMRs of all traits tested (261 and 198, respectively, false discovery rate < 0.05). DMRs were near genes involved in immune function and lipid metabolism. Triglyceride-associated CpGs were associated with child adiposity at 4-6 months (32 CpGs) and 4-6 years (2 CpGs). One, near CD226, was observed at both timepoints, mediating 10% and 22% of the relationship between maternal triglycerides and child adiposity at 4-6 months and 4-6 years, respectively.

CONCLUSIONS

DNA methylation may play a role in the association of maternal triglycerides and child adiposity.

Matrix Linear Models for connecting metabolite composition to individual characteristics

High-throughput metabolomics data provide a detailed molecular window into biological processes. We consider the problem of assessing how the association of metabolite levels with individual (sample) characteristics such as sex or treatment may depend on metabolite characteristics such as pathway. Typically this is one in a two-step process: In the first step we assess the association of each metabolite with individual characteristics. In the second step an enrichment analysis is performed by metabolite characteristics among significant associations. We combine the two steps using a bilinear model based on the matrix linear model (MLM) framework we have previously developed for high-throughput genetic screens. Our framework can estimate relationships in metabolites sharing known characteristics, whether categorical (such as type of lipid or pathway) or numerical (such as number of double bonds in triglycerides). We demonstrate how MLM offers flexibility and interpretability by applying our method to three metabolomic studies. We show that our approach can separate the contribution of the overlapping triglycerides characteristics, such as the number of double bonds and the number of carbon atoms. The proposed method have been implemented in the open-source Julia package, MatrixLM. Data analysis scripts with example data analyses are also available.

Accelerated epigenetic age at birth and child emotional and behavioura development in early childhood: a meta-analysis of four prospective cohort studies in ECHO

Background: 'Epigenetic clocks' have been developed to accurately predict chronologic gestational age and have been associated with child health outcomes in prior work.Methods: We meta-analysed results from four prospective U.S cohorts investigating the association between epigenetic age acceleration estimated using blood DNA methylation collected at birth and preschool age Childhood Behavior Checklist (CBCL) scores.Results: Epigenetic ageing was not significantly associated with CBCL total problem scores (β = 0.33, 95% CI: -0.95, 0.28) and DSM-oriented pervasive development problem scores (β = -0.23, 95% CI: -0.61, 0.15). No associations were observed for other DSM-oriented subscales.Conclusions: The meta-analysis results suggest that epigenetic gestational age acceleration is not associated with child emotional and behavioural functioning for preschool age group. These findings may relate to our study population, which includes two cohorts enriched for ASD and one preterm birth cohort.; future work should address the role of epigenetic age in child health in other study populations.Abbreviations: DNAm: DNA methylation; CBCL: Child Behavioral Checklist; ECHO: Environmental Influences on Child Health Outcomes; EARLI: Early Autism Risk Longitudinal Investigation; MARBLES: Markers of Autism Risk in Babies - Learning Early Signs; ELGAN: Extremely Low Gestational Age Newborns; ASD: autism spectrum disorder; BMI: body mass index; DSM: Diagnostic and Statistical Manual of Mental Disorders.

Sex-based differences in placental DNA methylation profiles related to gestational age: an NIH ECHO meta-analysis

The placenta undergoes many changes throughout gestation to support the evolving needs of the foetus. There is also a growing appreciation that male and female foetuses develop differently in utero, with unique epigenetic changes in placental tissue. Here, we report meta-analysed sex-specific associations between gestational age and placental DNA methylation from four cohorts in the National Institutes of Health (NIH) Environmental influences on Child Health Outcomes (ECHO) Programme (355 females/419 males, gestational ages 23-42 weeks). We identified 407 cytosine-guanine dinucleotides (CpGs) in females and 794 in males where placental methylation levels were associated with gestational age. After cell-type adjustment, 55 CpGs in females and 826 in males were significant. These were enriched for biological processes critical to the immune system in females and transmembrane transport in males. Our findings are distinct between the sexes: in females, associations with gestational age are largely explained by differences in placental cellular composition, whereas in males, gestational age is directly associated with numerous alterations in methylation levels.

TreeKernel: interpretable kernel machine tests for interactions between -omics and clinical predictors with applications to metabolomics and COPD phenotypes

BACKGROUND

In this paper, we are interested in interactions between a high-dimensional -omics dataset and clinical covariates. The goal is to evaluate the relationship between a phenotype of interest and a high-dimensional omics pathway, where the effect of the omics data depends on subjects' clinical covariates (age, sex, smoking status, etc.). For instance, metabolic pathways can vary greatly between sexes which may also change the relationship between certain metabolic pathways and a clinical phenotype of interest. We propose partitioning the clinical covariate space and performing a kernel association test within those partitions. To illustrate this idea, we focus on hierarchical partitions of the clinical covariate space and kernel tests on metabolic pathways.

RESULTS

We see that our proposed method outperforms competing methods in most simulation scenarios. It can identify different relationships among clinical groups with higher power in most scenarios while maintaining a proper Type I error rate. The simulation studies also show a robustness to the grouping structure within the clinical space. We also apply the method to the COPDGene study and find several clinically meaningful interactions between metabolic pathways, the clinical space, and lung function.

CONCLUSION

TreeKernel provides a simple and interpretable process for testing for relationships between high-dimensional omics data and clinical outcomes in the presence of interactions within clinical cohorts. The method is broadly applicable to many studies.

Quantifying the spatial clustering characteristics of radiographic emphysema explains variability in pulmonary function

Quantitative assessment of emphysema in CT scans has mostly focused on calculating the percentage of lung tissue that is deemed abnormal based on a density thresholding strategy. However, this overall measure of disease burden discards virtually all the spatial information encoded in the scan that is implicitly utilized in a visual assessment. This simplification is likely grouping heterogenous disease patterns and is potentially obscuring clinical phenotypes and variable disease outcomes. To overcome this, several methods that attempt to quantify heterogeneity in emphysema distribution have been proposed. Here, we compare three of those: one based on estimating a power law for the size distribution of contiguous emphysema clusters, a second that looks at the number of emphysema-to-emphysema voxel adjacencies, and a third that applies a parametric spatial point process model to the emphysema voxel locations. This was done using data from 587 individuals from Phase 1 of COPDGene that had an inspiratory CT scan and plasma protein abundance measurements. The associations between these imaging metrics and visual assessment with clinical measures (FEV[Formula: see text], FEV[Formula: see text]-FVC ratio, etc.) and plasma protein biomarker levels were evaluated using a variety of regression models. Our results showed that a selection of spatial measures had the ability to discern heterogeneous patterns among CTs that had similar emphysema burdens. The most informative quantitative measure, average cluster size from the point process model, showed much stronger associations with nearly every clinical outcome examined than existing CT-derived emphysema metrics and visual assessment. Moreover, approximately 75% more plasma biomarkers were found to be associated with an emphysema heterogeneity phenotype when accounting for spatial clustering measures than when they were excluded.

Prenatal exposures to per- and polyfluoroalkyl substances and epigenetic aging in umbilical cord blood: The Healthy Start study

BACKGROUND

Per- and polyfluoroalkyl substances (PFAS) are ubiquitous, environmentally persistent chemicals, and prenatal exposures have been associated with adverse child health outcomes. Prenatal PFAS exposure may lead to epigenetic age acceleration (EAA), defined as the discrepancy between an individual's chronologic and epigenetic or biological age.

OBJECTIVES

We estimated associations of maternal serum PFAS concentrations with EAA in umbilical cord blood DNA methylation using linear regression, and a multivariable exposure-response function of the PFAS mixture using Bayesian kernel machine regression.

METHODS

Five PFAS were quantified in maternal serum (median: 27 weeks of gestation) among 577 mother-infant dyads from a prospective cohort. Cord blood DNA methylation data were assessed with the Illumina HumanMethylation450 array. EAA was calculated as the residuals from regressing gestational age on epigenetic age, calculated using a cord-blood specific epigenetic clock. Linear regression tested for associations between each maternal PFAS concentration with EAA. Bayesian kernel machine regression with hierarchical selection estimated an exposure-response function for the PFAS mixture.

RESULTS

In single pollutant models we observed an inverse relationship between perfluorodecanoate (PFDA) and EAA (-0.148 weeks per log-unit increase, 95% CI: -0.283, -0.013). Mixture analysis with hierarchical selection between perfluoroalkyl carboxylates and sulfonates indicated the carboxylates had the highest group posterior inclusion probability (PIP), or relative importance. Within this group, PFDA had the highest conditional PIP. Univariate predictor-response functions indicated PFDA and perfluorononanoate were inversely associated with EAA, while perfluorohexane sulfonate had a positive association with EAA.

CONCLUSIONS

Maternal mid-pregnancy serum concentrations of PFDA were negatively associated with EAA in cord blood, suggesting a pathway by which prenatal PFAS exposures may affect infant development. No significant associations were observed with other PFAS. Mixture models suggested opposite directions of association between perfluoroalkyl sulfonates and carboxylates. Future studies are needed to determine the importance of neonatal EAA for later child health outcomes.

Intrauterine Smoke Exposure, microRNA Expression during Human Lung Development, and Childhood Asthma

Intrauterine smoke (IUS) exposure during early childhood has been associated with a number of negative health consequences, including reduced lung function and asthma susceptibility. The biological mechanisms underlying these associations have not been established. MicroRNAs regulate the expression of numerous genes involved in lung development. Thus, investigation of the impact of IUS on miRNA expression during human lung development may elucidate the impact of IUS on post-natal respiratory outcomes. We sought to investigate the effect of IUS exposure on miRNA expression during early lung development. We hypothesized that miRNA-mRNA networks are dysregulated by IUS during human lung development and that these miRNAs may be associated with future risk of asthma and allergy. Human fetal lung samples from a prenatal tissue retrieval program were tested for differential miRNA expression with IUS exposure (measured using placental cotinine concentration). RNA was extracted and miRNA-sequencing was performed. We performed differential expression using IUS exposure, with covariate adjustment. We also considered the above model with an additional sex-by-IUS interaction term, allowing IUS effects to differ by male and female samples. Using paired gene expression profiles, we created sex-stratified miRNA-mRNA correlation networks predictive of IUS using DIABLO. We additionally evaluated whether miRNAs were associated with asthma and allergy outcomes in a cohort of childhood asthma. We profiled pseudoglandular lung miRNA in n = 298 samples, 139 (47%) of which had evidence of IUS exposure. Of 515 miRNAs, 25 were significantly associated with intrauterine smoke exposure (q-value < 0.10). The IUS associated miRNAs were correlated with well-known asthma genes (e.g., ORM1-Like Protein 3, ORDML3) and enriched in disease-relevant pathways (oxidative stress). Eleven IUS-miRNAs were also correlated with clinical measures (e.g., Immunoglobulin E andlungfunction) in children with asthma, further supporting their likely disease relevance. Lastly, we found substantial differences in IUS effects by sex, finding 95 significant IUS-miRNAs in male samples, but only four miRNAs in female samples. The miRNA-mRNA correlation networks were predictive of IUS (AUC = 0.78 in males and 0.86 in females) and suggested that IUS-miRNAs are involved in regulation of disease-relevant genes (e.g., A disintegrin and metalloproteinase domain 19 (ADAM19), LBH regulator of WNT signaling (LBH)) and sex hormone signaling (Coactivator associated methyltransferase 1(CARM1)). Our study demonstrated differential expression of miRNAs by IUS during early prenatal human lung development, which may be modified by sex. Based on their gene targets and correlation to clinical asthma and atopy outcomes, these IUS-miRNAs may be relevant for subsequent allergy and asthma risk. Our study provides insight into the impact of IUS in human fetal lung transcriptional networks and on the developmental origins of asthma and allergic disorders.

Deep learning on graphs for multi-omics classification of COPD

Network approaches have successfully been used to help reveal complex mechanisms of diseases including Chronic Obstructive Pulmonary Disease (COPD). However despite recent advances, we remain limited in our ability to incorporate protein-protein interaction (PPI) network information with omics data for disease prediction. New deep learning methods including convolution Graph Neural Network (ConvGNN) has shown great potential for disease classification using transcriptomics data and known PPI networks from existing databases. In this study, we first reconstructed the COPD-associated PPI network through the AhGlasso (Augmented High-Dimensional Graphical Lasso Method) algorithm based on one independent transcriptomics dataset including COPD cases and controls. Then we extended the existing ConvGNN methods to successfully integrate COPD-associated PPI, proteomics, and transcriptomics data and developed a prediction model for COPD classification. This approach improves accuracy over several conventional classification methods and neural networks that do not incorporate network information. We also demonstrated that the updated COPD-associated network developed using AhGlasso further improves prediction accuracy. Although deep neural networks often achieve superior statistical power in classification compared to other methods, it can be very difficult to explain how the model, especially graph neural network(s), makes decisions on the given features and identifies the features that contribute the most to prediction generally and individually. To better explain how the spectral-based Graph Neural Network model(s) works, we applied one unified explainable machine learning method, SHapley Additive exPlanations (SHAP), and identified CXCL11, IL-2, CD48, KIR3DL2, TLR2, BMP10 and several other relevant COPD genes in subnetworks of the ConvGNN model for COPD prediction. Finally, Gene Ontology (GO) enrichment analysis identified glycosaminoglycan, heparin signaling, and carbohydrate derivative signaling pathways significantly enriched in the top important gene/proteins for COPD classifications.

Novel Metabolic Subtypes in Pregnant Women and Risk of Early Childhood Obesity in Offspring

Importance

The in utero metabolic milieu is associated with offspring adiposity. Standard definitions of maternal obesity (according to prepregnancy body mass index [BMI]) and gestational diabetes (GDM) may not be adequate to capture subtle yet important differences in the intrauterine environment that could be involved in programming.

Objectives

To identify maternal metabolic subgroups during pregnancy and to examine associations of subgroup classification with adiposity traits in their children.

Design, Setting, and Participants

This cohort study included mother-offspring pairs in the Healthy Start prebirth cohort (enrollment: 2010-2014) recruited from University of Colorado Hospital obstetrics clinics in Aurora, Colorado. Follow-up of women and children is ongoing. Data were analyzed from March to December 2022.

Exposures

Metabolic subtypes of pregnant women ascertained by applying k-means clustering on 7 biomarkers and 2 biomarker indices measured at approximately 17 gestational weeks: glucose, insulin, Homeostatic Model Assessment for Insulin Resistance, total cholesterol, high-density lipoprotein cholesterol (HDL-C), triglycerides, free fatty acids (FFA), HDL-C:triglycerides ratio, and tumor necrosis factor α.

Main Outcomes and Measures

Offspring birthweight z score and neonatal fat mass percentage (FM%). In childhood at approximately 5 years of age, offspring BMI percentile, FM%, BMI in the 95th percentile or higher, and FM% in the 95th percentile or higher.

Results

A total of 1325 pregnant women (mean [SD] age, 27.8 [6.2 years]; 322 [24.3%] Hispanic, 207 non-Hispanic Black [15.6%], and 713 [53.8%] non-Hispanic White), and 727 offspring with anthropometric data measured in childhood (mean [SD] age 4.81 [0.72] years, 48% female) were included. We identified the following 5 maternal metabolic subgroups: reference (438 participants), high HDL-C (355 participants), dyslipidemic-high triglycerides (182 participants), dyslipidemic-high FFA (234 participants), and insulin resistant (IR)-hyperglycemic (116 participants). Compared with the reference subgroup, women in the IR-hyperglycemic and dyslipidemic-high FFA subgroups had offspring with 4.27% (95% CI, 1.94-6.59) and 1.96% (95% CI, 0.45-3.47) greater FM% during childhood, respectively. There was a higher risk of high FM% among offspring of the IR-hyperglycemic (relative risk, 8.7; 95% CI, 2.7-27.8) and dyslipidemic-high FFA (relative risk, 3.4; 95% CI, 1.0-11.3) subgroups; this risk was of greater magnitude compared with prepregnancy obesity alone, GDM alone, or both conditions.

Conclusions and Relevance

In this cohort study, an unsupervised clustering approach revealed distinct metabolic subgroups of pregnant women. These subgroups exhibited differences in risk of offspring adiposity in early childhood. Such approaches have the potential to refine understanding of the in utero metabolic milieu, with utility for capturing variation in sociocultural, anthropometric, and biochemical risk factors for offspring adiposity.

Genome-wide analysis of oxylipins and oxylipin profiles in a pediatric population

Background

Oxylipins are inflammatory biomarkers derived from omega-3 and-6 fatty acids implicated in inflammatory diseases but have not been studied in a genome-wide association study (GWAS). The aim of this study was to identify genetic loci associated with oxylipins and oxylipin profiles to identify biologic pathways and therapeutic targets for oxylipins.

Methods

We conducted a GWAS of plasma oxylipins in 316 participants in the Diabetes Autoimmunity Study in the Young (DAISY). DNA samples were genotyped using the TEDDY-T1D Exome array, and additional variants were imputed using the Trans-Omics for Precision Medicine (TOPMed) multi-ancestry reference panel. Principal components analysis of 36 plasma oxylipins was used to capture oxylipin profiles. PC1 represented linoleic acid (LA)- and alpha-linolenic acid (ALA)-related oxylipins, and PC2 represented arachidonic acid (ARA)-related oxylipins. Oxylipin PC1, PC2, and the top five loading oxylipins from each PC were used as outcomes in the GWAS (genome-wide significance: p < 5×10-8).

Results

The SNP rs143070873 was associated with (p < 5×10-8) the LA-related oxylipin 9-HODE, and rs6444933 (downstream of CLDN11) was associated with the LA-related oxylipin 13 S-HODE. A locus between MIR1302-7 and LOC100131146, rs10118380 and an intronic variant in TRPM3 were associated with the ARA-related oxylipin 11-HETE. These loci are involved in inflammatory signaling cascades and interact with PLA2, an initial step to oxylipin biosynthesis.

Conclusion

Genetic loci involved in inflammation and oxylipin metabolism are associated with oxylipin levels.

Propensity scores as a novel method to guide sample allocation and minimize batch effects during the design of high throughput experiments

BACKGROUND

We developed a novel approach to minimize batch effects when assigning samples to batches. Our algorithm selects a batch allocation, among all possible ways of assigning samples to batches, that minimizes differences in average propensity score between batches. This strategy was compared to randomization and stratified randomization in a case-control study (30 per group) with a covariate (case vs control, represented as β1, set to be null) and two biologically relevant confounding variables (age, represented as β2, and hemoglobin A1c (HbA1c), represented as β3). Gene expression values were obtained from a publicly available dataset of expression data obtained from pancreas islet cells. Batch effects were simulated as twice the median biological variation across the gene expression dataset and were added to the publicly available dataset to simulate a batch effect condition. Bias was calculated as the absolute difference between observed betas under the batch allocation strategies and the true beta (no batch effects). Bias was also evaluated after adjustment for batch effects using ComBat as well as a linear regression model. In order to understand performance of our optimal allocation strategy under the alternative hypothesis, we also evaluated bias at a single gene associated with both age and HbA1c levels in the 'true' dataset (CAPN13 gene).

RESULTS

Pre-batch correction, under the null hypothesis (β1), maximum absolute bias and root mean square (RMS) of maximum absolute bias, were minimized using the optimal allocation strategy. Under the alternative hypothesis (β2 and β3 for the CAPN13 gene), maximum absolute bias and RMS of maximum absolute bias were also consistently lower using the optimal allocation strategy. ComBat and the regression batch adjustment methods performed well as the bias estimates moved towards the true values in all conditions under both the null and alternative hypotheses. Although the differences between methods were less pronounced following batch correction, estimates of bias (average and RMS) were consistently lower using the optimal allocation strategy under both the null and alternative hypotheses.

CONCLUSIONS

Our algorithm provides an extremely flexible and effective method for assigning samples to batches by exploiting knowledge of covariates prior to sample allocation.

An Oxylipin-Related Nutrient Pattern and Risk of Type 1 Diabetes in the Diabetes Autoimmunity Study in the Young (DAISY)

Oxylipins, pro-inflammatory and pro-resolving lipid mediators, are associated with the risk of type 1 diabetes (T1D) and may be influenced by diet. This study aimed to develop a nutrient pattern related to oxylipin profiles and test their associations with the risk of T1D among youth. The nutrient patterns were developed with a reduced rank regression in a nested case-control study (n = 335) within the Diabetes Autoimmunity Study in the Young (DAISY), a longitudinal cohort of children at risk of T1D. The oxylipin profiles (adjusted for genetic predictors) were the response variables. The nutrient patterns were tested in the case-control study (n = 69 T1D cases, 69 controls), then validated in the DAISY cohort using a joint Cox proportional hazards model (n = 1933, including 81 T1D cases). The first nutrient pattern (NP1) was characterized by low beta cryptoxanthin, flavanone, vitamin C, total sugars and iron, and high lycopene, anthocyanidins, linoleic acid and sodium. After adjusting for T1D family history, the HLA genotype, sex and race/ethnicity, NP1 was associated with a lower risk of T1D in the nested case-control study (OR: 0.44, p = 0.0126). NP1 was not associated with the risk of T1D (HR: 0.54, p-value = 0.1829) in the full DAISY cohort. Future studies are needed to confirm the nested case-control findings and investigate the modifiable factors for oxylipins.

RCFGL: Rapid Condition adaptive Fused Graphical Lasso and application to modeling brain region co-expression networks

Inferring gene co-expression networks is a useful process for understanding gene regulation and pathway activity. The networks are usually undirected graphs where genes are represented as nodes and an edge represents a significant co-expression relationship. When expression data of multiple (p) genes in multiple (K) conditions (e.g., treatments, tissues, strains) are available, joint estimation of networks harnessing shared information across them can significantly increase the power of analysis. In addition, examining condition-specific patterns of co-expression can provide insights into the underlying cellular processes activated in a particular condition. Condition adaptive fused graphical lasso (CFGL) is an existing method that incorporates condition specificity in a fused graphical lasso (FGL) model for estimating multiple co-expression networks. However, with computational complexity of O(p2K log K), the current implementation of CFGL is prohibitively slow even for a moderate number of genes and can only be used for a maximum of three conditions. In this paper, we propose a faster alternative of CFGL named rapid condition adaptive fused graphical lasso (RCFGL). In RCFGL, we incorporate the condition specificity into another popular model for joint network estimation, known as fused multiple graphical lasso (FMGL). We use a more efficient algorithm in the iterative steps compared to CFGL, enabling faster computation with complexity of O(p2K) and making it easily generalizable for more than three conditions. We also present a novel screening rule to determine if the full network estimation problem can be broken down into estimation of smaller disjoint sub-networks, thereby reducing the complexity further. We demonstrate the computational advantage and superior performance of our method compared to two non-condition adaptive methods, FGL and FMGL, and one condition adaptive method, CFGL in both simulation study and real data analysis. We used RCFGL to jointly estimate the gene co-expression networks in different brain regions (conditions) using a cohort of heterogeneous stock rats. We also provide an accommodating C and Python based package that implements RCFGL.

Fat content in infant mesenchymal stem cells prospectively associates with childhood adiposity and fasting glucose

OBJECTIVE

In human studies, new model systems are needed for improved mechanistic investigation of developmental predisposition for metabolic disease but also to serve as benchmarks in early life prevention or intervention efforts. In this regard, human infant umbilical cord-derived mesenchymal stem cells (MSCs) are an emerging tool. However, long-term clinical relevance to in vivo markers of metabolic disease is unknown.

METHODS

In a cohort of 124 mother/child dyads, this study tested the hypothesis that triglyceride content (TG) of infant MSCs undergoing adipogenesis in vitro (MSC-TG) is associated with in vivo adiposity (percent fat mass) from birth to early childhood and with fasting glucose and insulin in early childhood.

RESULTS

MSC-TG was positively associated with in vivo child adiposity at birth, age 4 to 6 months, and age 4 to 6 years. MSC-TG was associated with fasting glucose, but not insulin, at 4 to 6 years. Importantly, MSC-TG explained an additional 13% variance in child adiposity at 4 to 6 years, after accounting for other established birth predictors (weight and percent fat mass at birth) and other established covariates related to child adiposity (e.g., breastfeeding exposure, physical activity).

CONCLUSIONS

This work demonstrates the strength of the MSC model for predicting offspring metabolic phenotype into childhood, even when considering the important contribution of other early life risk factors.

Ambient air pollution during pregnancy and DNA methylation in umbilical cord blood, with potential mediation of associations with infant adiposity: The Healthy Start study

BACKGROUND

Prenatal exposure to ambient air pollution has been associated with adverse offspring health outcomes. Childhood health effects of prenatal exposures may be mediated through changes to DNA methylation detectable at birth.

METHODS

Among 429 non-smoking women in a cohort study of mother-infant pairs in Colorado, USA, we estimated associations between prenatal exposure to ambient fine particulate matter (PM2.5) and ozone (O3), and epigenome-wide DNA methylation of umbilical cord blood cells at delivery (2010-2014). We calculated average PM2.5 and O3 in each trimester of pregnancy and the full pregnancy using inverse-distance-weighted interpolation. We fit linear regression models adjusted for potential confounders and cell proportions to estimate associations between air pollutants and methylation at each of 432,943 CpGs. Differentially methylated regions (DMRs) were identified using comb-p. Previously in this cohort, we reported positive associations between 3rd trimester O3 exposure and infant adiposity at 5 months of age. Here, we quantified the potential for mediation of that association by changes in DNA methylation in cord blood.

RESULTS

We identified several DMRs for each pollutant and period of pregnancy. The greatest number of significant DMRs were associated with third trimester PM2.5 (21 DMRs). No single CpGs were associated with air pollutants at a false discovery rate <0.05. We found that up to 8% of the effect of 3rd trimester O3 on 5-month adiposity may be mediated by locus-specific methylation changes, but mediation estimates were not statistically significant.

CONCLUSIONS

Differentially methylated regions in cord blood were identified in association with maternal exposure to PM2.5 and O3. Genes annotated to the significant sites played roles in cardiometabolic disease, immune function and inflammation, and neurologic disorders. We found limited evidence of mediation by DNA methylation of associations between third trimester O3 exposure and 5-month infant adiposity.

Altered Immunoglobulin Repertoire and Decreased IgA Somatic Hypermutation in the Gut during Chronic HIV-1 Infection

Humoral immune perturbations contribute to pathogenic outcomes in persons with HIV-1 infection (PWH). Gut barrier dysfunction in PWH is associated with microbial translocation and alterations in microbial communities (dysbiosis), and IgA, the most abundant immunoglobulin (Ig) isotype in the gut, is involved in gut homeostasis by interacting with the microbiome. We determined the impact of HIV-1 infection on the antibody repertoire in the gastrointestinal tract by comparing Ig gene utilization and somatic hypermutation (SHM) in colon biopsies from PWH (n = 19) versus age and sex-matched controls (n = 13). We correlated these Ig parameters with clinical, immunological, microbiome and virological data. Gene signatures of enhanced B cell activation were accompanied by skewed frequencies of multiple Ig Variable genes in PWH. PWH showed decreased frequencies of SHM in IgA and possibly IgG, with a substantial loss of highly mutated IgA sequences. The decline in IgA SHM in PWH correlated with gut CD4+ T cell loss and inversely correlated with mucosal inflammation and microbial translocation. Diminished gut IgA SHM in PWH was driven by transversion mutations at A or T deoxynucleotides, suggesting a defect not at the AID/APOBEC3 deamination step but at later stages of IgA SHM. These results expand our understanding of humoral immune perturbations in PWH that could have important implications in understanding mucosal immune defects in individuals with chronic HIV-1 infection. IMPORTANCE The gut is a major site of early HIV-1 replication and pathogenesis. Extensive CD4+ T cell depletion in this compartment results in a compromised epithelial barrier that facilitates the translocation of microbes into the underlying lamina propria and systemic circulation, resulting in chronic immune activation. To date, the consequences of microbial translocation on the mucosal humoral immune response (or vice versa) remains poorly integrated into the panoply of mucosal immune defects in PWH. We utilized next-generation sequencing approaches to profile the Ab repertoire and ascertain frequencies of somatic hypermutation in colon biopsies from antiretroviral therapy-naive PWH versus controls. Our findings identify perturbations in the Ab repertoire of PWH that could contribute to development or maintenance of dysbiosis. Moreover, IgA mutations significantly decreased in PWH and this was associated with adverse clinical outcomes. These data may provide insight into the mechanisms underlying impaired Ab-dependent gut homeostasis during chronic HIV-1 infection.

Changes in the Coexpression of Innate Immunity Genes During Persistent Islet Autoimmunity Are Associated With Progression of Islet Autoimmunity: Diabetes Autoimmunity Study in the Young (DAISY)

Longitudinal changes in gene expression during islet autoimmunity (IA) may provide insight into biological processes that explain progression to type 1 diabetes (T1D). We identified individuals from Diabetes Autoimmunity Study in the Young (DAISY) who developed IA, autoantibodies present on two or more visits. Illumina's NovaSeq 6000 was used to quantify gene expression in whole blood. With linear mixed models we tested for changes in expression after IA that differed across individuals who progressed to T1D (progressors) (n = 25), reverted to an autoantibody-negative stage (reverters) (n = 47), or maintained IA positivity but did not develop T1D (maintainers) (n = 66). Weighted gene coexpression network analysis was used to identify coexpression modules. Gene Ontology pathway analysis of the top 150 differentially expressed genes (nominal P < 0.01) identified significantly enriched pathways including leukocyte activation involved in immune response, innate immune response, and regulation of immune response. We identified a module of 14 coexpressed genes with roles in the innate immunity. The hub gene, LTF, is known to have immunomodulatory properties. Another gene within the module, CAMP, is potentially relevant based on its role in promoting β-cell survival in a murine model. Overall, results provide evidence of alterations in expression of innate immune genes prior to onset of T1D.

Mechanism-aware imputation: a two-step approach in handling missing values in metabolomics

When analyzing large datasets from high-throughput technologies, researchers often encounter missing quantitative measurements, which are particularly frequent in metabolomics datasets. Metabolomics, the comprehensive profiling of metabolite abundances, are typically measured using mass spectrometry technologies that often introduce missingness via multiple mechanisms: (1) the metabolite signal may be smaller than the instrument limit of detection; (2) the conditions under which the data are collected and processed may lead to missing values; (3) missing values can be introduced randomly. Missingness resulting from mechanism (1) would be classified as Missing Not At Random (MNAR), that from mechanism (2) would be Missing At Random (MAR), and that from mechanism (3) would be classified as Missing Completely At Random (MCAR). Two common approaches for handling missing data are the following: (1) omit missing data from the analysis; (2) impute the missing values. Both approaches may introduce bias and reduce statistical power in downstream analyses such as testing metabolite associations with clinical variables. Further, standard imputation methods in metabolomics often ignore the mechanisms causing missingness and inaccurately estimate missing values within a data set. We propose a mechanism-aware imputation algorithm that leverages a two-step approach in imputing missing values. First, we use a random forest classifier to classify the missing mechanism for each missing value in the data set. Second, we impute each missing value using imputation algorithms that are specific to the predicted missingness mechanism (i.e., MAR/MCAR or MNAR). Using complete data, we conducted simulations, where we imposed different missingness patterns within the data and tested the performance of combinations of imputation algorithms. Our proposed algorithm provided imputations closer to the original data than those using only one imputation algorithm for all the missing values. Consequently, our two-step approach was able to reduce bias for improved downstream analyses.

Granzyme B+ CD4 T cells accumulate in the colon during chronic HIV-1 infection

Chronic HIV-1 infection results in the sustained disruption of gut homeostasis culminating in alterations in microbial communities (dysbiosis) and increased microbial translocation. Major questions remain on how interactions between translocating microbes and gut immune cells impact HIV-1-associated gut pathogenesis. We previously reported that in vitro exposure of human gut cells to enteric commensal bacteria upregulated the serine protease and cytotoxic marker Granzyme B (GZB) in CD4 T cells, and GZB expression was further increased in HIV-1-infected CD4 T cells. To determine if these in vitro findings extend in vivo, we evaluated the frequencies of GZB⁺ CD4 T cells in colon biopsies and peripheral blood of untreated, chronically infected people with HIV-1 (PWH). Colon and blood GZB⁺ CD4 T cells were found at significantly higher frequencies in PWH. Colon, but not blood, GZB⁺ CD4 T cell frequencies were associated with gut and systemic T cell activation and Prevotella species abundance. In vitro, commensal bacteria upregulated GZB more readily in gut versus blood or tonsil-derived CD4 T cells, particularly in inflammatory T helper 17 cells. Bacteria-induced GZB expression in gut CD4 T cells required the presence of accessory cells, the IL-2 pathway and in part, MHC Class II. Overall, we demonstrate that GZB⁺ CD4 T cells are prevalent in the colon during chronic HIV-1 infection and may emerge following interactions with translocated bacteria in an IL-2 and MHC Class II-dependent manner. Associations between GZB⁺ CD4 T cells, dysbiosis and T cell activation suggest that GZB⁺ CD4 T cells may contribute to gut HIV-1 pathogenesis.

Sex-Specific Differences in MicroRNA Expression During Human Fetal Lung Development

Background: Sex-specific differences in fetal lung maturation have been well described; however, little is known about the sex-specific differences in microRNA (miRNA) expression during human fetal lung development. Interestingly, many adult chronic lung diseases also demonstrate sex-specific differences in prevalence. The developmental origins of health and disease hypothesis suggests that these sex-specific differences in fetal lung development may influence disease susceptibility later in life. In this study, we performed miRNA sequencing on human fetal lung tissue samples to investigate differential expression of miRNAs between males and females in the pseudoglandular stage of lung development. We hypothesized that differences in miRNA expression are present between sexes in early human lung development and may contribute to the sex-specific differences seen in pulmonary diseases later in life.

Methods: RNA was isolated from human fetal lung tissue samples for miRNA sequencing. The count of each miRNA was modeled by sex using negative binomial regression models in DESeq2, adjusting for post-conception age, age², smoke exposure, batch, and RUV factors. We tested for differential expression of miRNAs by sex, and for the presence of sex-by-age interactions to determine if miRNA expression levels by age were distinct between males and females.

Results: miRNA expression profiles were generated on 298 samples (166 males and 132 females). Of the 809 miRNAs expressed in human fetal lung tissue during the pseudoglandular stage of lung development, we identified 93 autosomal miRNAs that were significantly differentially expressed by sex and 129 miRNAs with a sex-specific pattern of miRNA expression across the course of the pseudoglandular period.

Conclusion: Our study demonstrates differential expression of numerous autosomal miRNAs between the male and female developing human lung. Additionally, the expression of some miRNAs are modified by age across the pseudoglandular stage in a sex-specific way. Some of these differences in miRNA expression may impact susceptibility to pulmonary disease later in life. Our results suggest that sex-specific miRNA expression during human lung development may be a potential mechanism to explain sex-specific differences in lung development and may impact subsequent disease susceptibility.

A Metabolomic Severity Score for Airflow Obstruction and Emphysema

Chronic obstructive pulmonary disease (COPD) is a disease with marked metabolic disturbance. Previous studies have shown the association between single metabolites and lung function for COPD, but whether a combination of metabolites could predict phenotype is unknown. We developed metabolomic severity scores using plasma metabolomics from the Metabolon platform from two US cohorts of ever-smokers: the Subpopulations and Intermediate Outcome Measures in COPD Study (SPIROMICS) (n = 648; training/testing cohort; 72% non-Hispanic, white; average age 63 years) and the COPDGene Study (n = 1120; validation cohort; 92% non-Hispanic, white; average age 67 years). Separate adaptive LASSO (adaLASSO) models were used to model forced expiratory volume at one second (FEV₁) and MESA-adjusted lung density using 762 metabolites common between studies. Metabolite coefficients selected by the adaLASSO procedure were used to create a metabolomic severity score (metSS) for each outcome. A total of 132 metabolites were selected to create a metSS for FEV₁. The metSS-only models explained 64.8% and 31.7% of the variability in FEV₁ in the training and validation cohorts, respectively. For MESA-adjusted lung density, 129 metabolites were selected, and metSS-only models explained 59.0% of the variability in the training cohort and 17.4% in the validation cohort. Regression models including both clinical covariates and the metSS explained more variability than either the clinical covariate or metSS-only models (53.4% vs. 46.4% and 31.6%) in the validation dataset. The metabolomic pathways for arginine biosynthesis; aminoacyl-tRNA biosynthesis; and glycine, serine, and threonine pathway were enriched by adaLASSO metabolites for FEV₁. This is the first demonstration of a respiratory metabolomic severity score, which shows how a metSS can add explanation of variance to clinical predictors of FEV₁ and MESA-adjusted lung density. The advantage of a comprehensive metSS is that it explains more disease than individual metabolites and can account for substantial collinearity among classes of metabolites. Future studies should be performed to determine whether metSSs are similar in younger, and more racially and ethnically diverse populations as well as whether a metabolomic severity score can predict disease development in individuals who do not yet have COPD.

Metabolomic Profiles in Childhood and Adolescence Are Associated with Fetal Overnutrition

Fetal overnutrition predisposes offspring to increased metabolic risk. The current study used metabolomics to assess sustained differences in serum metabolites across childhood and adolescence among youth exposed to three typologies of fetal overnutrition: maternal obesity only, gestational diabetes mellitus (GDM) only, and obesity + GDM. We included youth exposed in utero to obesity only (BMI ≥ 30; n = 66), GDM only (n = 56), obesity + GDM (n = 25), or unexposed (n = 297), with untargeted metabolomics measured at ages 10 and 16 years. We used linear mixed models to identify metabolites across both time-points associated with exposure to any overnutrition, using a false-discovery-rate correction (FDR) <0.20. These metabolites were included in a principal component analysis (PCA) to generate profiles and assess metabolite profile differences with respect to overnutrition typology (adjusted for prenatal smoking, offspring age, sex, and race/ethnicity). Fetal overnutrition was associated with 52 metabolites. PCA yielded four factors accounting for 17−27% of the variance, depending on age of measurement. We observed differences in three factor patterns with respect to overnutrition typology: sphingomyelin-mannose (8−13% variance), skeletal muscle metabolism (6−10% variance), and 3-carboxy-4-methyl-5-propyl-2-furanpropanoic acid (CMPF; 3−4% variance). The sphingomyelin-mannose factor score was higher among offspring exposed to obesity vs. GDM. Exposure to obesity + GDM (vs. GDM or obesity only) was associated with higher skeletal muscle metabolism and CMPF scores. Fetal overnutrition is associated with metabolic changes in the offspring, but differences between typologies of overnutrition account for a small amount of variation in the metabolome, suggesting there is likely greater pathophysiological overlap than difference.

An Augmented High-Dimensional Graphical Lasso Method to Incorporate Prior Biological Knowledge for Global Network Learning

Biological networks are often inferred through Gaussian graphical models (GGMs) using gene or protein expression data only. GGMs identify conditional dependence by estimating a precision matrix between genes or proteins. However, conventional GGM approaches often ignore prior knowledge about protein-protein interactions (PPI). Recently, several groups have extended GGM to weighted graphical Lasso (wGlasso) and network-based gene set analysis (Netgsa) and have demonstrated the advantages of incorporating PPI information. However, these methods are either computationally intractable for large-scale data, or disregard weights in the PPI networks. To address these shortcomings, we extended the Netgsa approach and developed an augmented high-dimensional graphical Lasso (AhGlasso) method to incorporate edge weights in known PPI with omics data for global network learning. This new method outperforms weighted graphical Lasso-based algorithms with respect to computational time in simulated large-scale data settings while achieving better or comparable prediction accuracy of node connections. The total runtime of AhGlasso is approximately five times faster than weighted Glasso methods when the graph size ranges from 1,000 to 3,000 with a fixed sample size (n = 300). The runtime difference between AhGlasso and weighted Glasso increases when the graph size increases. Using proteomic data from a study on chronic obstructive pulmonary disease, we demonstrate that AhGlasso improves protein network inference compared to the Netgsa approach by incorporating PPI information.

Longitudinal analysis of SARS-CoV-2 spike and RNA-dependent RNA polymerase protein sequences reveals the emergence and geographic distribution of diverse mutations

Amid the ongoing COVID-19 pandemic, it has become increasingly important to monitor the mutations that arise in the SARS-CoV-2 virus, to prepare public health strategies and guide the further development of vaccines and therapeutics. The spike (S) protein and the proteins comprising the RNA-Dependent RNA Polymerase (RdRP) are key vaccine and drug targets, respectively, making mutation surveillance of these proteins of great importance. Full protein sequences were downloaded from the GISAID database, aligned, and the variants identified. 437,006 unique viral genomes were analyzed. Polymorphisms in the protein sequence were investigated and examined longitudinally to identify sequence and strain variants appearing between January 5th, 2020 and January 16th, 2021. A structural analysis was also performed to investigate mutations in the receptor binding domain and the N-terminal domain of the spike protein. Within the spike protein, there were 766 unique mutations observed in the N-terminal domain and 360 in the receptor binding domain. Four residues that directly contact ACE2 were mutated in more than 100 sequences, including positions K417, Y453, S494, and N501. Within the furin cleavage site of the spike protein, a high degree of conservation was observed, but the P681H mutation was observed in 10.47% of sequences analyzed. Within the RNA dependent RNA polymerase complex proteins, 327 unique mutations were observed in Nsp8, 166 unique mutations were observed in Nsp7, and 1157 unique mutations were observed in Nsp12. Only 4 sequences analyzed contained mutations in the 9 residues that directly interact with the therapeutic Remdesivir, suggesting limited mutations in drug interacting residues. The identification of new variants emphasizes the need for further study on the effects of the mutations and the implications of increased prevalence, particularly for vaccine or therapeutic efficacy.

Identifying miRNA-mRNA Networks Associated With COPD Phenotypes

Chronic obstructive pulmonary disease (COPD) is characterized by expiratory airflow limitation and symptoms such as shortness of breath. Although many studies have demonstrated dysregulated microRNA (miRNA) and gene (mRNA) expression in the pathogenesis of COPD, how miRNAs and mRNAs systematically interact and contribute to COPD development is still not clear. To gain a deeper understanding of the gene regulatory network underlying COPD pathogenesis, we used Sparse Multiple Canonical Correlation Network (SmCCNet) to integrate whole blood miRNA and RNA-sequencing data from 404 participants in the COPDGene study to identify novel miRNA-mRNA networks associated with COPD-related phenotypes including lung function and emphysema. We hypothesized that phenotype-directed interpretable miRNA-mRNA networks from SmCCNet would assist in the discovery of novel biomarkers that traditional single biomarker discovery methods (such as differential expression) might fail to discover. Additionally, we investigated whether adjusting -omics and clinical phenotypes data for covariates prior to integration would increase the statistical power for network identification. Our study demonstrated that partial covariate adjustment for age, sex, race, and CT scanner model (in the quantitative emphysema networks) improved network identification when compared with no covariate adjustment. However, further adjustment for current smoking status and relative white blood cell (WBC) proportions sometimes weakened the power for identifying lung function and emphysema networks, a phenomenon which may be due to the correlation of smoking status and WBC counts with the COPD-related phenotypes. With partial covariate adjustment, we found six miRNA-mRNA networks associated with COPD-related phenotypes. One network consists of 2 miRNAs and 28 mRNAs which had a 0.33 correlation (p = 5.40E-12) to forced expiratory volume in 1 s (FEV₁) percent predicted. We also found a network of 5 miRNAs and 81 mRNAs that had a 0.45 correlation (p = 8.80E-22) to percent emphysema. The miRNA-mRNA networks associated with COPD traits provide a systems view of COPD pathogenesis and complements biomarker identification with individual miRNA or mRNA expression data.

PaIRKAT: A pathway integrated regression-based kernel association test with applications to metabolomics and COPD phenotypes

High-throughput data such as metabolomics, genomics, transcriptomics, and proteomics have become familiar data types within the "-omics" family. For this work, we focus on subsets that interact with one another and represent these "pathways" as graphs. Observed pathways often have disjoint components, i.e., nodes or sets of nodes (metabolites, etc.) not connected to any other within the pathway, which notably lessens testing power. In this paper we propose the Pathway Integrated Regression-based Kernel Association Test (PaIRKAT), a new kernel machine regression method for incorporating known pathway information into the semi-parametric kernel regression framework. This work extends previous kernel machine approaches. This paper also contributes an application of a graph kernel regularization method for overcoming disconnected pathways. By incorporating a regularized or "smoothed" graph into a score test, PaIRKAT can provide more powerful tests for associations between biological pathways and phenotypes of interest and will be helpful in identifying novel pathways for targeted clinical research. We evaluate this method through several simulation studies and an application to real metabolomics data from the COPDGene study. Our simulation studies illustrate the robustness of this method to incorrect and incomplete pathway knowledge, and the real data analysis shows meaningful improvements of testing power in pathways. PaIRKAT was developed for application to metabolomic pathway data, but the techniques are easily generalizable to other data sources with a graph-like structure.

An effective processing pipeline for harmonizing DNA methylation data from Illumina's 450K and EPIC platforms for epidemiological studies

OBJECTIVE

Illumina BeadChip arrays are commonly used to generate DNA methylation data for large epidemiological studies. Updates in technology over time create challenges for data harmonization within and between studies, many of which obtained data from the older 450K and newer EPIC platforms. The pre-processing pipeline for DNA methylation is not trivial, and influences the downstream analyses. Incorporating different platforms adds a new level of technical variability that has not yet been taken into account by recommended pipelines. Our study evaluated the performance of various tools on different versions of platform data harmonization at each step of pre-processing pipeline, including quality control (QC), normalization, batch effect adjustment, and genomic inflation. We illustrate our novel approach using 450K and EPIC data from the Diabetes Autoimmunity Study in the Young (DAISY) prospective cohort.

RESULTS

We found normalization and probe filtering had the biggest effect on data harmonization. Employing a meta-analysis was an effective and easily executable method for accounting for platform variability. Correcting for genomic inflation also helped with harmonization. We present guidelines for studies seeking to harmonize data from the 450K and EPIC platforms, which includes the use of technical replicates for evaluating numerous pre-processing steps, and employing a meta-analysis.

Reproducibility of mass spectrometry based metabolomics data

BACKGROUND

Assessing the reproducibility of measurements is an important first step for improving the reliability of downstream analyses of high-throughput metabolomics experiments. We define a metabolite to be reproducible when it demonstrates consistency across replicate experiments. Similarly, metabolites which are not consistent across replicates can be labeled as irreproducible. In this work, we introduce and evaluate the use (Ma)ximum (R)ank (R)eproducibility (MaRR) to examine reproducibility in mass spectrometry-based metabolomics experiments. We examine reproducibility across technical or biological samples in three different mass spectrometry metabolomics (MS-Metabolomics) data sets.

RESULTS

We apply MaRR, a nonparametric approach that detects the change from reproducible to irreproducible signals using a maximal rank statistic. The advantage of using MaRR over model-based methods that it does not make parametric assumptions on the underlying distributions or dependence structures of reproducible metabolites. Using three MS Metabolomics data sets generated in the multi-center Genetic Epidemiology of Chronic Obstructive Pulmonary Disease (COPD) study, we applied the MaRR procedure after data processing to explore reproducibility across technical or biological samples. Under realistic settings of MS-Metabolomics data, the MaRR procedure effectively controls the False Discovery Rate (FDR) when there was a gradual reduction in correlation between replicate pairs for less highly ranked signals. Simulation studies also show that the MaRR procedure tends to have high power for detecting reproducible metabolites in most situations except for smaller values of proportion of reproducible metabolites. Bias (i.e., the difference between the estimated and the true value of reproducible signal proportions) values for simulations are also close to zero. The results reported from the real data show a higher level of reproducibility for technical replicates compared to biological replicates across all the three different datasets. In summary, we demonstrate that the MaRR procedure application can be adapted to various experimental designs, and that the nonparametric approach performs consistently well.

CONCLUSIONS

This research was motivated by reproducibility, which has proven to be a major obstacle in the use of genomic findings to advance clinical practice. In this paper, we developed a data-driven approach to assess the reproducibility of MS-Metabolomics data sets. The methods described in this paper are implemented in the open-source R package marr, which is freely available from Bioconductor at http://bioconductor.org/packages/marr .

A Mediation Approach to Discovering Causal Relationships between the Metabolome and DNA Methylation in Type 1 Diabetes

Environmental factors including viruses, diet, and the metabolome have been linked with the appearance of islet autoimmunity (IA) that precedes development of type 1 diabetes (T1D). We measured global DNA methylation (DNAm) and untargeted metabolomics prior to IA and at the time of seroconversion to IA in 92 IA cases and 91 controls from the Diabetes Autoimmunity Study in the Young (DAISY). Causal mediation models were used to identify seven DNAm probe-metabolite pairs in which the metabolite measured at IA mediated the protective effect of the DNAm probe measured prior to IA against IA risk. These pairs included five DNAm probes mediated by histidine (a metabolite known to affect T1D risk), one probe (cg01604946) mediated by phostidyl choline p-32:0 or o-32:1, and one probe (cg00390143) mediated by sphingomyelin d34:2. The top 100 DNAm probes were over-represented in six reactome pathways at the FDR <0.1 level (q = 0.071), including transport of small molecules and inositol phosphate metabolism. While the causal pathways in our mediation models require further investigation to better understand the biological mechanisms, we identified seven methylation sites that may improve our understanding of epigenetic protection against T1D as mediated by the metabolome.

The oxylipin profile is associated with development of type 1 diabetes: the Diabetes Autoimmunity Study in the Young (DAISY)

AIMS/HYPOTHESIS

Oxylipins are lipid mediators derived from polyunsaturated fatty acids. Some oxylipins are proinflammatory (e.g. those derived from arachidonic acid [ARA]), others are pro-resolving of inflammation (e.g. those derived from α-linolenic acid [ALA], docosahexaenoic acid [DHA] and eicosapentaenoic acid [EPA]) and others may be both (e.g. those derived from linoleic acid [LA]). The goal of this study was to examine whether oxylipins are associated with incident type 1 diabetes.

METHODS

We conducted a nested case-control analysis in the Diabetes Autoimmunity Study in the Young (DAISY), a prospective cohort study of children at risk of type 1 diabetes. Plasma levels of 14 ARA-derived oxylipins, ten LA-derived oxylipins, six ALA-derived oxylipins, four DHA-derived oxylipins and two EPA-related oxylipins were measured by ultra-HPLC-MS/MS at multiple timepoints related to autoantibody seroconversion in 72 type 1 diabetes cases and 71 control participants, which were frequency matched on age at autoantibody seroconversion (of the case), ethnicity and sample availability. Linear mixed models were used to obtain an age-adjusted mean of each oxylipin prior to type 1 diabetes. Age-adjusted mean oxylipins were tested for association with type 1 diabetes using logistic regression, adjusting for the high risk HLA genotype HLA-DR3/4,DQB1*0302. We also performed principal component analysis of the oxylipins and tested principal components (PCs) for association with type 1 diabetes. Finally, to investigate potential critical timepoints, we examined the association of oxylipins measured before and after autoantibody seroconversion (of the cases) using PCs of the oxylipins at those visits.

RESULTS

The ARA-related oxylipin 5-HETE was associated with increased type 1 diabetes risk. Five LA-related oxylipins, two ALA-related oxylipins and one DHA-related oxylipin were associated with decreased type 1 diabetes risk. A profile of elevated LA- and ALA-related oxylipins (PC1) was associated with decreased type 1 diabetes risk (OR 0.61; 95% CI 0.40, 0.94). A profile of elevated ARA-related oxylipins (PC2) was associated with increased diabetes risk (OR 1.53; 95% CI 1.03, 2.29). A critical timepoint analysis showed type 1 diabetes was associated with a high ARA-related oxylipin profile at post-autoantibody-seroconversion but not pre-seroconversion.

CONCLUSIONS/INTERPRETATION

The protective association of higher LA- and ALA-related oxylipins demonstrates the importance of both inflammation promotion and resolution in type 1 diabetes. Proinflammatory ARA-related oxylipins may play an important role once the autoimmune process has begun.

Phospholipid Levels at Seroconversion Are Associated With Resolution of Persistent Islet Autoimmunity: The Diabetes Autoimmunity Study in the Young

Reversion of islet autoimmunity (IA) may point to mechanisms that prevent IA progression. We followed 199 individuals who developed IA during the Diabetes Autoimmunity Study in the Young. Untargeted metabolomics was performed in serum samples following IA. Cox proportional hazards models were used to test whether the metabolites (2,487) predicted IA reversion: two or more consecutive visits negative for all autoantibodies. We conducted a principal components analysis (PCA) of the top metabolites; |hazard ratio (HR) >1.25| and nominal P < 0.01. Phosphatidylcholine (16:0_18:1(9Z)) was the strongest individual metabolite (HR per 1 SD 2.16, false discovery rate (FDR)-adjusted P = 0.0037). Enrichment analysis identified four clusters (FDR P < 0.10) characterized by an overabundance of sphingomyelin (d40:0), phosphatidylcholine (16:0_18:1(9Z)), phosphatidylcholine (30:0), and l-decanoylcarnitine. Overall, 63 metabolites met the criteria for inclusion in the PCA. PC1 (HR 1.4, P < 0.0001), PC2 (HR 0.85, P = 0.0185), and PC4 (HR 1.28, P = 0.0103) were associated with IA reversion. Given the potential influence of diet on the metabolome, we investigated whether nutrients were correlated with PCs. We identified 20 nutrients that were correlated with the PCs (P < 0.05). Total sugar intake was the top nutrient. Overall, we identified an association between phosphatidylcholine, sphingomyelin, and carnitine levels and reversion of IA.

Human Mobility Associated With Risk of Schistosoma japonicum Infection in Sichuan, China

Urbanization increases human mobility in ways that can alter the transmission of classically rural, vector-borne diseases like schistosomiasis. The impact of human mobility on individual-level Schistosoma risk is poorly characterized. Travel outside endemic areas may protect against infection by reducing exposure opportunities, whereas travel to other endemic regions may increase risk. Using detailed monthly travel- and water-contact surveys from 27 rural communities in Sichuan, China, in 2008, we aimed to describe human mobility and to identify mobility-related predictors of S. japonicum infection. Candidate predictors included timing, frequency, distance, duration, and purpose of recent travel as well as water-contact measures. Random forests machine learning was used to detect key predictors of individual infection status. Logistic regression was used to assess the strength and direction of associations. Key mobility-related predictors include frequent travel and travel during July-both associated with decreased probability of infection and less time engaged in risky water-contact behavior, suggesting travel may remove opportunities for schistosome exposure. The importance of July travel and July water contact suggests a high-risk window for cercarial exposure. The frequency and timing of human movement out of endemic areas should be considered when assessing potential drivers of rural infectious diseases.

Inverse probability weighting is an effective method to address selection bias during the analysis of high dimensional data

Omics studies frequently use samples collected during cohort studies. Conditioning on sample availability can cause selection bias if sample availability is nonrandom. Inverse probability weighting (IPW) is purported to reduce this bias. We evaluated IPW in an epigenome-wide analysis testing the association between DNA methylation (261,435 probes) and age in healthy adolescent subjects (n = 114). We simulated age and sex to be correlated with sample selection and then evaluated four conditions: complete population/no selection bias (all subjects), naïve selection bias (no adjustment), and IPW selection bias (selection bias with IPW adjustment). Assuming the complete population condition represented the "truth," we compared each condition to the complete population condition. Bias or difference in associations between age and methylation was reduced in the IPW condition versus the naïve condition. However, genomic inflation and type 1 error were higher in the IPW condition relative to the naïve condition. Postadjustment using bacon, type 1 error and inflation were similar across all conditions. Power was higher under the IPW condition compared with the naïve condition before and after inflation adjustment. IPW methods can reduce bias in genome-wide analyses. Genomic inflation is a potential concern that can be minimized using methods that adjust for inflation.

Comparing Statistical Tests for Differential Network Analysis of Gene Modules

Genes often work together to perform complex biological processes, and "networks" provide a versatile framework for representing the interactions between multiple genes. Differential network analysis (DiNA) quantifies how this network structure differs between two or more groups/phenotypes (e.g., disease subjects and healthy controls), with the goal of determining whether differences in network structure can help explain differences between phenotypes. In this paper, we focus on gene co-expression networks, although in principle, the methods studied can be used for DiNA for other types of features (e.g., metabolome, epigenome, microbiome, proteome, etc.). Three common applications of DiNA involve (1) testing whether the connections to a single gene differ between groups, (2) testing whether the connection between a pair of genes differs between groups, or (3) testing whether the connections within a "module" (a subset of 3 or more genes) differs between groups. This article focuses on the latter, as there is a lack of studies comparing statistical methods for identifying differentially co-expressed modules (DCMs). Through extensive simulations, we compare several previously proposed test statistics and a new p-norm difference test (PND). We demonstrate that the true positive rate of the proposed PND test is competitive with and often higher than the other methods, while controlling the false positive rate. The R package discoMod (differentially co-expressed modules) implements the proposed method and provides a full pipeline for identifying DCMs: clustering tools to derive gene modules, tests to identify DCMs, and methods for visualizing the results.

Soluble receptor for advanced glycation end products (sRAGE) as a biomarker of COPD

BACKGROUND

Soluble receptor for advanced glycation end products (sRAGE) is a proposed emphysema and airflow obstruction biomarker; however, previous publications have shown inconsistent associations and only one study has investigate the association between sRAGE and emphysema. No cohorts have examined the association between sRAGE and progressive decline of lung function. There have also been no evaluation of assay compatibility, receiver operating characteristics, and little examination of the effect of genetic variability in non-white population. This manuscript addresses these deficiencies and introduces novel data from Pittsburgh COPD SCCOR and as well as novel work on airflow obstruction. A meta-analysis is used to quantify sRAGE associations with clinical phenotypes.

METHODS

sRAGE was measured in four independent longitudinal cohorts on different analytic assays: COPDGene (n = 1443); SPIROMICS (n = 1623); ECLIPSE (n = 2349); Pittsburgh COPD SCCOR (n = 399). We constructed adjusted linear mixed models to determine associations of sRAGE with baseline and follow up forced expiratory volume at one second (FEV₁) and emphysema by quantitative high-resolution CT lung density at the 15th percentile (adjusted for total lung capacity).

RESULTS

Lower plasma or serum sRAGE values were associated with a COPD diagnosis (P < 0.001), reduced FEV₁ (P < 0.001), and emphysema severity (P < 0.001). In an inverse-variance weighted meta-analysis, one SD lower log₁₀-transformed sRAGE was associated with 105 ± 22 mL lower FEV₁ and 4.14 ± 0.55 g/L lower adjusted lung density. After adjusting for covariates, lower sRAGE at baseline was associated with greater FEV₁ decline and emphysema progression only in the ECLIPSE cohort. Non-Hispanic white subjects carrying the rs2070600 minor allele (A) and non-Hispanic African Americans carrying the rs2071288 minor allele (A) had lower sRAGE measurements compare to those with the major allele, but their emphysema-sRAGE regression slopes were similar.

CONCLUSIONS

Lower blood sRAGE is associated with more severe airflow obstruction and emphysema, but associations with progression are inconsistent in the cohorts analyzed. In these cohorts, genotype influenced sRAGE measurements and strengthened variance modelling. Thus, genotype should be included in sRAGE evaluations.

tidyMicro: a pipeline for microbiome data analysis and visualization using the tidyverse in R

Background

The drive to understand how microbial communities interact with their environments has inspired innovations across many fields. The data generated from sequence-based analyses of microbial communities typically are of high dimensionality and can involve multiple data tables consisting of taxonomic or functional gene/pathway counts. Merging multiple high dimensional tables with study-related metadata can be challenging. Existing microbiome pipelines available in R have created their own data structures to manage this problem. However, these data structures may be unfamiliar to analysts new to microbiome data or R and do not allow for deviations from internal workflows. Existing analysis tools also focus primarily on community-level analyses and exploratory visualizations, as opposed to analyses of individual taxa.

Results

We developed the R package “tidyMicro” to serve as a more complete microbiome analysis pipeline. This open source software provides all of the essential tools available in other popular packages (e.g., management of sequence count tables, standard exploratory visualizations, and diversity inference tools) supplemented with multiple options for regression modelling (e.g., negative binomial, beta binomial, and/or rank based testing) and novel visualizations to improve interpretability (e.g., Rocky Mountain plots, longitudinal ordination plots). This comprehensive pipeline for microbiome analysis also maintains data structures familiar to R users to improve analysts’ control over workflow. A complete vignette is provided to aid new users in analysis workflow.

Conclusions

tidyMicro provides a reliable alternative to popular microbiome analysis packages in R. We provide standard tools as well as novel extensions on standard analyses to improve interpretability results while maintaining object malleability to encourage open source collaboration. The simple examples and full workflow from the package are reproducible and applicable to external data sets.

A Practical Guide to Metabolomics Software Development

A growing number of software tools have been developed for metabolomics data processing and analysis. Many new tools are contributed by metabolomics practitioners who have limited prior experience with software development, and the tools are subsequently implemented by users with expertise that ranges from basic point-and-click data analysis to advanced coding. This Perspective is intended to introduce metabolomics software users and developers to important considerations that determine the overall impact of a publicly available tool within the scientific community. The recommendations reflect the collective experience of an NIH-sponsored Metabolomics Consortium working group that was formed with the goal of researching guidelines and best practices for metabolomics tool development. The recommendations are aimed at metabolomics researchers with little formal background in programming and are organized into three stages: (i) preparation, (ii) tool development, and (iii) distribution and maintenance.

Predictive Modeling for Metabolomics Data

In recent years, mass spectrometry (MS)-based metabolomics has been extensively applied to characterize biochemical mechanisms, and study physiological processes and phenotypic changes associated with disease. Metabolomics has also been important for identifying biomarkers of interest suitable for clinical diagnosis. For the purpose of predictive modeling, in this chapter, we will review various supervised learning algorithms such as random forest (RF), support vector machine (SVM), and partial least squares-discriminant analysis (PLS-DA). In addition, we will also review feature selection methods for identifying the best combination of metabolites for an accurate predictive model. We conclude with best practices for reproducibility by including internal and external replication, reporting metrics to assess performance, and providing guidelines to avoid overfitting and to deal with imbalanced classes. An analysis of an example data will illustrate the use of different machine learning methods and performance metrics.

Prenatal Exposure to Per- and Polyfluoroalkyl Substances, Umbilical Cord Blood DNA Methylation, and Cardio-Metabolic Indicators in Newborns: The Healthy Start Study

BACKGROUND

Per- and polyfluoroalkyl substances (PFAS) are environmentally persistent chemicals widely detected in women of reproductive age. Prenatal PFAS exposure is associated with adverse health outcomes in children. We hypothesized that DNA methylation changes may result from prenatal PFAS exposure and may be linked to offspring cardio-metabolic phenotype.

OBJECTIVES

We estimated associations of prenatal PFAS with DNA methylation in umbilical cord blood. We evaluated associations of methylation at selected sites with neonatal cardio-metabolic indicators.

METHODS

Among 583 mother-infant pairs in a prospective cohort, five PFAS were quantified in maternal serum (median 27 wk of gestation). Umbilical cord blood DNA methylation was evaluated using the Illumina HumanMethylation450 array. Differentially methylated positions (DMPs) were evaluated at a false discovery rate ( FDR ) < 0.05 and differentially methylated regions (DMRs) were identified using comb-p (Šidák-adjusted p < 0.05 ). We estimated associations between methylation at candidate DMPs and DMR sites and the following outcomes: newborn weight, adiposity, and cord blood glucose, insulin, lipids, and leptin.

RESULTS

Maternal serum PFAS concentrations were below the median for females in the U.S. general population. Moderate to high pairwise correlations were observed between PFAS concentrations (ρ = 0.28 - 0.76 ). Methylation at one DMP (cg18587484), annotated to the gene TJAP1, was associated with perfluorooctanoate (PFOA) at FDR <   0.05 . Comb-p detected between 4 and 15 DMRs for each PFAS. Associated genes, some common across multiple PFAS, were implicated in growth (RPTOR), lipid homeostasis (PON1, PON3, CIDEB, NR1H2), inflammation and immune activity (RASL11B, RNF39), among other functions. There was suggestive evidence that two PFAS-associated loci (cg09093485, cg09637273) were associated with cord blood triglycerides and birth weight, respectively (FDR <   0.1 ).

DISCUSSION

DNA methylation in umbilical cord blood was associated with maternal serum PFAS concentrations during pregnancy, suggesting potential associations with offspring growth, metabolism, and immune function. Future research should explore whether DNA methylation changes mediate associations between prenatal PFAS exposures and child health outcomes. https://doi.org/10.1289/EHP6888.

Associations between the activity of placental nutrient-sensing pathways and neonatal and postnatal metabolic health: the ECHO Healthy Start cohort

OBJECTIVE

Our hypothesis was that the activity of placental nutrient-sensing pathways is associated with adiposity and metabolic health in childhood.

RESEARCH DESIGN AND METHODS

Using placental villus samples from healthy mothers from the Healthy Start Study, we measured the abundance and phosphorylation of key intermediates in the mTOR, insulin, AMPK, and ER stress signaling pathways. Using multivariate multiple regression models, we tested the association between placental proteins and offspring adiposity (%fat mass) at birth (n = 109), 4-6 months (n = 104), and 4-6 years old (n = 64), adjusted for offspring sex and age.

RESULTS

Placental mTORC1 phosphorylation was positively associated with adiposity at birth (R2 = 0.13, P = 0.009) and 4-6 years (R2 = 0.15, P = 0.046). The mTORC2 target PKCα was positively associated with systolic blood pressure at 4-6 years (β = 2.90, P = 0.005). AMPK phosphorylation was positively associated with adiposity at birth (β = 2.32, P = 0.023), but the ratio of phosphorylated to total AMPK was negatively associated with skinfold thickness (β = -2.37, P = 0.022) and body weight (β = -2.92, P = 0.005) at 4-6 years.

CONCLUSIONS

This is the first report of associations between key placental protein activity measures and longitudinal child outcomes at various life stages. Our data indicate that AMPK and mTOR signaling are linked to cardiometabolic measures at birth and 4-6 years, providing novel insight into potential mechanisms underpinning how metabolic signaling in the placenta is associated with future risk of cardiovascular disease.

Metabolomics-related nutrient patterns at seroconversion and risk of progression to type 1 diabetes

OBJECTIVE

Our aim was to elucidate the role of diet in type 1 diabetes (T1D) by examining combinations of nutrient intake in the progression from islet autoimmunity (IA) to T1D.

METHODS

We measured 2457 metabolites and dietary intake at the time of seroconversion in 132 IA-positive children in the prospective Diabetes Autoimmunity Study in the Young. IA was defined as the first of two consecutive visits positive for at least one autoantibody (insulin, GAD, IA-2, or ZnT8). By December 2018, 40 children progressed to T1D. Intakes of 38 nutrients were estimated from semiquantitative food frequency questionnaires. We tested the association of each metabolite with progression to T1D using multivariable Cox regression. Nutrient patterns that best explained variation in candidate metabolites were identified using reduced rank regression (RRR), and their association with progression to T1D was tested using Cox regression adjusting for age at seroconversion and high-risk HLA genotype.

RESULTS

In stepwise selection, 22 nutrients significantly predicted at least two of the 13 most significant metabolites associated with progression to T1D, and were included in RRR. A nutrient pattern corresponding to intake lower in linoleic acid, niacin, and riboflavin, and higher in total sugars, explained 18% of metabolite variability. Children scoring higher on this metabolite-related nutrient pattern at seroconversion had increased risk for progressing to T1D (HR = 3.17, 95%CI = 1.42-7.05).

CONCLUSIONS

Combinations of nutrient intake reflecting candidate metabolites are associated with increased risk of T1D, and may help focus dietary prevention efforts.

Systems genetics analysis of the LXS recombinant inbred mouse strains:Genetic and molecular insights into acute ethanol tolerance

We have been using the Inbred Long- and Short-Sleep mouse strains (ILS, ISS) and a recombinant inbred panel derived from them, the LXS, to investigate the genetic underpinnings of acute ethanol tolerance which is considered to be a risk factor for alcohol use disorders (AUDs). Here, we have used RNA-seq to examine the transcriptome of whole brain in 40 of the LXS strains 8 hours after a saline or ethanol “pretreatment” as in previous behavioral studies. Approximately 1/3 of the 14,184 expressed genes were significantly heritable and many were unique to the pretreatment. Several thousand cis- and trans-eQTLs were mapped; a portion of these also were unique to pretreatment. Ethanol pretreatment caused differential expression (DE) of 1,230 genes. Gene Ontology (GO) enrichment analysis suggested involvement in numerous biological processes including astrocyte differentiation, histone acetylation, mRNA splicing, and neuron projection development. Genetic correlation analysis identified hundreds of genes that were correlated to the behaviors. GO analysis indicated that these genes are involved in gene expression, chromosome organization, and protein transport, among others. The expression profiles of the DE genes and genes correlated to AFT in the ethanol pretreatment group (AFT-Et) were found to be similar to profiles of HDAC inhibitors. Hdac1, a cis-regulated gene that is located at the peak of a previously mapped QTL for AFT-Et, was correlated to 437 genes, most of which were also correlated to AFT-Et. GO analysis of these genes identified several enriched biological process terms including neuron-neuron synaptic transmission and potassium transport. In summary, the results suggest widespread genetic effects on gene expression, including effects that are pretreatment-specific. A number of candidate genes and biological functions were identified that could be mediating the behavioral responses. The most prominent of these was Hdac1 which may be regulating genes associated with glutamatergic signaling and potassium conductance.

Qualitative Differences Between the IFNα subtypes and IFNβ Influence Chronic Mucosal HIV-1 Pathogenesis

The Type I Interferons (IFN-Is) are innate antiviral cytokines that include 12 different IFNα subtypes and IFNβ that signal through the IFN-I receptor (IFNAR), inducing hundreds of IFN-stimulated genes (ISGs) that comprise the 'interferome'. Quantitative differences in IFNAR binding correlate with antiviral activity, but whether IFN-Is exhibit qualitative differences remains controversial. Moreover, the IFN-I response is protective during acute HIV-1 infection, but likely pathogenic during the chronic stages. To gain a deeper understanding of the IFN-I response, we compared the interferomes of IFNα subtypes dominantly-expressed in HIV-1-exposed plasmacytoid dendritic cells (1, 2, 5, 8 and 14) and IFNβ in the earliest cellular targets of HIV-1 infection. Primary gut CD4 T cells from 3 donors were treated for 18 hours ex vivo with individual IFN-Is normalized for IFNAR signaling strength. Of 1,969 IFN-regulated genes, 246 'core ISGs' were induced by all IFN-Is tested. However, many IFN-regulated genes were not shared between the IFNα subtypes despite similar induction of canonical antiviral ISGs such as ISG15, RSAD2 and MX1, formally demonstrating qualitative differences between the IFNα subtypes. Notably, IFNβ induced a broader interferome than the individual IFNα subtypes. Since IFNβ, and not IFNα, is upregulated during chronic HIV-1 infection in the gut, we compared core ISGs and IFNβ-specific ISGs from colon pinch biopsies of HIV-1-uninfected (n = 13) versus age- and gender-matched, antiretroviral-therapy naïve persons with HIV-1 (PWH; n = 19). Core ISGs linked to inflammation, T cell activation and immune exhaustion were elevated in PWH, positively correlated with plasma lipopolysaccharide (LPS) levels and gut IFNβ levels, and negatively correlated with gut CD4 T cell frequencies. In sharp contrast, IFNβ-specific ISGs linked to protein translation and anti-inflammatory responses were significantly downregulated in PWH, negatively correlated with gut IFNβ and LPS, and positively correlated with plasma IL6 and gut CD4 T cell frequencies. Our findings reveal qualitative differences in interferome induction by diverse IFN-Is and suggest potential mechanisms for how IFNβ may drive HIV-1 pathogenesis in the gut.

MCMSeq: Bayesian hierarchical modeling of clustered and repeated measures RNA sequencing experiments

BACKGROUND

As the barriers to incorporating RNA sequencing (RNA-Seq) into biomedical studies continue to decrease, the complexity and size of RNA-Seq experiments are rapidly growing. Paired, longitudinal, and other correlated designs are becoming commonplace, and these studies offer immense potential for understanding how transcriptional changes within an individual over time differ depending on treatment or environmental conditions. While several methods have been proposed for dealing with repeated measures within RNA-Seq analyses, they are either restricted to handling only paired measurements, can only test for differences between two groups, and/or have issues with maintaining nominal false positive and false discovery rates. In this work, we propose a Bayesian hierarchical negative binomial generalized linear mixed model framework that can flexibly model RNA-Seq counts from studies with arbitrarily many repeated observations, can include covariates, and also maintains nominal false positive and false discovery rates in its posterior inference.

RESULTS

In simulation studies, we showed that our proposed method (MCMSeq) best combines high statistical power (i.e. sensitivity or recall) with maintenance of nominal false positive and false discovery rates compared the other available strategies, especially at the smaller sample sizes investigated. This behavior was then replicated in an application to real RNA-Seq data where MCMSeq was able to find previously reported genes associated with tuberculosis infection in a cohort with longitudinal measurements.

CONCLUSIONS

Failing to account for repeated measurements when analyzing RNA-Seq experiments can result in significantly inflated false positive and false discovery rates. Of the methods we investigated, whether they model RNA-Seq counts directly or worked on transformed values, the Bayesian hierarchical model implemented in the mcmseq R package (available at https://github.com/stop-pre16/mcmseq ) best combined sensitivity and nominal error rate control.