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.

Methods

Implementation

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

Figure 1.

Workflow of cophesim has three stages: (1) Input stage, where the input data (can be provided in one of the three formats: Plink, MS and GENOME, see the user manual - Supplementary File 1) along with the other input parameters (such as causal variants with size effects, output format, etc.) is prepared for phenotype simulation; (2) Phenotype simulation stage, where different types of phenotypic traits are simulated: dichotomous, continuous and time-to-event (‘survival’); (3) Output stage – the final stage, where simulated phenotype data are packed to various formats in order to be directly usable by six GWAS tools: EMMAX, BLOSSOC, Plink, QTDT, TASSEL and GenABEL. Summary statistics are generated at the output stage as well.

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 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 is the probability of a particular outcome. In cophesim, it is a probability of a “case” (cases are marked by “1”, and “0” are controls in simulated dichotomous phenotype) for i individual. If p is greater than the some threshold p 0 (we use p 0 ~ U(0, 1)), then the phenotype for i individual is set to “1” and to “0” otherwise. The variable z is determined with the following equation: E – effect size for j variant, user-defined; g – value of j genetic marker for i individual; α - effect size for j covariate and X is a value of j covariate for a i individual (the term is added to represent gene-environment iterations); ϵ – a standard normal residual, ϵ (0, 1), computed for i 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 is a weight and computed as follows: (a standardization procedure, and the matrix W containing element w is called a standardized genotype matrix [8]; MAF – a minor allele frequency for j 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 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: S( t| z) = e – , 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 ; 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 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 individual: where the term E 12 g 1 g 2 is the interaction term in which E 12 is the epistatic effect size (user-defined, zero by default); α is the effect size for j covariate X.

Output files

Output files are in the formats as the direct inputs for the following tools: EMMAX [14], Blossoc [4], Plink (.ped file), QTDT [15], TASSEL [16], GenABEL [17] (see Table 1).

Table 1.

Output file formats supported by phenotype simulator cophesim.

Applying one of the options shown below controls the output format. Each output format has a special suffix type, which defines the file format. These output formats are concordant to those used in phenosim.

Application	Option	Commentary
EMMAX	-emmax	Suffices .emma_geno, .emma_pheno
BlOSSOC	-blossoc	Suffices .blossoc_pos, .blossoc_geno
PLINK	-plink	Used by default across all phenotypes, except survival. Suffices .ped, .map, .pheno.
QTDT	-qtdt	Suffices .ped, .map, .dat
TASSEL	-tassel	Suffices .poly, .trait
GenABEL	-	This format is used in simulation of survival phenotype.

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.

Figure 2.

ROC curves constructed from results of association tests performed on a simulated dataset of N = 10,000 individuals, 100 causal and 1,000 of total SNP sites.

TPR: True Positive Ratio, FPR: False Positive Ratio. These results were calculated for dichotomous, continuous and survival traits. The dashed, 45 degrees line represents random guessing.

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

Archived source code as at time of publication: doi: 10.5281/zenodo.810195 [18]

License: MIT

The example script and output files for the software are available at: https://doi.org/10.5281/zenodo.804090 [19].

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.

Zhbannikov and co-workers present the cophesim software that they have developed for simulating phenotypic 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 Consider adding a short paragraph where you discuss limitations and the possibility to add further functionality in the future, including:

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.

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: https://bitbucket.org/izhbannikov/cophesim_data/ROC/roc.R (we were able to retrieve the code from https://bitbucket.org/izhbannikov/cophesim_data/src/)

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.

I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.