explodecomputer

for i in {1..22}
do
plink --bfile /mnt/storage/private/mrcieu/research/mr-eve/vcf-reference-datasets/1000g_filtered/data_maf0.01_rs --chr $i --make-bed --out ld/ref/data_$i
done

for i in {1..22}
do
bcftools query -r $i -f '%ID\t%POS\t%REF\t%ALT\t%B\t%SE\n' ../data.bcf | awk 'BEGIN {print "rsID\tpos\tA0\tA1\tZ"; OFS="\t"}; { print $1, $2, $3, $4, $5/$6}' > ld/ref/z_2_$i.txt
done

explodecomputer

library(TwoSampleMR)
library(RadialMR)
library(MendelianRandomization)

# Read in example data
a <- extract_instruments(2)
b <- extract_outcome_data(a$SNP, 7)
ab <- harmonise_data(a,b)
ab$SNP <- as.character(ab$SNP)
# IVW in TwoSampleMR package

explodecomputer

#!/bin/bash

# Download 1000 genomes

wget ALL.chr1.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz
wget ALL.chr1.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz.tbi
wget ALL.chr10.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz
wget ALL.chr10.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz.tbi
wget ALL.chr11.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz
wget ALL.chr11.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz.tbi

explodecomputer

#===================================#
#          GWAS Practical           #
#      iBSC Genomic Medicine        #
#       Lavinia Paternoster         #
#           08/11/2018              #
#===================================#


# Run the GWAS without PCs but using the cleaned genetic data

explodecomputer

# Get the list of arguments passed to R
args <- commandArgs(T)

# Get the first argument, make it numeric because by default it is read as a string
job_id <- as.numeric(args[1])

# Print the argument
message("job number ", job_id)

explodecomputer

n <- 5
a <- combn(c(1:(2*n)), 2)
b <- combn(1:ncol(a), n)
l <- list()
m <- list()
j <- 1
for(i in 1:ncol(b))
{
	x <- a[, b[,i]]
	if(!any(duplicated(c(x))))

explodecomputer

# How does adjusted r-square do with overfitting?

rm(list=ls())
set.seed(101)

library(ggplot2)
library(tidyr)

y <- rnorm(1000)
x <- matrix(rnorm(1000 * 1000), 1000, 1000)

explodecomputer

dat_3a$rsq.exposure <- with(dat_3a, 2 * beta.exposure^2 * eaf.exposure * (1-eaf.exposure))
dat_3a$rsq.outcome <- with(dat_3a, 2 * beta.outcome^2 * eaf.outcome * (1-eaf.outcome))
dat_3a$steiger_dir <- dat_3a$rsq.exposure > dat_3a$rsq.outcome

dat_3b$rsq.exposure <- with(dat_3b, 2 * beta.exposure^2 * eaf.exposure * (1-eaf.exposure))
dat_3b$rsq.outcome <- with(dat_3b, 2 * beta.outcome^2 * eaf.outcome * (1-eaf.outcome))
dat_3b$steiger_dir <- dat_3b$rsq.exposure > dat_3b$rsq.outcome

explodecomputer

---
title: Calculating HWE
---

Choose a set of parameters. For each allele frequency and sample size, run 100 simulations

```{r}
library(tidyr)
library(ggplot2)
library(dplyr)

explodecomputer

library(devtools)
install_github("explodecomputer/simulateGP")
library(simulateGP)
library(ggplot2)

param <- rbind(
	expand.grid(
		ncase=100,
		ncontrol=100,
		prev=0.001,