BAYESIAN STATISTICAL METHODS IN GENE-ENVIRONMENT AND GENE-GENE INTERACTION STUDIES

Complex diseases such as cancer result from multiple genetic changes and environmental exposures. Due to the rapid development of genotyping and sequencing technologies, we are now able to more accurately assess causal effects of many genetic and environmental factors. Genome-wide association studies have been able to localize many causal genetic variants predisposing to certain diseases. However, these studies only explain a small portion of variations in the heritability of diseases. More advanced statistical models are urgently needed to identify and characterize some additional genetic and environmental factors and their interactions, which will enable us to better understand the causes of complex diseases. In the past decade, thanks to the increasing computational capabilities and novel statistical developments, Bayesian methods have been widely applied in the genetics/genomics researches and demonstrating superiority over some regular approaches in certain research areas. Gene-environment and gene-gene interaction studies are among the areas where Bayesian methods may fully exert its functionalities and advantages. This dissertation focuses on developing new Bayesian statistical methods for data analysis with complex gene-environment and gene-gene interactions, as well as extending some existing methods for gene-environment interactions to other related areas. It includes three sections: (1) Deriving the Bayesian variable selection framework for the hierarchical gene-environment and gene-gene interactions; (2) Developing the Bayesian Natural and Orthogonal Interaction (NOIA) models for gene-environment interactions; and (3) extending the applications of two Bayesian statistical methods which were developed for gene-environment interaction studies, to other related types of studies such as adaptive borrowing historical data. We propose a Bayesian hierarchical mixture model framework that allows us to investigate the genetic and environmental effects, gene by gene interactions (epistasis) and gene by environment interactions in the same model. It is well known that, in many practical situations, there exists a natural hierarchical structure between the main effects and interactions in the linear model. Here we propose a model that incorporates this hierarchical structure into the Bayesian mixture model, such that the irrelevant interaction effects can be removed more efficiently, resulting in more robust, parsimonious and powerful models. We evaluate both of the 'strong hierarchical' and 'weak hierarchical' models, which specify that both or one of the main effects between interacting factors must be present for the interactions to be included in the model. The extensive simulation results show that the proposed strong and weak hierarchical mixture models control the proportion of false positive discoveries and yield a powerful approach to identify the predisposing main effects and interactions in the studies with complex gene-environment and gene-gene interactions. We also compare these two models with the 'independent' model that does not impose this hierarchical constraint and observe their superior performances in most of the considered situations. The proposed models are implemented in the real data analysis of gene and environment interactions in the cases of lung cancer and cutaneous melanoma case-control studies. The Bayesian statistical models enjoy the properties of being allowed to incorporate useful prior information in the modeling process. Moreover, the Bayesian mixture model outperforms the multivariate logistic model in terms of the performances on the parameter estimation and variable selection in most cases. Our proposed models hold the hierarchical constraints, that further improve the Bayesian mixture model by reducing the proportion of false positive findings among the identified interactions and successfully identifying the reported associations. This is practically appealing for the study of investigating the causal factors from a moderate number of candidate genetic and environmental factors along with a relatively large number of interactions. The natural and orthogonal interaction (NOIA) models of genetic effects have previously been developed to provide an analysis framework, by which the estimates of effects for a quantitative trait are statistically orthogonal regardless of the existence of Hardy-Weinberg Equilibrium (HWE) within loci. Ma et al. (2012) recently developed a NOIA model for the gene-environment interaction studies and have shown the advantages of using the model for detecting the true main effects and interactions, compared with the usual functional model. In this project, we propose a novel Bayesian statistical model that combines the Bayesian hierarchical mixture model with the NOIA statistical model and the usual functional model. The proposed Bayesian NOIA model demonstrates more power at detecting the non-null effects with higher marginal posterior probabilities. Also, we review two Bayesian statistical models (Bayesian empirical shrinkage-type estimator and Bayesian model averaging), which were developed for the gene-environment interaction studies. Inspired by these Bayesian models, we develop two novel statistical methods that are able to handle the related problems such as borrowing data from historical studies. The proposed methods are analogous to the methods for the gene-environment interactions on behalf of the success on balancing the statistical efficiency and bias in a unified model. By extensive simulation studies, we compare the operating characteristics of the proposed models with the existing models including the hierarchical meta-analysis model. The results show that the proposed approaches adaptively borrow the historical data in a data-driven way. These novel models may have a broad range of statistical applications in both of genetic/genomic and clinical studies.

Evolution of cis-regulatory elements in a repetitive sequence adjacent to a sea urchin aboral ectoderm -specific gene

The creation, preservation, and degeneration of cis-regulatory elements controlling developmental gene expression are fundamental genome-level evolutionary processes about which little is known. In this study, critical differences in cis-regulatory elements controlling the expression of the sea urchin aboral ectoderm-specific spec genes were identified and explored. In genomes of species within the Strongylocentrotidae family, multiple copies of a repetitive sequence element termed RSR were present, but RSRs were not detected in genomes of species outside Strongylocentrotidae. RSRs are invariably associated with spec genes, and in Strongylocentrotus purpuratus, the spec2a RSR functioned as a transcriptional enhancer displaying greater activity than RSRs from the spec1 or spec2c paralogs. Single base-pair differences at two cis-regulatory elements within the spec2a RSR greatly increased the binding affinities of four transcription factors: SpCCAAT-binding factor at one element and SpOtx, SpGoosecoid, and SpGATA-E at another. The cis-regulatory elements to which SpCCAAT-binding factor, SpOtx, SpGoosecoid, and SpGATA-E bound were recent evolutionary acquisitions that could act either to activate or repress transcription, depending on the cell type. These elements were found in the spec2a RSR ortholog in Strongylocentrotus pallidus but not in the RSR orthologs of Strongylocentrotus droebachiensis or Hemicentrotus pulcherrimus. These results indicate that spec genes exhibit a dynamic pattern of cis-regulatory element evolution while stabilizing selection preserves their aboral ectoderm expression domain. ^

Detecting reliable gene interactions by a hierarchy of Bayesian network classifiers

The main purpose of a gene interaction network is to map the relationships of the genes that are out of sight when a genomic study is tackled. DNA microarrays allow the measure of gene expression of thousands of genes at the same time. These data constitute the numeric seed for the induction of the gene networks. In this paper, we propose a new approach to build gene networks by means of Bayesian classifiers, variable selection and bootstrap resampling. The interactions induced by the Bayesian classifiers are based both on the expression levels and on the phenotype information of the supervised variable. Feature selection and bootstrap resampling add reliability and robustness to the overall process removing the false positive findings. The consensus among all the induced models produces a hierarchy of dependences and, thus, of variables. Biologists can define the depth level of the model hierarchy so the set of interactions and genes involved can vary from a sparse to a dense set. Experimental results show how these networks perform well on classification tasks. The biological validation matches previous biological findings and opens new hypothesis for future studies

Metagenomic Anlaysis of microsymbiont selection by the legume plant host

Rhizobium leguminosarum bv.viciae is able to establish nitrogen-fixing symbioses with legumes of the genera Pisum, Lens, Lathyrus and Vicia. Classic studies using trap plants (Laguerre et al., Young et al.) provided evidence that different plant hosts are able to select different rhizobial genotypes among those available in a given soil. However, these studies were necessarily limited by the paucity of relevant biodiversity markers. We have now reappraised this problem with the help of genomic tools. A well-characterized agricultural soil (INRA Bretennieres) was used as source of rhizobia. Plants of Pisum sativum, Lens culinaris, Vicia sativa and V. faba were used as traps. Isolates from 100 nodules were pooled, and DNA from each pool was sequenced (BGI-Hong Kong; Illumina Hiseq 2000, 500 bp PE libraries, 100 bp reads, 12 Mreads). Reads were quality filtered (FastQC, Trimmomatic), mapped against reference R. leguminosarum genomes (Bowtie2, Samtools), and visualized (IGV). An important fraction of the filtered reads were not recruited by reference genomes, suggesting that plant isolates contain genes that are not present in the reference genomes. For this study, we focused on three conserved genomic regions: 16S-23S rDNA, atpD and nodDABC, and a Single Nucleotide Polymorphism (SNP) analysis was carried out with meta / multigenomes from each plant. Although the level of polymorphism varied (lowest in the rRNA region), polymorphic sites could be identified that define the specific soil population vs. reference genomes. More importantly, a plant-specific SNP distribution was observed. This could be confirmed with many other regions extracted from the reference genomes (data not shown). Our results confirm at the genomic level previous observations regarding plant selection of specific genotypes. We expect that further, ongoing comparative studies on differential meta / multigenomic sequences will identify specific gene components of the plant-selected genotypes

The FEZ1 gene at chromosome 8p22 encodes a leucine-zipper protein, and its expression is altered in multiple human tumors

Alterations of human chromosome 8p occur frequently in many tumors. We identified a 1.5-Mb common region of allelic loss on 8p22 by allelotype analysis. cDNA selection allowed isolation of several genes, including FEZ1. The predicted Fez1 protein contained a leucine-zipper region with similarity to the DNA-binding domain of the cAMP-responsive activating-transcription factor 5. RNA blot analysis revealed that FEZ1 gene expression was undetectable in more than 60% of epithelial tumors. Mutations were found in primary esophageal cancers and in a prostate cancer cell line. Transcript analysis from several FEZ1-expressing tumors revealed truncated mRNAs, including a frameshift. Alteration and inactivation of the FEZ1 gene may play a role in various human tumors.

Functional copies of a human gene can be directly isolated by transformation-associated recombination cloning with a small 3′ end target sequence

Unique, small sequences (sequence tag sites) have been identified at the 3′ ends of most human genes that serve as landmarks in genome mapping. We investigated whether a single copy gene could be isolated directly from total human DNA by transformation-associated recombination (TAR) cloning in yeast using a short, 3′ unique target. A TAR cloning vector was constructed that, when linearized, contained a small amount (381 bp) of 3′ hypoxanthine phosphoribosyltransferase (HPRT) sequence at one end and an 189-bp Alu repeat at the other end. Transformation with this vector along with human DNA led to selective isolations of the entire HPRT gene as yeast artificial chromosomes (YACs) that extended from the 3′ end sequence to various Alu positions as much as 600 kb upstream. These YACs were retrofitted with a NeoR and a bacterial artificial chromosome (BAC) sequence to transfer the YACs to bacteria and subsequently the BACs to mouse cells by using a Neo selection. Most of the HPRT isolates were functional, demonstrating that TAR cloning retains the functional integrity of the isolated material. Thus, this modified version of TAR cloning, which we refer to as radial TAR cloning, can be used to isolate large segments of the human genome accurately and directly with only a small amount of sequence information.

Cross-lineage expression of Ig-β (B29) in thymocytes: Positive and negative gene regulation to establish T cell identity

Developmental commitment involves activation of lineage-specific genes, stabilization of a lineage-specific gene expression program, and permanent inhibition of inappropriate characteristics. To determine how these processes are coordinated in early T cell development, the expression of T and B lineage-specific genes was assessed in staged subsets of immature thymocytes. T lineage characteristics are acquired sequentially, with germ-line T cell antigen receptor-β transcripts detected very early, followed by CD3ɛ and terminal deoxynucleotidyl transferase, then pTα, and finally RAG1. Only RAG1 expression coincides with commitment. Thus, much T lineage gene expression precedes commitment and does not depend on it. Early in the course of commitment to the T lineage, thymocytes lose the ability to develop into B cells. To understand how this occurs, we also examined expression of well defined B lineage-specific genes. Although λ5 and Ig-α are not expressed, the μ0 and Iμ transcripts from the unrearranged IgH locus are expressed early, in distinct patterns, then repressed just before RAG1 expression. By contrast, RNA encoding the B cell receptor component Ig-β was found to be transcribed in all immature thymocyte subpopulations and throughout most thymocyte differentiation. Ig-β expression is down-regulated only during positive selection of CD4+CD8– cells. Thus several key participants in the B cell developmental program are expressed in non-B lineage-committed cells, and one is maintained even through commitment to an alternative lineage, and repressed only after extensive T lineage differentiation. The results show that transcriptional activation of “lymphocyte-specific” genes can occur in uncommitted precursors, and that T lineage commitment is a composite of distinct positive and negative regulatory events.

Amyrel, a paralogous gene of the amylase gene family in Drosophila melanogaster and the Sophophora subgenus

We describe a gene from Drosophila melanogaster related to the alpha-amylase gene Amy. This gene, which exists as a single copy, was named Amyrel. It is strikingly divergent from Amy because the amino acid divergence is 40%. The coding sequence is interrupted by a short intron at position 655, which is unusual in amylase genes. Amyrel has also been cloned in Drosophila ananassae, Drosophila pseudoobscura, and Drosophila subobscura and is likely to be present throughout the Sophophora subgenus, but, to our knowledge, it has not been detected outside. Unexpectedly, there is a strong conservation of 5′ and 3′ flanking regions between Amyrel genes from different species, which is not the case for Amy and which suggests that selection acts on these regions. In contrast to the Amy genes, Amyrel is transcribed in larvae of D. melanogaster but not in adults. However, the protein has not been detected yet. Amyrel evolves about twice as fast as Amy in the several species studied. We suggest that this gene could result from a duplication of Amy followed by accelerated and selected divergence toward a new adaptation.

The atypical codon usage of the plant psbA gene may be the remnant of an ancestral bias

The psbA gene of the chloroplast genome has a codon usage that is unusual for plant chloroplast genes. In the present study the evolutionary status of this codon usage is tested by reconstructing putative ancestral psbA sequences to determine the pattern of change in codon bias during angiosperm divergence. It is shown that the codon biases of the ancestral genes are much stronger than all extant flowering plant psbA genes. This is related to previous work that demonstrated a significant increase in synonymous substitution in psbA relative to other chloroplast genes. It is suggested, based on the two lines of evidence, that the codon bias of this gene currently is not being maintained by selection. Rather, the atypical codon bias simply may be a remnant of an ancestral codon bias that now is being degraded by the mutation bias of the chloroplast genome, in other words, that the psbA gene is not at equilibrium. A model for the evolution of selective pressure on the codon usage of plant chloroplast genes is discussed.

Somatic mutations in the p53 tumor suppressor gene in rheumatoid arthritis synovium

The factors that regulate the perpetuation and invasiveness of rheumatoid synovitis have been the subject of considerable inquiry, and the possibility that nonimmunologic defects can contribute to the disease has not been rigorously addressed. Using a mismatch detection system, we report that synovial tissue from the joints of severe chronic rheumatoid arthritis patients contain mutant p53 transcripts, which were not found in skin samples from the same patients or in joints of patients with osteoarthritis. Mutant p53 transcripts also were identified in synoviocytes cultured from rheumatoid joints. The predicted amino acid substitutions in p53 were identical or similar to those commonly observed in a variety of tumors and might influence growth and survival of rheumatoid synoviocytes. Thus, mutations in p53 and subsequent selection of the mutant cells may occur in the joints of patients as a consequence of inflammation and contribute to the pathogenesis of the disease.

Evolution of foraging behavior in Drosophila by density-dependent selection

One of the rare examples of a single major gene underlying a naturally occurring behavioral polymorphism is the foraging locus of Drosophila melanogaster. Larvae with the rover allele, forR, have significantly longer foraging path lengths on a yeast paste than do those homozygous for the sitter allele, fors. These variants do not differ in general activity in the absence of food. The evolutionary significance of this polymorphism is not as yet understood. Here we examine the effect of high and low animal rearing densities on the larval foraging path-length phenotype and show that density-dependent natural selection produces changes in this trait. In three unrelated base populations the long path (rover) phenotype was selected for under high-density rearing conditions, whereas the short path (sitter) phenotype was selected for under low-density conditions. Genetic crosses suggested that these changes resulted from alterations in the frequency of the fors allele in the low-density-selected lines. Further experiments showed that density-dependent selection during the larval stage rather than the adult stage of development was sufficient to explain these results. Density-dependent mechanisms may be sufficient to maintain variation in rover and sitter behavior in laboratory populations.

Recurrent invasion and extinction of a selfish gene

Homing endonuclease genes show super-Mendelian inheritance, which allows them to spread in populations even when they are of no benefit to the host organism. To test the idea that regular horizontal transmission is necessary for the long-term persistence of these genes, we surveyed 20 species of yeasts for the ω-homing endonuclease gene and associated group I intron. The status of ω could be categorized into three states (functional, nonfunctional, or absent), and status was not clustered on the host phylogeny. Moreover, the phylogeny of ω differed significantly from that of the host, strong evidence of horizontal transmission. Further analyses indicate that horizontal transmission is more common than transposition, and that it occurs preferentially between closely related species. Parsimony analysis and coalescent theory suggest that there have been 15 horizontal transmission events in the ancestry of our yeast species, through simulations indicate that this value is probably an underestimate. Overall, the data support a cyclical model of invasion, degeneration, and loss, followed by reinvasion, and each of these transitions is estimated to occur about once every 2 million years. The data are thus consistent with the idea that frequent horizontal transmission is necessary for the long-term persistence of homing endonuclease genes, and further, that this requirement limits these genes to organisms with easily accessible germ lines. The data also show that mitochondrial DNA sequences are transferred intact between yeast species; if other genes do not show such high levels of horizontal transmission, it would be due to lack of selection, rather than lack of opportunity.

Somatic hypermutation of the new antigen receptor gene (NAR) in the nurse shark does not generate the repertoire: Possible role in antigen-driven reactions in the absence of germinal centers

The new antigen receptor (NAR) gene in the nurse shark diversifies extensively by somatic hypermutation. It is not known, however, whether NAR somatic hypermutation generates the primary repertoire (like in the sheep) or rather is used in antigen-driven immune responses. To address this issue, the sequences of NAR transmembrane (Tm) and secretory (Sec) forms, presumed to represent the primary and secondary repertoires, respectively, were examined from the peripheral blood lymphocytes of three adult nurse sharks. More than 40% of the Sec clones but fewer than 11% of Tm clones contained five mutations or more. Furthermore, more than 75% of the Tm clones had few or no mutations. Mutations in the Sec clones occurred mostly in the complementarity-determining regions (CDR) with a significant bias toward replacement substitutions in CDR1; in Tm clones there was no significant bias toward replacements and only a low level of targeting to the CDRs. Unlike the Tm clones where the replacement mutational pattern was similar to that seen for synonymous changes, Sec replacements displayed a distinct pattern of mutations. The types of mutations in NAR were similar to those found in mouse Ig genes rather than to the unusual pattern reported for shark and Xenopus Ig. Finally, an oligoclonal family of Sec clones revealed a striking trend toward acquisition of glutamic/aspartic acid, suggesting some degree of selection. These data strongly suggest that hypermutation of NAR does not generate the repertoire, but instead is involved in antigen-driven immune responses.

Identification of semaphorin E as a non-MDR drug resistance gene of human cancers

To improve cancer chemotherapy, a better understanding of the molecular mechanisms of drug resistance is essential. To identify the molecules responsible for drug resistance that is unrelated to MDR1 or MRP gene products, a eukaryotic expression cDNA library of cis-diamminedichloroplatinum(II) (CDDP)-resistant ovarian cancer TYKnuR cells was introduced into Cos-7 cells. After repeated CDDP selection, cDNA homologous to murine semaphorin E was isolated from surviving cells. Human semaphorin E (H-sema E) was overexpressed in CDDP-resistant cell lines and was readily induced not only by diverse chemotherapeutic drugs but also by x-ray and UV irradiation. Transfection of H-sema E conferred a drug-resistant phenotype to CDDP-sensitive cells. In addition, the aberrant expression of H-sema E protein was detected immunohistochemically in 14 of 42 (33.3%) recurrent squamous cell carcinomas removed at autopsy after extensive radiochemotherapy. Recently, another member of the semaphorin family, CD100, was shown to significantly improve the viability of B lymphocytes. These results suggest the involvement of semaphorins in diverse cell survival mechanisms.

Effects of Genome Position and the DNA Damage Checkpoint on the Structure and Frequency of sod2 Gene Amplification in Fission Yeast

The Schizosaccharomyces pombe sod2 gene, located near the telomere on the long arm of chromosome I, encodes a Na+ (or Li+)/H+ antiporter. Amplification of sod2 has previously been shown to confer resistance to LiCl. We analyzed 20 independent LiCl-resistant strains and found that the only observed mechanism of resistance is amplification of sod2. The amplicons are linear, extrachromosomal elements either 225 or 180 kb long, containing both sod2 and telomere sequences. To determine whether proximity to a telomere is necessary for sod2 amplification, a strain was constructed in which the gene was moved to the middle of the same chromosomal arm. Selection of LiCl-resistant strains in this genetic background also yielded amplifications of sod2, but in this case the amplified DNA was exclusively chromosomal. Thus, proximity to a telomere is not a prerequisite for gene amplification in S. pombe but does affect the mechanism. Relative to wild-type cells, mutants with defects in the DNA damage aspect of the rad checkpoint control pathway had an increased frequency of sod2 amplification, whereas mutants defective in the S-phase completion checkpoint did not. Two models for generating the amplified DNA are presented.