explodecomputer

m <- data.frame(matrix(rnorm(1000*20), 1000, 20))

y = rnorm(1000)

performScan <- function(y, dat)
{
	n <- ncol(dat)
	res <- array(0, n)
	for(i in 1:n)
	{

explodecomputer

# Create some SNPs:

px <- 0.5
py <- 0.3
pz <- 0.2

snpx <- rbinom(1000, 2, px)
snpy <- rbinom(1000, 2, py)
snpz <- rbinom(1000, 2, pz)

explodecomputer

#!/bin/bash

#PBS -N sleepgwas
#PBS -o sleepsnptest-output
#PBS -e sleepsnptest-error
#PBS -t 1-22
#PBS -l walltime=24:00:00
#PBS -l nodes=1:ppn=2
#PBS -S /bin/bash

explodecomputer

# SNP-disease association - need effect size and standard error, can use log(OR) instead of effect size
y_or <- 1.244
y_upperconf <- 1.375
y_lowerconf <- 1.126

y_eff <- log(y_or)
y_se <- (log(y_upperconf) - log(y_lowerconf)) / (2*1.96)

# Alternative method from Shinn 2000, but uncertain what scale this transforms to
# y_eff <- log(y_or) / 1.81

explodecomputer

#!/bin/bash

# name of binary plink filename (excluding .bed/.bim/.fam suffix)
plinkfile=""

# name of phenotype file in plink format (i.e. col1=fid, col2=iid, col3=phenotype
phenfile=""

# make up a name for the grm for all individuals
allfile=""

explodecomputer

#' Read binary GRM files into R
#'
#' @param rootname
#' @export
#' @return List of GRM and id data frames
readGRM <- function(rootname)
{
	bin.file.name <- paste(rootname, ".grm.bin", sep="")
	n.file.name <- paste(rootname, ".grm.N.bin", sep="")
	id.file.name <- paste(rootname, ".grm.id", sep="")

explodecomputer

#!/bin/bash
 
set -e


# This script extracts a list of SNPs from gzipped .gen files that are split into 23 chromosomes
# and outputs as a single gzipped .gen file containing all extracted SNPs that it can find
# To run use the command
#
# $ ./extractSnpsGenfiles.sh \

explodecomputer

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

explodecomputer

#!/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

explodecomputer

## ---- 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)])
}