chr1swallace

##' Fit a specific two Guassian mixture distribution
##'
##' Assumes Z ~ N(0,1) with prob pi0,
##' Z ~ N(0,1+sigma^2) with prob 1-pi0
##'
##' Aims to estimate sigma^2 and pi0.
##'
##' @title fit.em
##' @param Z numeric vector of observed data
##' @param s2 initial value for sigma^2

chr1swallace

# replace "gantt.csv" with path to wherever your .csv file is stored
x <- read.table("gantt.csv",sep=",",as.is=TRUE, header=TRUE) # read data

## look at what we read, check it makes sense
head(x)

## make sure dates are recognised as such
x$Start <- as.Date(x$Start)
x$End <- as.Date(x$End)

chr1swallace

## download https://www.t1dbase.org/downloads/protected_data/GWAS_Data/hg19_gwas_t1d_barrett_4_18_0.gz
## you will need to create a login!
## Then...
x <- read.table("hg19_gwas_t1d_barrett_4_18_0.gz",sep="\t",skip=1,as.is=TRUE)
colnames(x) <- c("snp","Chr","Pos","PValue","OR","kk")

## libraries
library(ggplot2)
library(cowplot) # nicer ggplot theme

chr1swallace

coloc.bayes.3t <- function(df1,snps=setdiff(colnames(df1),response),response="Y",priors=list(rep(1,14)),r2.trim=0.99,thr=0.005,quiet=TRUE) {
    #we consider all models which contain at most 1 snp for each of the three traits
    snps <- unique(snps)
    n.orig <- length(snps)
    if(n.orig<2)
        return(1)
    prep <- prepare.df(df1, snps, r2.trim=r2.trim, dataset=1, quiet=quiet)
    df1 <- prep$df
    snps <- prep$snps
    

chr1swallace

library(MASS)
data<-mvrnorm(n=100, mu=c(1,2), Sigma=matrix(c(1,0.5,0.5,1),2))

colnames(data) <- c("x","y")
var(data)
cor(data)

data[,"x"] <- data[,"x"] + 3
var(data) # unchanged
cor(data) # unchanged

chr1swallace

library(snpStats)
  best <- cond.best(X,Y, family=family)
  while(length(newbest <- cond.best(X, Y, best, family=family)))
    best <- c(best,newbest)

cond.best <- function(X,Y,best=NULL,p.thr=1e-6, ...) {
  if(is.null(best)) {
    cond <- snp.rhs.tests(Y ~ 1, snp.data=X, data=data.frame(Y=Y,
row.names=rownames(X)), ...)
    p.thr <- 1

chr1swallace

library(ggdendro)
library(ggplot2)
library(reshape)
library(grid)
library(gtable)

## colours, generated by
## library(RColorBrewer)
## rev(brewer.pal(11,name="RdYlBu"))
my.colours <- c("#313695", "#4575B4", "#74ADD1", "#ABD9E9", "#E0F3F8", "#FFFFBF", 

chr1swallace

## sharing parameter, largish, so we expect to see a big enough difference to check if working
S <- 100
prior.bin.fn <- function(nT1,s=100,mT1=3) {
    dbinom(nT1,size=s,prob=mT1/s)/choose(s,nT1)
}

bothpp <- function(data) {
    ss1 <- strsplit(data$t1.str,"%")
    ss2 <- strsplit(data$t2.str,"%")
    data$overlap <- sapply(1:nrow(data), function(i)

chr1swallace

## sharing parameter, largish, so we expect to see a big enough difference to check if working
S <- 100
prior.bin.fn <- function(nT1,s=100,mT1=3) {
    dbinom(nT1,size=s,prob=mT1/s)/choose(s,nT1)
}

bothpp <- function(data) {
    ss1 <- strsplit(data$t1.str,"%")
    ss2 <- strsplit(data$t2.str,"%")
    data$overlap <- sapply(1:nrow(data), function(i)

chr1swallace

library(data.table)
library(parallel)
library(annotSnpStats)
library(snpStats)
library(GUESSFM)

##
inv.logit.fn <-function(x) return(exp(x)/(1+exp(x)))
##