Bayesian modeling of factorial time-course data with applications to a bone aging gene expression study.

Many scientific studies, especially in the biomedical sciences, generate data measured simultaneously over a multitude of units, over a period of time, and under different conditions or combinations of factors. Often, an important question of interest asked relates to which units behave similarly under different conditions, but measuring the variation over time complicates the analysis significantly. In this article we address such a problem arising from a gene expression study relating to bone aging, and develop a Bayesian statistical method that can simultaneously detect and uncover signals on three levels within such data: factorial, longitudinal, and transcriptional. Our model framework considers both cluster and time-point-specific parameters and these parameters uniquely determine the shapes of the temporal gene expression profiles, allowing the discovery and characterization of latent gene clusters based on similar underlying biological mechanisms. Our methodology was successfully applied to discover transcriptional networks in a microarray data set comparing the transcriptomic changes that occurred during bone aging in male and female mice expressing one or both copies of the bromodomain (Brd2) gene, a transcriptional regulator which exhibits an age-dependent sex-linked bone loss phenotype.

The impact of prior information on estimates of disease transmissibility using Bayesian tools.

The basic reproductive number (R0) and the distribution of the serial interval (SI) are often used to quantify transmission during an infectious disease outbreak. In this paper, we present estimates of R0 and SI from the 2003 SARS outbreak in Hong Kong and Singapore, and the 2009 pandemic influenza A(H1N1) outbreak in South Africa using methods that expand upon an existing Bayesian framework. This expanded framework allows for the incorporation of additional information, such as contact tracing or household data, through prior distributions. The results for the R0 and the SI from the influenza outbreak in South Africa were similar regardless of the prior information (R0 = 1.36-1.46, µ = 2.0-2.7, µ = mean of the SI). The estimates of R0 and µ for the SARS outbreak ranged from 2.0-4.4 and 7.4-11.3, respectively, and were shown to vary depending on the use of contact tracing data. The impact of the contact tracing data was likely due to the small number of SARS cases relative to the size of the contact tracing sample.

Strategies to design and analyze targeted sequencing data: cohorts for Heart and Aging Research in Genomic Epidemiology (CHARGE) Consortium Targeted Sequencing Study.

BACKGROUND: Genome-wide association studies have identified thousands of genetic variants that influence a variety of diseases and health-related quantitative traits. However, the causal variants underlying the majority of genetic associations remain unknown. Cohorts for Heart and Aging Research in Genomic Epidemiology (CHARGE) Consortium Targeted Sequencing Study aims to follow up genome-wide association study signals and identify novel associations of the allelic spectrum of identified variants with cardiovascular-related traits. METHODS AND RESULTS: The study included 4231 participants from 3 CHARGE cohorts: the Atherosclerosis Risk in Communities Study, the Cardiovascular Health Study, and the Framingham Heart Study. We used a case-cohort design in which we selected both a random sample of participants and participants with extreme phenotypes for each of 14 traits. We sequenced and analyzed 77 genomic loci, which had previously been associated with ≥1 of 14 phenotypes. A total of 52 736 variants were characterized by sequencing and passed our stringent quality control criteria. For common variants (minor allele frequency ≥1%), we performed unweighted regression analyses to obtain P values for associations and weighted regression analyses to obtain effect estimates that accounted for the sampling design. For rare variants, we applied 2 approaches: collapsed aggregate statistics and joint analysis of variants using the sequence kernel association test. CONCLUSIONS: We sequenced 77 genomic loci in participants from 3 cohorts. We established a set of filters to identify high-quality variants and implemented statistical and bioinformatics strategies to analyze the sequence data and identify potentially functional variants within genome-wide association study loci.

Associations of NINJ2 sequence variants with incident ischemic stroke in the Cohorts for Heart and Aging in Genomic Epidemiology (CHARGE) consortium.

BACKGROUND: Stroke, the leading neurologic cause of death and disability, has a substantial genetic component. We previously conducted a genome-wide association study (GWAS) in four prospective studies from the Cohorts for Heart and Aging Research in Genomic Epidemiology (CHARGE) consortium and demonstrated that sequence variants near the NINJ2 gene are associated with incident ischemic stroke. Here, we sought to fine-map functional variants in the region and evaluate the contribution of rare variants to ischemic stroke risk. METHODS AND RESULTS: We sequenced 196 kb around NINJ2 on chromosome 12p13 among 3,986 European ancestry participants, including 475 ischemic stroke cases, from the Atherosclerosis Risk in Communities Study, Cardiovascular Health Study, and Framingham Heart Study. Meta-analyses of single-variant tests for 425 common variants (minor allele frequency [MAF] ≥ 1%) confirmed the original GWAS results and identified an independent intronic variant, rs34166160 (MAF = 0.012), most significantly associated with incident ischemic stroke (HR = 1.80, p = 0.0003). Aggregating 278 putatively-functional variants with MAF≤ 1% using count statistics, we observed a nominally statistically significant association, with the burden of rare NINJ2 variants contributing to decreased ischemic stroke incidence (HR = 0.81; p = 0.026). CONCLUSION: Common and rare variants in the NINJ2 region were nominally associated with incident ischemic stroke among a subset of CHARGE participants. Allelic heterogeneity at this locus, caused by multiple rare, low frequency, and common variants with disparate effects on risk, may explain the difficulties in replicating the original GWAS results. Additional studies that take into account the complex allelic architecture at this locus are needed to confirm these findings.

Targeted sequencing in candidate genes for atrial fibrillation: the Cohorts for Heart and Aging Research in Genomic Epidemiology (CHARGE) Targeted Sequencing Study.

BACKGROUND: Genome-wide association studies (GWAS) have identified common genetic variants that predispose to atrial fibrillation (AF). It is unclear whether rare and low-frequency variants in genes implicated by such GWAS confer additional risk of AF. OBJECTIVE: To study the association of genetic variants with AF at GWAS top loci. METHODS: In the Cohorts for Heart and Aging Research in Genomic Epidemiology (CHARGE) Targeted Sequencing Study, we selected and sequenced 77 target gene regions from GWAS loci of complex diseases or traits, including 4 genes hypothesized to be related to AF (PRRX1, CAV1, CAV2, and ZFHX3). Sequencing was performed in participants with (n = 948) and without (n = 3330) AF from the Atherosclerosis Risk in Communities Study, the Cardiovascular Health Study, the Framingham Heart Study, and the Massachusetts General Hospital. RESULTS: One common variant (rs11265611; P = 1.70 × 10(-6)) intronic to IL6R (interleukin-6 receptor gene) was significantly associated with AF after Bonferroni correction (odds ratio 0.70; 95% confidence interval 0.58-0.85). The variant was not genotyped or imputed by prior GWAS, but it is in linkage disequilibrium (r(2) = .69) with the single-nucleotide polymorphism, with the strongest association with AF so far at this locus (rs4845625). In the rare variant joint analysis, damaging variants within the PRRX1 region showed significant association with AF after Bonferroni correction (P = .01). CONCLUSIONS: We identified 1 common single-nucleotide polymorphism and 1 gene region that were significantly associated with AF. Future sequencing efforts with larger sample sizes and more comprehensive genome coverage are anticipated to identify additional AF-related variants.

Differential expression analysis with global network adjustment.

BACKGROUND: Large-scale chromosomal deletions or other non-specific perturbations of the transcriptome can alter the expression of hundreds or thousands of genes, and it is of biological interest to understand which genes are most profoundly affected. We present a method for predicting a gene's expression as a function of other genes thereby accounting for the effect of transcriptional regulation that confounds the identification of genes differentially expressed relative to a regulatory network. The challenge in constructing such models is that the number of possible regulator transcripts within a global network is on the order of thousands, and the number of biological samples is typically on the order of 10. Nevertheless, there are large gene expression databases that can be used to construct networks that could be helpful in modeling transcriptional regulation in smaller experiments. RESULTS: We demonstrate a type of penalized regression model that can be estimated from large gene expression databases, and then applied to smaller experiments. The ridge parameter is selected by minimizing the cross-validation error of the predictions in the independent out-sample. This tends to increase the model stability and leads to a much greater degree of parameter shrinkage, but the resulting biased estimation is mitigated by a second round of regression. Nevertheless, the proposed computationally efficient "over-shrinkage" method outperforms previously used LASSO-based techniques. In two independent datasets, we find that the median proportion of explained variability in expression is approximately 25%, and this results in a substantial increase in the signal-to-noise ratio allowing more powerful inferences on differential gene expression leading to biologically intuitive findings. We also show that a large proportion of gene dependencies are conditional on the biological state, which would be impossible with standard differential expression methods. CONCLUSIONS: By adjusting for the effects of the global network on individual genes, both the sensitivity and reliability of differential expression measures are greatly improved.

A comparison of gene region simulation methods.

BACKGROUND: Accurately modeling LD in simulations is essential to correctly evaluate new and existing association methods. At present, there has been minimal research comparing the quality of existing gene region simulation methods to produce LD structures similar to an existing gene region. Here we compare the ability of three approaches to accurately simulate the LD within a gene region: HapSim (2005), Hapgen (2009), and a minor extension to simple haplotype resampling. METHODOLOGY/PRINCIPAL FINDINGS: In order to observe the variation and bias for each method, we compare the simulated pairwise LD measures and minor allele frequencies to the original HapMap data in an extensive simulation study. When possible, we also evaluate the effects of changing parameters. HapSim produces samples of haplotypes with lower LD, on average, compared to the original haplotype set while both our resampling method and Hapgen do not introduce this bias. The variation introduced across the replicates by our resampling method is quite small and may not provide enough sampling variability to make a generalizable simulation study. CONCLUSION: We recommend using Hapgen to simulate replicate haplotypes from a gene region. Hapgen produces moderate sampling variation between the replicates while retaining the overall unique LD structure of the gene region.

A generalized hidden Markov model for determining sequence-based predictors of nucleosome positioning.

Chromatin structure, in terms of positioning of nucleosomes and nucleosome-free regions in the DNA, has been found to have an immense impact on various cell functions and processes, ranging from transcriptional regulation to growth and development. In spite of numerous experimental and computational approaches being developed in the past few years to determine the intrinsic relationship between chromatin structure (nucleosome positioning) and DNA sequence features, there is yet no universally accurate approach to predict nucleosome positioning from the underlying DNA sequence alone. We here propose an alternative approach to predicting nucleosome positioning from sequence, making use of characteristic sequence differences, and inherent dependencies in overlapping sequence features. Our nucleosomal positioning prediction algorithm, based on the idea of generalized hierarchical hidden Markov models (HGHMMs), was used to predict nucleosomal state based on the DNA sequence in yeast chromosome III, and compared with two other existing methods. The HGHMM method performed favorably among the three models in terms of specificity and sensitivity, and provided estimates that were largely consistent with predictions from the method of Yuan and Liu (2008). However, all the methods still give higher than desirable misclassification rates, indicating that sequence-based features may provide only limited information towards understanding positioning of nucleosomes. The method is implemented in the open-source statistical software R, and is freely available from the authors' website.

Identification of homogeneous genetic architecture of multiple genetically correlated traits by block clustering of genome-wide associations.

Genome-wide association studies (GWAS) using high-density genotyping platforms offer an unbiased strategy to identify new candidate genes for osteoporosis. It is imperative to be able to clearly distinguish signal from noise by focusing on the best phenotype in a genetic study. We performed GWAS of multiple phenotypes associated with fractures [bone mineral density (BMD), bone quantitative ultrasound (QUS), bone geometry, and muscle mass] with approximately 433,000 single-nucleotide polymorphisms (SNPs) and created a database of resulting associations. We performed analysis of GWAS data from 23 phenotypes by a novel modification of a block clustering algorithm followed by gene-set enrichment analysis. A data matrix of standardized regression coefficients was partitioned along both axes--SNPs and phenotypes. Each partition represents a distinct cluster of SNPs that have similar effects over a particular set of phenotypes. Application of this method to our data shows several SNP-phenotype connections. We found a strong cluster of association coefficients of high magnitude for 10 traits (BMD at several skeletal sites, ultrasound measures, cross-sectional bone area, and section modulus of femoral neck and shaft). These clustered traits were highly genetically correlated. Gene-set enrichment analyses indicated the augmentation of genes that cluster with the 10 osteoporosis-related traits in pathways such as aldosterone signaling in epithelial cells, role of osteoblasts, osteoclasts, and chondrocytes in rheumatoid arthritis, and Parkinson signaling. In addition to several known candidate genes, we also identified PRKCH and SCNN1B as potential candidate genes for multiple bone traits. In conclusion, our mining of GWAS results revealed the similarity of association results between bone strength phenotypes that may be attributed to pleiotropic effects of genes. This knowledge may prove helpful in identifying novel genes and pathways that underlie several correlated phenotypes, as well as in deciphering genetic and phenotypic modularity underlying osteoporosis risk.

A continuous-index Bayesian hidden Markov model for prediction of nucleosome positioning in genomic DNA.

Nucleosomes are units of chromatin structure, consisting of DNA sequence wrapped around proteins called "histones." Nucleosomes occur at variable intervals throughout genomic DNA and prevent transcription factor (TF) binding by blocking TF access to the DNA. A map of nucleosomal locations would enable researchers to detect TF binding sites with greater efficiency. Our objective is to construct an accurate genomic map of nucleosome-free regions (NFRs) based on data from high-throughput genomic tiling arrays in yeast. These high-volume data typically have a complex structure in the form of dependence on neighboring probes as well as underlying DNA sequence, variable-sized gaps, and missing data. We propose a novel continuous-index model appropriate for non-equispaced tiling array data that simultaneously incorporates DNA sequence features relevant to nucleosome formation. Simulation studies and an application to a yeast nucleosomal assay demonstrate the advantages of using the new modeling framework, as well as its robustness to distributional misspecifications. Our results reinforce the previous biological hypothesis that higher-order nucleotide combinations are important in distinguishing nucleosomal regions from NFRs.

The RNA editor gene ADAR1 is induced in myoblasts by inflammatory ligands and buffers stress response.

Muscle atrophy remains a significant concern in multiple inflammatory conditions, including injury, sepsis, cachexia, and HIV-associated wasting. Herein, we show that inflammatory stressors, including TNF-alpha, IFN-gamma, or lipopolysaccharide, potently induced the novel expression of the RNA editor ADAR1, an observation not previously described in muscle cells. We also observed that cytokine stimulation suppressed muscle-associated microRNAs, an observation also not previously demonstrated. To map potential effects of ADAR1 induction in the muscle program, we conducted knockdown and overexpression studies in the mouse C2C12 muscle precursor cell (MPC) line and in primary human MPCs. We show that knockdown of stress-induced ADAR1 increased inflammation-mediated declines in the muscle differentiation markers Myogenin and myosin heavy chain, and knockdown reduced levels of active phosphorylated Akt (phospho-Akt), but had no effect on microRNA transcript levels, suggesting a role for ADAR1 in buffering inflammatory stress effects on myogenic transcription and protein synthesis pathways. In addition, overexpression of recombinant ADAR1 suppressed active phosphorylated double-stranded RNA (dsRNA)-dependent protein kinase (phospho-PKR), consistent with a role for ADAR1 in limiting inflammation-driven catabolic atrophy pathways. Collectively, these data identify a novel regulatory role for ADAR1 activation under inflammatory stress to both promote muscle protein synthesis pathways and limit atrophy pathways.

A Bayesian hidden Markov model for motif discovery through joint modeling of genomic sequence and ChIP-chip data.

We propose a unified framework for the analysis of chromatin (Ch) immunoprecipitation (IP) microarray (ChIP-chip) data for detecting transcription factor binding sites (TFBSs) or motifs. ChIP-chip assays are used to focus the genome-wide search for TFBSs by isolating a sample of DNA fragments with TFBSs and applying this sample to a microarray with probes corresponding to tiled segments across the genome. Present analytical methods use a two-step approach: (i) analyze array data to estimate IP-enrichment peaks then (ii) analyze the corresponding sequences independently of intensity information. The proposed model integrates peak finding and motif discovery through a unified Bayesian hidden Markov model (HMM) framework that accommodates the inherent uncertainty in both measurements. A Markov chain Monte Carlo algorithm is formulated for parameter estimation, adapting recursive techniques used for HMMs. In simulations and applications to a yeast RAP1 dataset, the proposed method has favorable TFBS discovery performance compared to currently available two-stage procedures in terms of both sensitivity and specificity.

A hierarchical model for incomplete alignments in phylogenetic inference.

MOTIVATION: Full-length DNA and protein sequences that span the entire length of a gene are ideally used for multiple sequence alignments (MSAs) and the subsequent inference of their relationships. Frequently, however, MSAs contain a substantial amount of missing data. For example, expressed sequence tags (ESTs), which are partial sequences of expressed genes, are the predominant source of sequence data for many organisms. The patterns of missing data typical for EST-derived alignments greatly compromise the accuracy of estimated phylogenies. RESULTS: We present a statistical method for inferring phylogenetic trees from EST-based incomplete MSA data. We propose a class of hierarchical models for modeling pairwise distances between the sequences, and develop a fully Bayesian approach for estimation of the model parameters. Once the distance matrix is estimated, the phylogenetic tree may be constructed by applying neighbor-joining (or any other algorithm of choice). We also show that maximizing the marginal likelihood from the Bayesian approach yields similar results to a profile likelihood estimation. The proposed methods are illustrated using simulated protein families, for which the true phylogeny is known, and one real protein family. AVAILABILITY: R code for fitting these models are available from: http://people.bu.edu/gupta/software.htm.

An Information Matrix Prior for Bayesian Analysis in Generalized Linear Models with High Dimensional Data.

An important challenge in analyzing high dimensional data in regression settings is that of facing a situation in which the number of covariates p in the model greatly exceeds the sample size n (sometimes termed the "p > n" problem). In this article, we develop a novel specification for a general class of prior distributions, called Information Matrix (IM) priors, for high-dimensional generalized linear models. The priors are first developed for settings in which p < n, and then extended to the p > n case by defining a ridge parameter in the prior construction, leading to the Information Matrix Ridge (IMR) prior. The IM and IMR priors are based on a broad generalization of Zellner's g-prior for Gaussian linear models. Various theoretical properties of the prior and implied posterior are derived including existence of the prior and posterior moment generating functions, tail behavior, as well as connections to Gaussian priors and Jeffreys' prior. Several simulation studies and an application to a nucleosomal positioning data set demonstrate its advantages over Gaussian, as well as g-priors, in high dimensional settings.

Model selection and sensitivity analysis for sequence pattern models.

In this article we propose a maximal a posteriori (MAP) criterion for model selection in the motif discovery problem and investigate conditions under which the MAP asymptotically gives a correct prediction of model size. We also investigate robustness of the MAP to prior specification and provide guidelines for choosing prior hyper-parameters for motif models based on sensitivity considerations.

Accumulation of M1dG DNA adducts after chronic exposure to PCBs, but not from acute exposure to polychlorinated aromatic hydrocarbons.

Oxidative DNA damage is one of the key events thought to be involved in mutation and cancer. The present study examined the accumulation of M1dG, 3-(2'-deoxy-beta-D-erythro-pentofuranosyl)-pyrimido[1,2-a]-purin-10(3H)-one, DNA adducts after single dose or 1-year exposure to polyhalogenated aromatic hydrocarbons (PHAH) in order to evaluate the potential role of oxidative DNA damage in PHAH toxicity and carcinogenicity. The effect of PHAH exposure on the number of M1dG adducts was explored initially in female mice exposed to a single dose of either 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD) or a PHAH mixture. This study demonstrated that a single exposure to PHAH had no significant effect on the number of M1dG adducts compared to the corn oil control group. The role of M1dG adducts in polychlorinated biphenyl (PCB)-induced toxicity and carcinogenicity was further investigated in rats exposed for a year to PCB 153, PCB 126, or a mixture of the two. PCB 153, at doses up to 3000 microg/kg/day, had no significant effect on the number of M1dG adducts in liver and brain tissues from the exposed rats compared to controls. However, 1000 ng/kg/day of PCB 126 resulted in M1dG adduct accumulation in the liver. More importantly, coadministration of equal proportions of PCB 153 and PCB 126 resulted in dose-dependent increases in M1dG adduct accumulation in the liver from 300 to 1000 ng/kg/day of PCB 126 with 300-1000 microg/kg/day of PCB 153. Interestingly, the coadministration of different amounts of PCB 153 with fixed amounts of PCB 126 demonstrated more M1dG adduct accumulation with higher doses of PCB 153. These results are consistent with the results from cancer bioassays that demonstrated a synergistic effect between PCB 126 and PCB 153 on toxicity and tumor development. In summary, the results from the present study support the hypothesis that oxidative DNA damage plays a key role in toxicity and carcinogenicity following long-term PCB exposure.

Generalized hierarchical markov models for the discovery of length-constrained sequence features from genome tiling arrays.

A generalized hierarchical Markov model for sequences that contain length-restricted features is introduced. This model is motivated by the recent development of high-density tiling array data for determining genomic elements of functional importance. Due to length constraints on certain features of interest, as well as variability in probe behavior, usual hidden Markov-type models are not always applicable. A robust Bayesian framework that can incorporate length constraints, probe variability, and bias is developed. Moreover, a novel recursion-based Monte Carlo algorithm is proposed to estimate the parameters and impute hidden states under length constraints. Application of this methodology to yeast chromosomal arrays demonstrate substantial improvement over currently existing methods in terms of sensitivity as well as biological interpretability.

A temporal hidden Markov regression model for the analysis of gene regulatory networks.

We propose a novel hierarchical hidden Markov regression model for determining gene regulatory networks from genomic sequence and temporally collected gene expression microarray data. The statistical challenge is to simultaneously determine the groupings of genes and subsets of motifs involved in their regulation, when the groupings may vary over time, and a large number of potential regulators are available. We devise a hybrid Monte Carlo methodology to estimate parameters under 2 classes of latent structure, one arising due to the unobservable state identity of genes and the other due to the unknown set of covariates influencing the response within a state. The effectiveness of this method is demonstrated through a simulation study and an application on an yeast cell-cycle data set.

Predictive power of biomarkers of oxidative stress and inflammation in patients with hepatitis C virus-associated hepatocellular carcinoma.

BACKGROUND: This study evaluated the relationship between inflammation, intra-hepatic oxidative stress, oxidative DNA damage and the progression of liver carcinogenesis in hepatitis C virus (HCV)-infected humans. METHODS: Non-cancerous liver tissues were collected from 30 patients with an HCV-associated solitary hepatocellular carcinoma (HCC) who received curative tumor removal. After surgery, the patients were followed at monthly intervals at the outpatient clinic. Distribution of the inflammatory cells (CD68+), the number of 8-hydroxydeoxyguanosine (8-OHdG) DNA adducts and 4-hydroxynonenal (HNE) protein adducts and the expression of apurinic/apyrimidinic endonuclease (APE) were determined by immunohistochemical analysis in serial liver sections from tumor-free parenchyma at the surgical margin around the tumor. RESULTS: Significant positive correlations were observed between the number of CD68+ cells, the amount of HNE protein adducts, and the number of 8-OHdG adducts in liver tissue of patients with HCC and HCV. The cumulative disease-free survival was significantly shorter in patients with the highest percentage of 8-OHdG-positive hepatocytes. Using a Cox proportional hazard model, 8-OHdG, HNE and CD68 were determined to be good biomarkers for predicting disease-free survival in patients with HCC and HCV. CONCLUSIONS: These results support the hypothesis that HCV-induced inflammation causes oxidative DNA damage and promotes hepatocarcinogenesis which directly affects the clinical outcome. Since patients with greater intra-hepatic oxidative stress had a higher incidence of HCC recurrence, we suggest that oxidative stress biomarkers could potentially be used as a useful clinical diagnostic tool to predict the duration of disease-free survival in patients with HCV-associated HCC.

Regulation of nucleosome stability as a mediator of chromatin function.

Interactions among regulatory proteins, enzymes and DNA are required to use the information encoded in genomes. In eukaryotes, the location and timing of interactions between these proteins and DNA is determined in large part by whether DNA at a given locus is packaged into a nucleosome. Given the central role of nucleosome formation in regulating genome function, many mechanisms have evolved to control nucleosome stability at loci across the genome. New conclusions about the role of the DNA sequence itself in the regulation of nucleosome stability have come from two new types of experiment: high-resolution mapping of in vivo nucleosome occupancy on a genomic scale; and in vivo versus in vitro comparisons of nucleosome stability on natural DNA templates. These new types of data raise intriguing questions about the evolutionary constraints on DNA sequence with regard to nucleosome formation, and at long last might enable the derivation of DNA sequence-based rules that produce reliable predictions of nucleosome behavior.