explodecomputer

library(TwoSampleMR)
library(MRInstruments)
library(ggplot2)
library(dplyr)
library(reshape2)
library(classInt)


calculate_ios <- function(exposure_dat, outcome_dat)
{

explodecomputer

library(dplyr)

# X causes DEATH
# P causes DEATH
# AGE causes DEATH
# G causes X
# AGE causes Y
# - IF cases and controls are AGE matched then NO SURVIVAL BIAS

n <- 1000000

explodecomputer

library(data.table)

arguments <- commandArgs(T)

phenfile <- arguments[1]
covfile <- arguments[2]
outfile <- arguments[3]
lmodel <- arguments[4]

# expect lmodel to be "logistic" or "linear"

explodecomputer

#!/bin/bash

# Load plink2 module
module add apps/plink-1.90b3v

# Location of the biobank genotype data
datadir="/panfs/panasas01/shared-biobank/data/"

# Location of the user supplied phenotype file
phenfile="phenfile.txt"

explodecomputer

n <- 1000
nsnp <- 10

r <- matrix(0.8, nsnp, nsnp)
diag(r) <- 1

x1 <- rmvnorm(n, rep(0, nsnp), sigma=r)
x2 <- rmvnorm(n, rep(0, nsnp), sigma=r)

x1[1:500, nsnp] <- x2[1:500, 1]

explodecomputer

# source("https://bioconductor.org/biocLite.R")
# biocLite("KEGGREST")
# biocLite("org.Hs.eg.db")
library(biomaRt)
library(KEGGREST)
library(org.Hs.eg.db)

# example snp list
snp <- c("rs1007648", "rs112275031", "rs1124639", "rs11731570", "rs11743154", "rs12439909")

explodecomputer

n <- 10000

gm <- rbinom(n, 2, 0.5)
gp <- rbinom(n, 2, 0.5)

# approximate child genotypes
gc <- round((gm + gp) / 2)

cor(gm, gc)^2
cor(gp, gc)^2

explodecomputer

---
title: "Bivariate heritabilities greater than 1"
date: "`r format(Sys.time(), '%d %B %Y')`"
output: pdf_document
---

the model -

$$
\begin{aligned}

explodecomputer

1M = 16mb

a <- read.table('results_9.gz', he=T)

snplist <- as.character(unique(a$SNP))
snplist <- data.frame(SNP=snplist, snpid=1:length(snplist), stringsAsFactors=FALSE)

cpglist <- as.character(unique(a$PHENOTYPE))
cpglist <- data.frame(CPG=cpglist, cpgid=1:length(cpglist), stringsAsFactors=FALSE)

explodecomputer

library(dplyr)

# Generate some LD blocks for 1000000 SNPs
# Some blocks have more SNPs than others
ld <- sample(1:1000, 1000000, replace=TRUE, prob=runif(1000))
range(table(ld))

# Generate some xsq stats for each SNP
xsq <- rchisq(1000000, df=1)