- Proceedings
- Open Access
- Published:

# Simulation of a medication and methylation effects on triglycerides in the Genetic Analysis Workshop 20

*BMC Proceedings***volume 12**, Article number: 25 (2018)

## Abstract

The GAW20 simulation data set is based upon the companion Genetics of Lipid Lowering Drugs and Diet Network (GOLDN) study fenofibrate clinical trial data set that forms the real data example for GAW20. The simulated data problem consists of 200 simulated replications of what might happen if we were to repeat the GOLDN clinical trial 200 independent times, for these exact same subjects, but using a new fictitious drug (called “genomethate”) that has a pharmaco-epigenetic effect on triglyceride response. For each replication, the pre-genomethate values at visits 1 and 2 are constant (ie, pedigree structures, age, sex, all phenotypes, covariates, genome-wide association study (GWAS) genotypes, and visit 2 methylation values), the same as the real GOLDN data across all 200 replications. Only the post-genomethate treatment data (ie, methylation and triglyceride levels for visits 3 and 4) change across the 200 replications. We postulate a growth curve pharmaco-epigenetic response model, in which each patient’s response to genomethate treatment is individualized, and is dependent upon their genotype as well as the methylation state for key genes.

## Background

The companion Genetics of Lipid Lowering Drugs and Diet Network (GOLDN) study fenofibrate clinical trial data set [1,2,3] was the foundation of our Genetic Analysis Workshop 20 (GAW20) simulation. The general simulation strategy was to first simulate visit 4 methylation array data for each subject (which measures the individual epigenetic responses to genomethate treatment), and then use this plus the genome-wide association study (GWAS) genotypes to produce the simulated triglycerides for visits 3 and 4 post-treatment values. The main simulated effect of genomethate is on the phenotype of the individual subject’s triglyceride (TG) values measured as slope in response to treatment (change in mg/dL per unit time of treatment).

## Methods, results and discussion

Figure 1 illustrates the graphical design of the simulations.

The *j* index in the figure represents the subject (*j* = 1, 2, …, *N* = 717). The *i* index is noting the single-nucleotide polymorphisms (SNPs) chosen to be causal in the simulating model (*i* = 1, 2, …, G = 105), where *i* = 1, 2, 3, 4, 5 also indexes the 5 main effects of the corresponding nearby cytosine-phosphate-guanine (CpG) sites, while beyond main effects, the sites from 6 to 105 are 100 SNPs with background genetic effects. The *k* index indicates replications (*k* = 1, 2, …, *R* = 200).

The first 5 causal SNPs are “major” effects (summarized in Table 1), and the last 100 SNPs are polygenic background effects (Table 2). Note that only the first 5 CpG sites are relevant to the model, the polygenic background effects do *not* depend upon CpG states.

We first defined a series of subjects’ triglyceride values from the original (real) Genetics of Lipid Lowering Drugs and Diet Network (GOLDN) data [1], which was used to generate the simulations. Because triglycerides were approximately log-normally distributed, we worked with log-transformed triglyceride values in all calculations, only transforming back to the measured triglyceride scale at the end of the simulations. In particular, for the *j*th subject, the average log triglycerides pre-treatment (average of visits 1 and 2, which are 1 day apart) and post-treatment (average of visits 3 and 4, which are also 1 day apart) in the original (real) GOLDN data are:

where ** O** –stands for “Observed / Original”,

**stands for “pre-medication treatment,”**

*preRx***stands for “after medication treatment,” and**

*postRx***labels “triglycerides” which were log transformed to ensure a normal distribution of the trait. The**

*TG***of person**

*TG**j*is measured in visits 1, 2, 3 and 4 and averaged as above for each individual as

**and**

*preRx***. The corresponding change in log triglycerides pre-treatment to post-treatment for subject**

*postRx**j*is given by:

where ** delta** is the “change”. The individual time on treatment (less than 30 days) for each subject (in days), is given by the following formula:

where ** daysRx** is “days after medication treatment,”

**_**

*draw***is “blood draw date” at a particular**

*date*

*v***- “**visit.” Thus, the observed slope (change in log triglycerides over the treatment period) is:

If ** mean**_

**_**

*O***_**

*PreRx***and**

*TG***_**

*sd***_**

*O***_**

*preRx***are the mean and standard deviations, respectively, of all the**

*TG***_**

*O***_**

*preRx*

*TG*_{j}across the

*j*= 1, …, N individuals, then the standardized original

**of**

*preRx*

*TG*_{j}are given by:

where ** O**_

**-is a standardized normally distributed variable with**

*preZ**N(0,1)*.

Tables 1 and 2 summarize the epigenetic model in our simulation. We chose 5 “major gene” causal variants (ranging from modest to small effect sizes corresponding to expected “heritabilities” of 0.125, 0.10, 0.075, 0.05, and 0.025), which, in the absence of any epigenetic effects, should govern individual genomethate treatment response along with 100 polygene variants (each of tiny effect size corresponding to “heritabilities” of 0.001 each). These were chosen randomly from chromosomes 1–20 of the GWAS Affymetrix Genome-wide Human SNP Array 6.0, which had 718,544 autosomal SNPs.

For the epigenetic component, we choose 5 CpG sites on the Illumina Infinium HumanMethylation450 BeadChip array (which had 463,995 CpG sites) that are physically closest to the 5 “major gene” causal SNPs, while the methylation sites near the 100 polygenes have no effect. The genomethate response model is based upon the idea that these CpG sites need to be sufficiently unmethylated for the corresponding causal SNPs to express their influence on each individual’s phenotype. If the nearby CpG site is totally methylated (=1), then the corresponding causal SNP actually has *no* effect on the phenotype. If the CpG site is totally unmethylated (=0), then the corresponding causal SNP carries its full effect size impact on the phenotype. If the CpG site is partially methylated (between 0 and 1), then the effect size of the causal SNP is proportionally attenuated.

Specifically, for the *k*th simulation, we first generated the simulated visit 4 methylation array results for all subjects, based upon their corresponding visit 2 and/ or visit 4 methylation values. For each subject *j* = 1, …, 717, and each CpG methylation site *i* = 1, 2, 3, 4, 5 (corresponding to 5 major effect CpGs)

where ** sim**_

**stands for “simulated methylation” at visit 4,**

*meth***_**

*real***is the**

*meth**j*th subject’s “real methylation” array data at visit 2 for the

*i*th CpG site,

*sd*_{i}= 0.4 represents the standard deviation of individual subject methylation responses to treatment, and

*Z***1**

*jik**~ N(0, 1)*is a pseudo-random standard normal variable drawn independently for each

*jik*.

For the remaining, non-causal CpG sites, if the subject *j* had real visit 4 methylation array data then

Otherwise, if the subject *j* only had visit 2 methylation array data, then

where ** real**_

**_**

*meth*

*v***2**

_{ji}and

**_**

*real***_**

*meth*

*v***4**

_{ji}are the real visit 2 and visit 4 methylation array data, respectively, for subject

*j*and CpG site

*i*,

*sd*_{i}represents the standard deviation of individual subject methylation responses to treatment for the

*i*th CpG site, and again,

*Z***1**

*jik**~ N(0,1)*is a pseudo-random variable drawn independently for each

*jik*.

We selected five random non-causal (red-herrings) CpG sites also (shown in Table 3). We set for them the *sd*_{i} = 0.4, to be similar to the simulated causal CpG sites. For the remaining non-causal CpG sites, we set the corresponding *sd*_{i} = 0.03, which is closer to that seen in the real visit 4 methylation data CpG sites, essentially at the measurement error level.

In all cases, *all* simulated visit 4 methylation values were then truncated to be strictly in the [0,1] interval, that is,

if (** sim**_

**_**

*meth*

*v***4**

_{jik}

*>***1**) then

**_**

*sim***_**

*meth*

*v***4**

_{jik}

*=***1**

if (** sim**_

**_**

*meth*

*v***4**

_{jik}

*<***0**) then

**_**

*sim***_**

*meth*

*v***4**

_{jik}

*=***0**

for all subjects *j*, CpG sites *i*, and simulation replications *k*.

Note that the model is such that, on average, the genomethate treatment has no effect on the amount of methylation increase/decrease from visit 2 to visit 4, however, there is variability across subjects. To reiterate, the variability is quite high (*sd*_{i} = 0.4) for the five CpG regions controlling the expression of the major causal variants and 5 other non-causal CpG (red-herrings) sites. The variability is low (*sd*_{i} = 0.03) for all other CpGs, at the level of measurement error.

Using these simulated visit 4 methylation data, we then generated the simulated slope change in triglyceride response for each individual *j* in each replication *k* as follows:

In the above formula, *zenv*_{jk} is an independently drawn pseudo-random normal deviate distributed *N(0,1)* for each subject *j* and each replication *k,* and it represents unexplained residual variation in the phenotype. *SSNP*_{ji} is the standardized *i*th SNP additive genotype-dosage (i.e., coded such that mean = 0 and *sd*_{i} = 1 in the sample), and the *i* = 1, 2, …, 105 regression coefficients in this linear model are given in terms of constants ** sqrt**(

*hg***2**

_{i}

*)*, in Tables 1 and 3. Note that if the five causal CpG sites were completely unmethylated for all subjects (i.e., no epigenetic effects), then (

*1 –*

**_**

*sim***_**

*meth*

*v***4**

_{jik}) would be = 1 for all

*j*= 1,…, N and

*i*= 1,…, 5, and

*k*= 1,…, 200, so that the regression coefficients would be interpreted as the square root of the locus specific heritability of the associated SNPs. Conversely, when the causal CpG site is totally methylated for that subject, (1

*-*

**_**

*sim***_**

*meth*

*v***4**

_{jik}) = 0, so that the corresponding major effect SNP

_{i}will not express its effect on the phenotype. Similarly, if the CpG site is partially methylated (between 0 and 1), the effect size of the causal SNP is proportionally attenuated.

To carry forward these simulated relationships in eq. (1), we must address the fact that the observed slope responses for each subject are correlated to their baseline values of triglyceride (i.e., lower baseline values should produce less dramatic declines with treatment, whereas higher baseline values can experience greater slope change with treatment). In the real GOLDN data, the correlation between slope change in response to fenofibrate treatment and baseline log triglycerides is − 0.41881, and we used this constant value in our genomethate simulation to introduce a correlation between slope change and baseline values:

Because the simulated individual slopes are generated on the standardized scale, we needed to rescale to that of the original scale of triglyceride changes per day of treatment, by working backwards. The mean and standard deviation of ** O**_

**_**

*slope*

*TG*_{j}over all subjects

*j*, are denoted by

**_**

*mean***_**

*O***_**

*slope***and**

*TG***_**

*sd***_**

*O***_**

*slope***, respectively. We used the above observed mean and standard deviation of slopes seen in the original GOLDN data, to rescale as follows:**

*TG*Then the expected response to genomethate treatment of the *j*th subject, after ** O**_

*DaysRx*_{j}original days of treatment, is given by:

Finally, we used the simulated individual responses to produce the simulated values of triglyceride at visits 3 and 4, based upon the variability we see between those visits in the real GOLDN fenofibrate data:

If only 1 replicate of the GAW20 simulated data was to be analyzed, we recommend the 84th replication, which was provided in a separate directory, as a “representative” of the 200 replicated simulations. Chromosomes 21 and 22 datasets were not used in the simulation, so an analyst can use the corresponding data for building a *NULL* hypothesis. The simulated GAW20 data are accompanied by README and Data Dictionary files.

## References

- 1.
Tintle NL, Fardo DW, deAndrade M, Aslibekyan S, Bailey JN, Bermejo JL, Cantor RM, Ghosh S, Melton P, Wang X, MacCluer JW, Almasy L. GAW20: methods and strategies for the new frontiers of epigenetics and pharmacogenomics. BMC Proc. 2018;12(Suppl 9) https://doi.org/10.1186/s12919-018-0113-1.

- 2.
Irvin MR, Zhi D, Joehanes R, Mendelson M, Aslibekyan S, Claas SA, Thibeault KS, Patel N, Day K, Jones LW, et al. Epigenome-wide association study of fasting blood lipids in the genetics of lipid-lowering drugs and diet network study. Circulation. 2014;130(7):565–72.

- 3.
Aslibekyan S, Almasy L, Province MA, Absher DM, Arnett DK. Data for GAW20: genome-wide DNA sequence variation and epigenome-wide DNA methylation before and after fenofibrate treatment in a family study of metabolic phenotypes. BMC Proc. 2018;12(Suppl 9) https://doi.org/10.1186/s12919-018-0114-0.

## Acknowledgements

This study was supported in part by the NHLBI grant R01HL117078

### Funding

Publication of this article was supported by NIH R01 GM031575.

### Availability of data and materials

The data that support the findings of this study are available from the Genetic Analysis Workshop (GAW), but restrictions apply to the availability of these data, which were used under license for the current study. Qualified researchers’ may request these data directly from GAW.

### About this supplement

This article has been published as part of BMC Proceedings Volume 12 Supplement 9, 2018: Genetic Analysis Workshop 20: envisioning the future of statistical genetics by exploring methods for epigenetic and pharmacogenomic data. The full contents of the supplement are available online at https://bmcproc.biomedcentral.com/articles/supplements/volume-12-supplement-9.

## Author information

### Affiliations

### Contributions

MAP commenced the idea of GAW20 simulation and wrote the equations. ATK selected main and background SNPs and CpGs for the causative model. ATK, PL and JEH programmed the GAW20 simulations. EWD, SJL and CW performed QC analyses of all replications for GAW20 simulation. PA facilitated the retrieval of GOLDN real data, on which the GAW20 simulation was built. ATK and MAP wrote the manuscript with contributions of EWD, PA, PL, JEH, SJL, and CW. All authors participated in all meetings of GAW20 simulation working group and have read and approved this manuscript.

### Corresponding authors

Correspondence to Aldi T. Kraja or Michael A. Province.

## Ethics declarations

### Ethics approval and consent to participate

Not applicable.

### Consent for publication

Not applicable.

### Competing interests

The authors declare that they have no competing interests.

### Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Rights and permissions

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

## About this article

#### Published

#### DOI