| Title: | Functions for Simulating Genotype-Phenotype Relationships |
|---|---|
| Description: | Functions for simulating SNP effects on phenotypes. Functions for simulating disease risk on probability and liability scales. Simulations to quickly obtain estimates of effect sizes for simulated genotype-phenotypes. |
| Authors: | Gibran Hemani [aut, cre] (ORCID: <https://orcid.org/0000-0003-0920-1055>), John Ferguson [aut], Rita Rasteiro [aut] (ORCID: <https://orcid.org/0000-0002-4217-3060>) |
| Maintainer: | Gibran Hemani <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.1.3 |
| Built: | 2026-10-08 15:55:47 UTC |
| Source: | https://github.com/explodecomputer/simulateGP |
This package is a collection of utilities relating to simulating genotype phenotype maps The simulations are largely used for Mendelian randomisation, and link up with the TwoSampleMR R package
Full documentation available here: https://explodecomputer.github.io/simulateGP
Maintainer: Gibran Hemani [email protected] (ORCID)
Authors:
Gibran Hemani [email protected] (ORCID)
John Ferguson [email protected]
Rita Rasteiro [email protected] (ORCID)
Useful links:
Report bugs at https://github.com/explodecomputer/simulateGP/issues
Add confounder variables and their instruments
add_u(parameters, nsnp_u, var_u.x, var_u.y, var_gu.u)add_u(parameters, nsnp_u, var_u.x, var_u.y, var_gu.u)
parameters |
Output from |
nsnp_u |
nsnp_u |
var_u.x |
var_u.x |
var_u.y |
var_u.y |
var_gu.u |
var_gu.u |
List of model parameters
Estimate allele frequency from SNP
allele_frequency(g)allele_frequency(g)
g |
Vector of 0/1/2 |
Allele frequency
Create an arbitrary map
arbitrary_map(af)arbitrary_map(af)
af |
vectory of allele frequencies |
data frame with snp name, chr, pos, ea, oa
Ascertain some proportion of cases and controls from binary phenotype
ascertain_samples(d, prop_cases)ascertain_samples(d, prop_cases)
d |
Vector of 1/0 |
prop_cases |
Proportion of 1s to retain |
Array of IDs
Get vector of effects that explain some amount of variance
choose_effects(nsnp, totvar, sqrt = TRUE, mua = 0)choose_effects(nsnp, totvar, sqrt = TRUE, mua = 0)
nsnp |
Number of SNPs |
totvar |
Total variance explained by all SNPs |
sqrt |
Output effect sizes in terms of standard deviations. Default=TRUE |
mua |
Constant term to be added to effects. Default = 0 |
Vector of effects
Columns are the case and control frequencies Rows are the frequencies for allele 1 and allele 2
contingency(af, prop, odds_ratio, eps = 1e-15)contingency(af, prop, odds_ratio, eps = 1e-15)
af |
Allele frequency of effect allele |
prop |
Proportion of cases |
odds_ratio |
Odds ratio |
eps |
tolerance. Default = 1e-15 |
2x2 contingency table as matrix
Based on the parameters specified this function will call init_parameters, add_u, sample_system_effects, simulate_population and estimate_system_effects. A separate population is generated for each phenotype (x, y and each of u) to allow 2SMR and PRS analyses
create_system( nidx, nidy, nidu = 0, nu = 0, na = 0, nb = 0, var_x.y, nsnp_x, var_gx.x, var_gx.y = 0, mu_gx.y = 0, prop_gx.y = 1, nsnp_y = 0, var_gy.y = 0, var_gy.x = 0, mu_gy.x = 0, prop_gy.x = 1 )create_system( nidx, nidy, nidu = 0, nu = 0, na = 0, nb = 0, var_x.y, nsnp_x, var_gx.x, var_gx.y = 0, mu_gx.y = 0, prop_gx.y = 1, nsnp_y = 0, var_gy.y = 0, var_gy.x = 0, mu_gy.x = 0, prop_gy.x = 1 )
nidx |
number of individuals for x |
nidy |
number of individuals for y |
nidu |
Default=0 If 0 then don't simulate separate populations for the u variables |
nu |
Default=0 Number of confounders influencing x and y |
na |
Default=0 Number of traits upstream of x |
nb |
Default=0 Number of traits upstream of y |
var_x.y |
var in y explained by x |
nsnp_x |
nsnp influencing x |
var_gx.x |
variance in x explained by instruments for x |
var_gx.y |
Default=0 |
mu_gx.y |
Default=0 |
prop_gx.y |
Default=1 |
nsnp_y |
Default=0 |
var_gy.y |
Default=0 |
var_gy.x |
Default=0 |
mu_gy.x |
Default=0 |
prop_gy.x |
Default=1 |
List of effects across system
Estimate the effects of all SNPs on all phenotypes
estimate_system_effects(sim)estimate_system_effects(sim)
sim |
Output from |
Lists of SNP-trait effect estimates
Calculate expected MSE
expected_mse(beta, af, vy)expected_mse(beta, af, vy)
beta |
array of effect sizes |
af |
array of allele frequencies |
vy |
variance of y |
Numeric
se = sqrt(sigma_e^2 / ss(x))
expected_se(beta, af, n, vy)expected_se(beta, af, n, vy)
beta |
array of effect sizes |
af |
array of allele frequencies |
n |
sample size |
vy |
variance of y |
array of standard errors
Calculate expected SSX
expected_ssx(af, n)expected_ssx(af, n)
af |
array of allele frequencies |
n |
sample size |
Numeric
Get summary statistics in simple linear regression
fast_assoc(y, x)fast_assoc(y, x)
y |
Vector of dependent variable |
x |
Vector of independent variable |
List
Generate SNP effects given MAF, h2 and selection
generate_gwas_params(map, h2, S = 0, Pi = 1)generate_gwas_params(map, h2, S = 0, Pi = 1)
map |
Data frame containing at least |
h2 |
Variance explained by all SNPs |
S |
Selection coefficient on trait. Default = 0 |
Pi |
Proportion of variants that have an effect - sampled randomly. Default=1 |
data frame
After generating the set of causal effects at each SNP, use an LD correlation matrix to transform the effects to reflect the correlation structure at the SNPs. Note if running many repeats, only need to generate the LD-modified params once and then can repeatedly re-sample using generate_gwas_ss
generate_gwas_ss( params, nid, vy = 1, minmaf = 0.001, ld = NULL, ldobj = NULL, ldobjlist = NULL, ldobjfiles = NULL, ldobjdir = NULL, gwasglue2 = FALSE, nthreads = 1 )generate_gwas_ss( params, nid, vy = 1, minmaf = 0.001, ld = NULL, ldobj = NULL, ldobjlist = NULL, ldobjfiles = NULL, ldobjdir = NULL, gwasglue2 = FALSE, nthreads = 1 )
params |
Output from |
nid |
sample size |
vy |
Variance of trait |
minmaf |
minimum allowed maf. default=0.01 to prevent instability |
ld |
LD correlation matrix. Must be same dimension as params |
ldobj |
LD objects (e.g. see test_ldobj) |
ldobjlist |
List of LD objects |
ldobjfiles |
Array of filenames containing LD object files (e.g. see |
ldobjdir |
Directory containing output from |
gwasglue2 |
Logical, creates a gwasglue2 SummarySet when @param ldobj is not NULL (default FALSE). |
nthreads |
Number of threads (can be slow for complete GWAS and large LD regions) |
Updated params
Determines SE and generates effect estimates given input parameters
generate_gwas_ss_1(params, nid, vy = 1, minmaf = 0.01, r = NULL)generate_gwas_ss_1(params, nid, vy = 1, minmaf = 0.01, r = NULL)
params |
Output from |
nid |
sample size |
vy |
Variance of trait |
minmaf |
minimum allowed maf. default=0.01 to prevent instability |
r |
LD correlation matrix. If NULL (default) then creates LD unaware sampling errors |
list of data frames
Creates a set of ldobj files, each corresponding to a single independent LD region from the reference panel. It also generates a map file.
generate_ldobj( outdir, bfile, regions, plink_bin = genetics.binaRies::get_plink_binary(), nthreads = 1 )generate_ldobj( outdir, bfile, regions, plink_bin = genetics.binaRies::get_plink_binary(), nthreads = 1 )
outdir |
Directory in which to store the ldobj files |
bfile |
Binary plink dataset |
regions |
A data frame containing the independent regions (see |
plink_bin |
Plink executable. Default=genetics.binaRies::get_plink_binary() |
nthreads |
How many threads. Default=1 |
map file
Get effs for two traits and make dat format
get_effs(x, y, g, xname = "X", yname = "Y")get_effs(x, y, g, xname = "X", yname = "Y")
x |
Vector of exposure phenotype |
y |
Vector of outcome phenotype |
g |
Matrix of genotypes |
xname |
xname |
yname |
yname |
Data frame
Get LD matrix for a specified region from bfile reference panel
get_ld( chr, from, to, bfile, plink_bin = genetics.binaRies::get_plink_binary(), nref = NULL )get_ld( chr, from, to, bfile, plink_bin = genetics.binaRies::get_plink_binary(), nref = NULL )
chr |
Chromosome |
from |
from bp |
to |
to bp |
bfile |
LD reference panel |
plink_bin |
Plink binary default=genetics.binaRies::get_plink_binary() |
nref |
Sample size of reference panel |
List of LD matrix and map info including MAF
Estimate the allele frequency in population from case/control summary data
get_population_allele_frequency(af, prop, odds_ratio, prevalence)get_population_allele_frequency(af, prop, odds_ratio, prevalence)
af |
Effect allele frequency (or MAF) |
prop |
Proportion of samples that are cases |
odds_ratio |
Odds ratio |
prevalence |
Population disease prevalence |
Population allele frequency
Determine regions from LD file
get_regions_from_ldobjdir(ldobjdir)get_regions_from_ldobjdir(ldobjdir)
ldobjdir |
Directory containing output from |
Data frame
Perform association of many SNPs against phenotype
gwas(y, g, logistic = FALSE)gwas(y, g, logistic = FALSE)
y |
Vector of phenotypes |
g |
Matrix of genotypes |
logistic |
Use logistic regression (much slower)? Default=FALSE |
Data frame
Translate risk from liability to probability scale
gx_to_gp(gx, h2x, prev)gx_to_gp(gx, h2x, prev)
gx |
score on liability scale |
h2x |
Disease heritability on liability scale |
prev |
Prevalence |
Vector of disease probabilities
Estimate haplotype frequencies for two loci
hap_freqs(r, p1, p2)hap_freqs(r, p1, p2)
r |
Required LD r |
p1 |
Freq 1 |
p2 |
Freq 2 |
vector
Choose initial parameters for direct effects on X and Y
init_parameters( nsnp_x, var_gx.x, var_x.y, var_gx.y = 0, nsnp_y = 0, var_gy.y = 0, mu_gx.y = 0, var_gy.x = 0, mu_gy.x = 0, prop_gy.x = 1, prop_gx.y = 1 )init_parameters( nsnp_x, var_gx.x, var_x.y, var_gx.y = 0, nsnp_y = 0, var_gy.y = 0, mu_gx.y = 0, var_gy.x = 0, mu_gy.x = 0, prop_gy.x = 1, prop_gx.y = 1 )
nsnp_x |
nsnp_x |
var_gx.x |
var_gx.x |
var_x.y |
var_x.y |
var_gx.y |
var_gx.y=0 |
nsnp_y |
nsnp_y=0 |
var_gy.y |
var_gy.y=0 |
mu_gx.y |
mu_gx.y=0 |
var_gy.x |
var_gy.x=0 |
mu_gy.x |
mu_gy.x=0 |
prop_gy.x |
prop_gy.x=1 |
prop_gx.y |
prop_gx.y=1 |
List of model parameters
Taken from Berisa and Pickrell (2016) - https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4731402/ - a list of independent LD regions for Africans, Europeans and Asians. All in build hg19
ldetectldetect
A data frame with 5731 rows and 4 columns
chromosome
start position
stop position
Population
This uses equation 10 in Genetic Epidemiology 36 : 214–224 (2012)
lor_to_rsq(b, af, ncase, ncontrol, prevalence, model = "logit")lor_to_rsq(b, af, ncase, ncontrol, prevalence, model = "logit")
b |
Log odds ratio |
af |
allele frequency |
ncase |
Number of cases |
ncontrol |
number of controls |
prevalence |
prevalence |
model |
Is the effect size estiamted in "logit" (default) or "probit" model |
Rsq
Create genotype matrix
make_geno(nid, nsnp, af)make_geno(nid, nsnp, af)
nid |
Number of samples |
nsnp |
Number of SNPs |
af |
Allele frequency of all SNPs (all SNPs have the same allele frequency) |
Matrix of genotypes, rows = individuals, columns = snps
Take several exposures and one outcome and make the data required for multivariable MR
make_mvdat(exposures, y, g)make_mvdat(exposures, y, g)
exposures |
List of exposure vectors |
y |
Vector of outcomes |
g |
Matrix of genotypes |
mv_harmonise_data output
I need a function that takes some variables, and an effect size for each of them, and this creates a new variable where y = Xb + e, where X is the matrix of independent variables, b is their effect sizes, and e is a noise term. By default it's easier if b relates to the variance explained by each vector in X, and e has variance of 1 - sum(b).
make_phen(effs, indep, vy = 1, vx = rep(1, length(effs)), my = 0)make_phen(effs, indep, vy = 1, vx = rep(1, length(effs)), my = 0)
effs |
An array of effect sizes that the variables in indep have on the variable that you are simulating |
indep |
A matrix of variables, rows = Samples and columns = variables that have an influence on the variable that you are simulating |
vy |
What variance the output should have. Default = 1, meaning effs relate to the signed rsq of the influence of the indep variables on the outcome |
vx |
What variance the indep variables should have. Default is to set to 1, meaning that the effects are the signed variance explained |
my |
mean value of y to be output |
Vector of y values
gwas into harmonised dat formatOrganise outputs from gwas into harmonised dat format
merge_exp_out(gwasx, gwasy, xname = "X", yname = "Y")merge_exp_out(gwasx, gwasy, xname = "X", yname = "Y")
gwasx |
Output from |
gwasy |
Output from |
xname |
exposure name |
yname |
outcome name |
data frame
Scale variable to have range between 0 and 1
range01(x)range01(x)
x |
Vector |
Vector
Read in LD objects into list
read_ldobjdir(ldobjdir, nthreads = 1)read_ldobjdir(ldobjdir, nthreads = 1)
ldobjdir |
Directory containing output from |
nthreads |
Number of threads. Default=1 |
List of ldobj
Can use simple method or by pivoting around intercept
recode_dat(dat, method = "intercept")recode_dat(dat, method = "intercept")
dat |
Output from get_effs |
method |
Default 'intercept'. Alternatively can specify 'simple' |
Data frame
Tries to avoid issue of recoding by finding intercept and pivoting negative g-x associations around intercept
recode_dat_intercept(dat)recode_dat_intercept(dat)
dat |
Output from get_effs |
Data frame
Simple recoding to have every effect on x positive
recode_dat_simple(dat)recode_dat_simple(dat)
dat |
Output from get_effs |
Data frame
Plot liability vs probability disease risk
risk_cross_plot( x, o, xlab = "Values (low to high)", ylab = "", title = "", xname = "GRS", oname = "Disease" )risk_cross_plot( x, o, xlab = "Values (low to high)", ylab = "", title = "", xname = "GRS", oname = "Disease" )
x |
Disease risk on liability scale |
o |
Disease risk on probability scale |
xlab |
default="Values (low to high)" |
ylab |
default="" |
title |
default="" |
xname |
Name of liability. Default="GRS" |
oname |
Name of disease. Default="Disease" |
ggplot
Compares the liability and probability of disease under two scenarios - where all SNPs are known, and where only some SNPs are known
risk_simulation(G, eff, prevalence, prop_discovered)risk_simulation(G, eff, prevalence, prop_discovered)
G |
Matrix of genotypes |
eff |
SNP effects on liability scale |
prevalence |
Disease prevalence |
prop_discovered |
Proportion of SNPs discovered |
Data frame
Sample beta values given standard errors
sample_beta(beta, se, r = NULL, af = NULL)sample_beta(beta, se, r = NULL, af = NULL)
beta |
array of beta values - i.e. the true coefficients. If using the correlation matrix r then this should be an LD-aware set of expected beta values. |
se |
array of se values |
r |
matrix of LD correlations amongst the SNPs. If NULL (default) then assumes no LD |
af |
array of allele frequencies. Must be non-null if r is non-null. |
array of beta hats
Sample the actual effects based on initial parameters
sample_system_effects(parameters)sample_system_effects(parameters)
parameters |
Output from |
List of effect parameters
Simulate genotypes from haplotypes
simulate_geno(nid, r, p1, p2)simulate_geno(nid, r, p1, p2)
nid |
Number of samples |
r |
Desired LD r |
p1 |
Freq 1 |
p2 |
Freq 2 |
Matrix
Simulate haplotypes of two loci
simulate_haplotypes(nid, r, p1, p2)simulate_haplotypes(nid, r, p1, p2)
nid |
Number of samples |
r |
Desired LD r |
p1 |
Freq 1 |
p2 |
Freq 2 |
Matrix
Simulate individual level data from initial parameters
simulate_population(parameters, nid)simulate_population(parameters, nid)
parameters |
Output from |
nid |
Sample size to generate |
List of matrices and vectors that represent individual level data
Allows arbitrary sample overlap
summary_set( beta_gx, beta_gy, af, n_gx, n_gy, n_overlap, cor_xy, prev_y = NA, sigma_x = 1, sigma_y = 1 )summary_set( beta_gx, beta_gy, af, n_gx, n_gy, n_overlap, cor_xy, prev_y = NA, sigma_x = 1, sigma_y = 1 )
beta_gx |
Array of true effects on x |
beta_gy |
Array of true effects on y |
af |
Array of effect allele frequencies |
n_gx |
Sample size of g-x association |
n_gy |
Sample size of g-y association |
n_overlap |
Number of overlapping samples |
cor_xy |
Observational correlation between x and y |
prev_y |
Disease prevalence of y. Default = NA (in which case treated as continuous) |
sigma_x |
SD of x. Default=1 |
sigma_y |
SD of y. Default=1 |
Data frame of summary statistics for x and y
Create test LD object
test_ldobj(nsnp, chunksize)test_ldobj(nsnp, chunksize)
nsnp |
Number of SNPs |
chunksize |
Chunksize for splitting |
list of chunks, which each contain map and LD matrix
Apply MR tests to system
test_system(ss, id = "test")test_system(ss, id = "test")
ss |
Output from create_syste, |
id |
string denoting simulation ID |
List
Convert continuous trait to binary
y_to_binary(y, prevalence = NULL, threshold = NULL)y_to_binary(y, prevalence = NULL, threshold = NULL)
y |
Phenotype vector |
prevalence |
Disease prevalence. Default = NULL |
threshold |
Disease threshold Default = NULL |
Vector of binary trait