cophesim: a comprehensive phenotype simulator for testing novel association methods

Simulation is important in evaluating novel methods when input data is not easily obtainable or specific assumptions are needed. We present cophesim, a software to add the phenotype to generated genotype data prepared with a genetic simulator. The output of cophesim can be used as a direct input for different genome wide association study tools. cophesim is available from https://bitbucket.org/izhbannikov/cophesim.


Introduction
Genome-wide association studies (GWAS) are routine in population research. New methods are being developed for better accessing complex associations between genotypes and phenotypes, uncovering genotype structures or testing evolutionary hypotheses. Testing the novel methods requires experimental data, which may not be easily obtainable. One solution is to use artificial data simulated with specific assumptions.
The best existing phenotype simulators, such as: GENOME 1 , Plink 2 , phenosim 3 , CoaSim 4 , Fregene 5 , ForSim 6 , QuantiNemo 7 , GCTA 8 , HapGen 9 , SeqSimla 10 , and SimRare 11 offer qualitative and dichotomous simulated phenotype. But the known phenotype simulation software tools have some limitations, which may prevent customers from using them: (i) the majority, if not all, of the phenotype simulation software tools do not offer simulation of survival traits/ time-to-event outcome, making it impossible to test respective hypotheses of associations; (ii) some of the tools are not easy to use, due to wide range of parameters, which the user has to provide and control (rather than calculate them automatically), making them unnecessarily difficult to use and preventing the user from future use of the tool; (iii) phenotype simulation is often offered as an auxiliary part of the genetic simulation routine, and therefore the user first has to perform a time-consuming unavoidable genetic simulation in order to obtain the phenotype; (iv) in situations when the genetic data is already simulated from other tools, only phenosim and GCTA offer adding simulated phenotype to such data. Consequently, it is necessary to have a new, simple and flexible phenotype simulation tool with plain algorithmic assumptions.
Consequently, we present cophesim, a comprehensive phenotype simulation tool that was developed to add a phenotype to corresponding genotypes simulated by other simulation tool (Table S1). cophesim offers simulation of continuous, dichotomous and survival traits, with different (user-provided) effect sizes of causal variants, with the ability to simulate epistatic interactions. It also can simulate phenotype within gene-environment interaction assumptions using up to 10 covariates.

Implementation
The workflow (see Figure 1) includes the following stages: (i) Input data pre-processing; (ii) phenotype simulation; (iii) generation of final output files.

Input data
Currently cophesim accepts the genotype output data from Plink, MS 12 and GENOME software applications. Phenotypes (dichotomous, continuous and survival) are then added according to the following simulation scenarios.

Dichotomous phenotype
Dichotomous phenotype for i th individual (i = 1...N, where N is the total number of individuals in a dataset) is simulated according to the logistic model (if the user provided effect sizes for causal variants): where p i is the probability of a particular outcome. In cophesim, it is a probability of a "case" (cases are m~arked by "1", and "0" are controls in simulated dichotomous phenotype) for i th individual. If p i is greater than the some threshold p 0 (we use p 0 ~ U(0, 1)), then the phenotype for i th individual is set to "1" and to "0" otherwise. The variable z is determined with the following equation: E j -effect size for j th variant, user-defined; g ij -value of j th genetic marker for i th individual; α j -effect size for j th covariate and X ij is a value of j th covariate for a i th individual (the term 1 K j ij j X = α ∑ is added to represent gene-environment iterations); ε i -a standard normal residual, ε i N (0, 1), computed for i th individual, M is a total number of genetic variants and K is a total number of covariates used.
If the user did not provide the effect sizes for causal variants, the following strategy is then used: Here w ij is a weight and computed as follows: (a standardization procedure, and the matrix W containing element w ij is called a standardized genotype matrix 8 ; MAF j -a minor allele frequency for j th genetic variant, and the other values are the same as described above. This strategy allows using defined genetic architecture in a simulated population.

Continuous phenotype
Qualitative (continuous) phenotype for i th individual is simulated according to the linear regression scenario according to the equations (2) or (3) (in case if effect sizes were not supplied).

Inverse Probability method
We model a survival phenotype from the proportional hazards model using the inverse probability method 13 : if U is uniform in (0, 1) and if S(·|z) is the conditional survival function derived from the proportional hazards model: , then the random variable has survival function S(·|z). In this equation, H 0 (t) is a cumulative baseline hazard. By default, we use the Weibull cumulative baseline hazard: H 0 (t) = λρt ρ-1 ; z is the same parameter that defined above, for each individual, and depends on whether the user provided effect sizes for causal variants or not. We also implemented exponential and Gompertz hazards.

Linkage Disequilibrium
The simplest way to simulate collinearity between two SNPs, g 1 and g 2 , with effect sizes E 1 and E 2 is to replace some portion of g 2 with g 1 values according to provided 2 12 r coefficient, which reflects a correlation between two SNPs. We also consider applying other techniques, such as copulas, in order to simulate LD.

Epistatic interactions
These are modeled with the following equation for i th individual: where the term E 12 g 1i g 2i is the interaction term in which E 12 is the epistatic effect size (user-defined, zero by default); α j is the effect size for j th covariate X.
Operation cophesim is freely available for download from the following link: https://bitbucket.org/izhbannikov/cophesim. Requirements: Python v2.7.10 and newer, plinkio v0.9.6, R v3.2.4 and newer, Plink v1.07, -in order to run the examples. The user manual is provided in a separate file "cophesim.pdf" located in the program directory and is also available as Supplementary File 1.

Use case
Below we present an example that shows simulation of genetic data and then simulation of three different phenotypic traits. Other examples and installation instructions are provided at the program website and also in the user manual. Refer to the user manual for description of input parameters. In this example, we first (Step 1) simulate genetic data using Plink. We simulate N.cases = N.control = 5,000 cases and controls and 1,000 SNPs (defined in wgas.sim file, refer to the Plink website to see documentation for this type of file). Then (Step 2) we convert a binary sim.plink.bed file to sim.plink.ped (option --recode in Plink). This step is not required since cophesim can handle binary Plink files (.bed, .bim, .fam), but we provide this step in order to show the ability of the program to deal with Plink PED format. Finally (Step 3), we simulate dichotomous (by default), continuous (option -c) and survival (option -s) traits from previously simulated data stored in files sim.plink.ped and sim.plink.map. Note that we simulate survival trait with Gompertz hazard function (option -gomp); effect sizes for causal variants are provided in file effects.txt (to include this file we use option -ce).

ROC curves
We provide Receiver-Operating Characteristic (ROC) curves ( Figure 2) constructed from association tests performed on a simulated dataset. Simulation and association testing were performed with Plink suite. The following parameters were used: N = 10,000 individuals, N.snp.c = 100 causal, with total N.snp = 1,000 variants. Causal variants were labeled with '1' and the other (neutral) variants were labeled with '0'. These labels are then used later as true identifiers during calculation of TPR (true positive rate) and FPR (false positive rate). Dichotomous, continuous and survival phenotypic traits were simulated with cophesim. Then association tests were performed with Plink for dichotomous and continuous traits (using Plink flags -logistic and -regression, respectively). Association tests for survival trait were performed with the R package GenABEL. Then calculated p-values provided by association tests for each variant were compared to the significance threshold. Those variants passed the threshold were recognized as causal and associated with simulated phenotype. These classification results later were compared to the true identifiers (defined above) in order to obtain TPR and FPR. For all these tests, we varied the significance threshold from 0 to 1 with the increment of 0.001.
The R code to construct ROC curves is provided in the file "roc.R". This file is attached to this computer note and also in the data repository: https://bitbucket.org/izhbannikov/cophesim_data/ROC/roc.R

Conclusion
In this work we presented the cophesim for phenotype simulation from genetic data obtained either from simulation or real data collecting. cophesim makes it possible to simulate various demographic models under user-defined scenarios.

Software and data availability
Tool and source code available from: https://bitbucket.org/ izhbannikov/cophesim To test the cophesim we provided a repository "cophesim_data": https://bitbucket.org/izhbannikov/cophesim_data. Download or clone this repository to be able to run tests.

Competing interests
No competing interests were disclosed.

Grant information
This work was supported by the National Institute on Aging of the National Institutes of Health (NIA/NIH) under award numbers P01AG043352, R01AG046860, and P30AG034424.

The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIA/NIH.
The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Zhbannikov and co-workers present the software that they have developed for simulating phenotypic cophesim data using genotype information. The input and output file formats are compatible with many of the most commonly used computer programs for genome wide association studies. The software is flexible, well documented and fills a gap in existing tools, especially for simulating time-to-event phenotypes. The paper is well written and easy to follow, and we only have some minor comments and suggestions.

Minor comments
Closing bracket missing in the sentence below equation (3) In equation (4), if the user does not provide gene effects then the phenotype is built by the sum of the standardized genotypes for each individual. Could you motivate this choice a bit and explain why it would be useful? In equation (5), the subscripts look wrong. a_i should be a_j In the Linkage Disequilibrium section the term "copula" is used. We do not think most readers of this paper can be expected to be acquainted with copulas and a reference is needed. Consider adding a short paragraph where you discuss limitations and the possibility to add further functionality in the future, including: Dominance effects Probit link for binary data Simulation of correlated traits Alternative ways to simulate LD including a copula approach Check that the following link (at the end of the paper) works: (we were able to retrieve the code from https://bitbucket.org/izhbannikov/cophesim_data/ROC/roc.R ) https://bitbucket.org/izhbannikov/cophesim_data/src/

Is sufficient information provided to allow interpretation of the expected output datasets and any results generated using the tool? Yes
Are the conclusions about the tool and its performance adequately supported by the findings presented in the article? Yes No competing interests were disclosed.

Competing Interests:
We have read this submission. We believe that we have an appropriate level of expertise to confirm that it is of an acceptable scientific standard. In the manuscript "Cophesim: a comprehensive phenotype simulator for testing novel association methods", I. Zhbannikov and colleagues from the Duke University, presented a software that allowed to generate genotype data prepared with a genetic simulator for the use in the investigations of the genome wide association study (GWAS) tools. The rational for development of the software is clearly explained. The idea of the study is to use computer simulations to model data with specific assumptions. Similar simulators are known but all of them do not allow simulate survival. There are several other disadvantage with the existing simulators reviewed by the authors.
The description of the software is technically sound. The methods section is clearly presented. Dichotomous phenotype are simulated according to the logistic model with the covariates being genetic variants and covariates. Continuous phenotypes are simulated using the linear regression. Survival phenotype is modeled using the proportional hazards with inverse probability method.
The details of the code, methods and analysis allow replication of the software and its use by the others. The methods section is clearly presented. Dichotomous phenotype are simulated according to the logistic model with the covariates being genetic variants and covariates. The output formats are compatible with the other applications (Table 1). It is useful example if using the simulator and the other examples are available in the manual. The ROC curve example is also very useful. The information provided is quite sufficient to allow interpretation of the expected results.
In short, the Cophesim is a useful tool that can be helpful in the genetic analyses. The article is scientifically sound, the methods are described with details -this article will greatly help the researcher interested in the application genetic analyses.
Is the rationale for developing the new software tool clearly explained?