In mdfortune/simGWAS: Simulating the expected output from a GWAS with a given causal model

Load the library

library(simGWAS)

simGWAS needs some reference haplotype frequencies from control subjcets. These can be found by taking phased haplotypes from public data sources, or by phasing genotype data you may already have, for example using snphap.

For the purpose of this vignette, we will simulate some reference haplotypes. The final format is a data.frame with n columns of 0s and 1s indicating alleles at each of n SNPs, and collections of alleles in a row being a haplotype. A final column, named "Probability", contains the fractional frequency of each haplotype. Note that haplotypes need not be unique, you can have one row per haplotype in a sample, and Probability set to 1/[number of haplotypes] = 1/(2*[number of samples]). The object we are creating will be called freq.

nsnps <- 100
nhaps <- 1000
lag <- 5 # genotypes are correlated between neighbouring variants
maf <- runif(nsnps+lag,0.05,0.5) # common SNPs
laghaps <- do.call("cbind", lapply(maf, function(f) rbinom(nhaps,1,f)))
haps <- laghaps[,1:nsnps]
for(j in 1:lag) 
    haps <- haps + laghaps[,(1:nsnps)+j]
haps <- round(haps/(lag+1))

snps <- colnames(haps) <- paste0("s",1:nsnps)
freq <- as.data.frame(haps+1)
freq$Probability <- 1/nrow(freq)
sum(freq$Probability)

Next, we need to specify the causal variants, and their effects on disease, as odds ratios. We create a vector CV with snp names that are a subset of column names in freq and a vector of odds ratios. In our simulated data, we pick two causal variants at random, with odds ratios of 1.4 and 1.2.

CV=sample(snps,2)
g1 <- c(1.4,1.2)

Now we simulate the results of a GWAS. There are two key functions, makeGenoProbList calculates genotype probabilities at each SNP conditional on genotypes at the causal variants, then est_statistic generates the vector of Z scores across all SNPs, conditional on the causal variants and their odds ratios.

FP <- make_GenoProbList(snps=snps,W=CV,freq=freq)
z <- est_statistic(N0=3000, # number of controls
              N1=2000, # number of cases
              snps=snps, # column names in freq of SNPs for which Z scores should be generated
              W=CV, # causal variants, subset of snps
              gamma1=g1, # odds ratios
              freq=freq, # reference haplotypes
              GenoProbList=FP) # FP above

Ignoring spacing, this would produce results like, with red lines indicating where the causal variants are.

plot(1:nsnps,z); abline(v=which(snps %in% CV),col="red"); abline(h=0)

mdfortune/simGWAS documentation built on May 22, 2019, 3:25 p.m.