-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathDYexperiments.Rmd
132 lines (112 loc) · 4.74 KB
/
DYexperiments.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
---
title: "DYexperiments"
author: "Daisy Yu"
date: "11/8/2020"
output: pdf_document
---
## Case control simulation with population as a confounder
```{r}
library(dplyr)
library(ggplot2)
simSNPCov_popl=function(n,Beta=log(rf(K,m,m)),MAF1,MAF2,beta_s){
ncase <- ncon <- n/2
conX <- matrix(NA,ncol=K+1,nrow=ncon)
caseX <- matrix(NA,ncol=K+1,nrow=ncase)
# Simulate population status first (population is the K+1 column)
f_0 <- f_1 <- 0.5
conX[,K+1] <- rbinom(ncon,1,f_0)
caseX[,K+1] <- sample(0:1,size=ncase,replace=TRUE,prob=c(f_0,f_1*exp(beta_s)))
# Simulate SNVs conditional on population status
for (i in 1:K){
maf0 <- MAF1[i]
maf1 <- MAF2[i]
beta <- Beta[i]
## control
con_pop0 <- which(conX[,K+1]==0)
con_pop1 <- which(conX[,K+1]==1)
conX[con_pop0,i] <- rbinom(length(con_pop0),size=2,prob=maf0)
conX[con_pop1,i] <- rbinom(length(con_pop1),size=2,prob=maf1)
## case
case_pop0 <- which(caseX[,K+1]==0)
case_pop1 <- which(caseX[,K+1]==1)
pp0 <- c((1-maf0)^2,2*maf0*(1-maf0)*exp(beta),maf0^2*exp(2*beta))
pp1 <- c((1-maf1)^2,2*maf1*(1-maf1)*exp(beta),maf1^2*exp(2*beta))
caseX[case_pop0,i] <- sample(0:2,size=length(case_pop0),replace=TRUE,prob=pp0)
caseX[case_pop1,i] <- sample(0:2,size=length(case_pop1),replace=TRUE,prob=pp1)
}
X <- rbind(caseX,conX)
colnames(X) <- c(paste0("X",1:K),"Population");rownames(X) <- NULL
case <- c(rep(1,ncase),rep(0,ncon))
return (data.frame(X,case))
}
```
### Sample MAF from the joint SFS of CEU and YRI populations.
```{r}
set.seed(8888)
load("~/Desktop/Simulation/1000GenomesProject/Input_dataset_1000GP_CEU.RData")
load("~/Desktop/Simulation/1000GenomesProject/Input_dataset_1000GP_YRI.RData")
m=4
K=50
index <- intersect(colnames(X_CEU),colnames(X_YRI))
X_CEU <- X_CEU[,index]
X_YRI <- X_YRI[,index]
MAF_CEU <- colSums(X_CEU)/(2*dim(X_CEU)[1])
MAF_YRI <- colSums(X_YRI)/(2*dim(X_YRI)[1])
Snp_data <- data.frame(CEU=MAF_CEU,YRI=MAF_YRI)
Snp_data <- Snp_data[-which(Snp_data$CEU<=0.05 & Snp_data$YRI<=0.05),]
maf <- Snp_data[sample(nrow(Snp_data),K,replace=F),]
Beta <- log(rf(K,m,m))
sim_data <- simSNPCov_popl(200,Beta,MAF1 = maf[,1],MAF2 = maf[,2],beta_s=1)
betahat <- c() # with popl
betahat0 <- c() # without popl
for (i in 1:K) {
XX <- sim_data %>% dplyr::select(c(paste0("X",i),Population,case)) %>% dplyr::rename(X=paste0("X",i))
betahat[i] <- coef(glm(case~.,data=XX,family=binomial(),maxit=100))[[2]]
XX0 <- sim_data %>% dplyr::select(c(paste0("X",i),case)) %>% dplyr::rename(X=paste0("X",i))
betahat0[i] <- coef(glm(case~.,data=XX0,family=binomial(),maxit=100))[[2]]
}
Confounder=c()
for (i in 1:K){
if (abs(betahat[i]-betahat0[i])/abs(betahat[i]) > 0.1) {Confounder[i]=T}
else {Confounder[i]=F}
}
sum(Confounder)
maf <- cbind(maf,Confounder)
X1 <- data.frame(value=betahat)
X2 <- data.frame(value=betahat0)
X1$var <- "with_popl"
X2$var <- "without_popl"
data <- rbind(X1, X2)
ggplot(data, aes(value, fill = var)) + geom_density(alpha = 0.2) + xlim(c(-5,5))
```
26 out of 50 coefficients are changed by more than 10\% when population is excluded. The distribution of estimated coefficients doesn't change a lot with or without population included in the model, so we should get a similar $m$.
### What if we only restrict to SNPs whose MAF differs by more than 0.3 between the two populations ?
```{r}
set.seed(8888)
Snp_data <- Snp_data[which(abs(Snp_data$CEU-Snp_data$YRI) >= 0.3),]
maf <- Snp_data[sample(nrow(Snp_data),K,replace=F),]
Beta <- log(rf(K,m,m))
sim_data <- simSNPCov_popl(200,Beta,MAF1=maf[,1],MAF2=maf[,2],beta_s=1)
betahat <- c() # with popl
betahat0 <- c() # without popl
for (i in 1:K) {
XX <- sim_data %>% dplyr::select(c(paste0("X",i),Population,case)) %>% dplyr::rename(X=paste0("X",i))
betahat[i] <- coef(glm(case~.,data=XX,family=binomial(),maxit=100))[[2]]
XX0 <- sim_data %>% dplyr::select(c(paste0("X",i),case)) %>% dplyr::rename(X=paste0("X",i))
betahat0[i] <- coef(glm(case~.,data=XX0,family=binomial(),maxit=100))[[2]]
}
Confounder=c()
for (i in 1:K){
if (abs(betahat[i]-betahat0[i])/abs(betahat[i]) > 0.1) {Confounder[i]=T}
else {Confounder[i]=F}
}
sum(Confounder)
maf <- cbind(maf,Confounder)
X1 <- data.frame(value=betahat)
X2 <- data.frame(value=betahat0)
X1$var <- "with_popl"
X2$var <- "without_popl"
data <- rbind(X1, X2)
ggplot(data, aes(value, fill = var)) + geom_density(alpha = 0.2) + xlim(c(-5,5))
```
Now 42 out of 50 coefficients are changed by more than 10\% when population is excluded. The distribution of estimated coefficients with population included in the model has a heavier tail than the distribution without population (i.e. coefficients are more extreme when we include population as a confounder). In this case, I guess we should get a much smaller $m$ when we include population.