GenomicLayers: sequence-based simulation of epi-genomes

Abstract

Background: Cellular development and differentiation in Eukaryotes depends upon sequential gene regulatory decisions that allow a single genome to encode many hundreds of distinct cellular phenotypes. Decisions are stored in the regulatory state of each cell, an important part of which is the epi-genome-the collection of proteins, RNA and their specific associations with the genome. Additionally, further cellular responses are, in part, determined by this regulatory state. To date, models of regulatory state have failed to include the contingency of incoming regulatory signals on the current epi-genetic state and none have done so at the whole-genome level.

Results: Here we introduce GenomicLayers, a new R package to run rules-based simulations of epigenetic state changes genome-wide in Eukaryotes. Simulations model the accumulation of changes to genome-wide layers by user-specified binding factors. As a first exemplar, we show two versions of a simple model of the recruitment and spreading of epigenetic marks near telomeres in the yeast Saccharomyces cerevisiae. By combining the output from 100 runs of the simulation, we generate whole genome predictions of epigenetic state at 1 bp resolution. The example yeast models are included within a 'vignette' with the GenomicLayers package, which is available at https://github.com/davetgerrard/GenomicLayers . To demonstrate the use of GenomicLayers on the full human reference genome (hg38), we show the results from parameter refinement on a simplistic model of the action of pluripotency factors against a self-spreading repressor seeded at CpG islands. The human genome model is included in supplementary information as an R script.

Conclusions: GenomicLayers enables scientists working on diverse eukaryotic organisms to test models of gene regulation in silico. Applications include epigenetic silencing, activation by combinatorial binding of transcription factors and the sink effects caused by down-regulation of components of epigenetic regulators. The software is intended to be used to parameterise, refine and combine models and thereby capitalise on data from the thousands of studies of Eukaryotic epigenomes.

Keywords: Development; Epigenome; Genome; R; Simulation.

Fig. 1

Core components of the GenomicLayers package. a Each simulation requires one ‘layerSet’ comprising a genome sequence as a BSgenome object and one or more layers as GRanges objects. The GRanges ‘layers’ are modified during the simulation and can be arbitrarily named. b One or more ‘bindingFactors’ each contain a ‘profile’ to match with combinations of layer states and/or the genome sequence. Sequence matches can be represented either as consensus sequences (using IUPAC codes), as position weight matrices (shown here as a seqLogo), or as regular expressions. The profile shown here only fully matches (both sequence and layers) at the position in a) shown by the red line (e.g. a sequence match AND LAYER.1 = present AND LAYER.3 = absent; LAYER.2 is ignored). Once a set of matches have been selected, the ‘mods’ of a binding factor dictate changes to the layerSet at that position (e.g. LAYER.2 set to absent AND LAYER.3 set to present). The extent of the modifications is determined by the statewidth parameter, which is independent of the width of the matching pattern. c A simulation comprises one or more cycles during which one or more bindingFactors are matched to the current state of the layerSet. The number of modifications performed by each bindingFactor each cycle is limited by a vector of ‘abundance’ values, which may be kept constant or varied between cycles. The layerSet state after modification becomes the starting state for the next cycle (green arrows) or can be measured against validation data (red arrow)

Fig. 2

Output of whole genome Sir3 recruitment and spreading models M1 and M2 on Saccharomyces cerevisiae chromosome III (SacCer3: chrIII). Output from 100 replicated simulations were stored every 10 cycles (e.g. M2.c50 = Model 2, cycle 50). The outputs were summed using the function coverage() (GenomicRanges package) representing the number of simulations in which each base of the chromosome was marked by Sir3p. Note the high frequency of early Sir3p recruitment at the very tips of the telomeres in both models (denoted by filled arrowheads). Model 1, which allows recruitment of Sir3p from monomeric Rap1 binding, generates Sir3p binding all across the chromosome. Model 2 requires two adjacent Rap1 sites to seed Sir3p binding and has fewer seed sites from which Sir3p binding spreads out along the chromosome. The bottom two plots show Sir3p binding signal for Log-phase and stationary phase samples (light and dark green, respectively) measured using chromatin immunoprecipitation on a genome-wide tiling array (ChIP-chip) [34]. ‘MATALPHA’ marks the mating type-locus, which is repressed by Sir3p binding in the log-phase sample. Plot generated using plotSignal function (plotGardener package [17])

Fig. 3

Proportion of Transcription Start Sites (TSS) covered by ‘active’ (left) or ‘repressed’ (right) states during 100 successive simulation cycles. Sub-plots a-d represent separate simulations using decreasing values of the shared statewidth parameter (shown above each plot). The symbols represent proportions from the three gene lists extracted from Tanaka et al. [33]: Open circles: Type I,Crosses: Type II; Open triangles: Type III. Filled grey circles: All other genes