| # Sample size of PAPA replication | |
| n <- 149 | |
| # Pvals from first 14 CpGs in PAPA replication | |
| pval <- c(9.6e-6, 7.2e-5, 2.8e-6, 3.7e-5, 4.8e-4, 6.2e-3, 4.6e-6, 1.3e-2, 2.6e-3, 3.8e-4, 3.2e-5, 9.0e-4, 8.8e-5, 8.3e-5) | |
| # Estimate F statistics | |
| qval <- qf(pval, 1, n-1, low=F) |
| calc.proprare <- function(nsnp, shape1, eff.func, nid) | |
| { | |
| # Generate allele frequency distribution | |
| mafs <- rbeta(nsnp, shape1, 1) / 2 | |
| # Limit min maf to be function of hypothetical sample size | |
| mafs <- pmax(mafs, 1/(nid*2)) | |
| # Create distribution of effect sizes | |
| eff <- eff.func(mafs) |
A lot of servers use Linux and we need to use a server to run a few specific programmes so let's get familiar with Linux.
First of all, to connect to the server we will use a SSH connection. The idea is that we log in to the remote server from our local computer, so when we are running anything on the server we are using its hardware and infrastructure and our local computer is being used to simply view what is happening remotely and to input commands.
| #!/bin/bash | |
| # name of binary plink filename (excluding .bed/.bim/.fam suffix) | |
| plinkfile="" | |
| # output filename | |
| outname="" | |
| # number of causal variants | |
| nvar=1000 |
| # Simulate the smoking thing | |
| # SNP influences smoking status | |
| # Smoking status influences income | |
| # Test for association of SNP on income in smokers and non smokers | |
| nsim <- 100 | |
| res <- matrix(0, nsim, 4) | |
| n <- 100000 | |
| for(i in 1:nsim) |
| # Null is true | |
| nsim <- 1000 | |
| p <- matrix(0, nsim, 3) | |
| n <- 1000 | |
| g1 <- rbinom(n, 2, 0.2) | |
| g2 <- rbinom(n, 2, 0.3) | |
| g3 <- rbinom(n, 2, 0.4) | |
| g <- cbind(g1, g2, g3) |
| nsnp <- 48 | |
| n <- 1000 | |
| mothers1 <- matrix(rbinom(n*nsnp*2, 1, 0.5), n) | |
| kids1 <- matrix(rbinom(n*nsnp, 1, 0.5), n) | |
| kids2 <- mothers1[, seq(1, 2*nsnp, 2)] | |
| kids <- kids1 + kids2 |
| ## ---- functions ---- | |
| getAlleleFreq <- function(n, ahat, Vp, pval) | |
| { | |
| stopifnot(ahat != 0) | |
| Fval <- qf(pval, 1, n-1, low=F) | |
| p <- 0.5 + c(1,-1) * -2 * sqrt(ahat^2 - 2 * Fval * Vp / (n - 1 - Fval)) / (4 * ahat) | |
| return(p[ifelse(sign(ahat) == 1, 1, 2)]) | |
| } |
| #!/bin/bash | |
| #PBS -N test | |
| #PBS -o test-output | |
| #PBS -e test-error | |
| #PBS -t 1-500 | |
| #PBS -l walltime=5:00:00 | |
| #PBS -l nodes=1:ppn=1 | |
| #PBS -S /bin/bash |
| # Assume that we have 2 traits and large sample size | |
| # Each trait is a function of g and e | |
| # g1 and g2 are correlated | |
| # e1 and e2 are correlated | |
| # g and e are independent | |
| # y1 and y2 are correlated due to correlations in g and e combined | |
| n <- 100000 | |
| g1 <- rnorm(n) | |
| e1 <- rnorm(n) |