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

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

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): p i = 1 1 + e − z i ( 1 ) where p _i 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 ^th individual. If p _i is greater than the some threshold p ₀ (we use p ₀ ~ 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: z i = ∑ j = 1 M E j g i j + ∑ j = 1 K α j X i j + ϵ i ( 2 ) 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 ∑ j = 1 K α j X i j 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: z i = ∑ j = 1 M w i j + ∑ j = 1 K α j X i j + ϵ i ( 3 )

Here w _ij is a weight and computed as follows: w i j = g i j − 2 M A F j ( 2 M A F j ( 1 − M A F j ) ) 1 / 2 (a standardization procedure, and the matrix W containing element w _ij is called a standardized genotype matrix ⁸ ; 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.

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).

We model a survival phenotype from the proportional hazards model using the inverse probability method ¹³ : 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 ^{– H ₀( t) e z} , then the random variable T = S − 1 ( ⋅ | z ) = H 0 ( t ) − 1 ( − log ( U ) e z ) ( 4 ) has survival function S(·| z). In this equation, H ₀( t) is a cumulative baseline hazard. By default, we use the Weibull cumulative baseline hazard: H ₀( 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.

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

These are modeled with the following equation for i ^th individual: z i = E 1 g 1 i + E 2 g 2 i + E 12 g 1 i g 2 i + ∑ j = 1 k α i X j i + ϵ i ( 5 ) where the term E ₁₂ g _{1 i} g _{2 i} is the interaction term in which E ₁₂ is the epistatic effect size (user-defined, zero by default); α _j is the effect size for j ^th covariate X.

Output files are in the formats as the direct inputs for the following tools: EMMAX ¹⁴ , Blossoc ⁴ , Plink (.ped file), QTDT ¹⁵ , TASSEL ¹⁶ , GenABEL ¹⁷ (see Table 1).

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.

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