Newer
Older
Yunfeng Huang
committed
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
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
#######################################################
######### UKBB lipids PTV-PRS interaction #############
######### Author: Yunfeng Huang #############
#######################################################
library(data.table)
library(dplyr)
library(MASS)
# Read phenotype data
lipids_pheno <- as.data.frame(fread("~/lipids_heiko/LoF_UKB/2021/lipids_lofs_pheno_300k.txt"))
# Test for PTV-PRS interaction through robust linear model fitting and extract BIC as well as interaction model coefficients (beta and t-value)
#HDL
genes <- colnames(lipids_pheno)[8:37]
lof_prs_results_hdl <- as.data.frame(matrix(NA,30,12))
for(i in 1:length(genes)){
PRS <- fread(paste("/camhpc/home/cchen11/lipid_interaction/hdl.pst_eff_a1_b05_phiauto_prscs.no_",genes[i],".totalscore",sep = ""),sep = " ")
data_sub <- merge(select(lipids_pheno,eid = eid.x,hdl_res,starts_with(genes[i])),select(PRS,eid = IID,totalscore = TOTALSCORE_AVG),by = "eid")
data_sub$prs_std <- scale(data_sub$totalscore)[,1]
data_sub <- na.omit(data_sub)
colnames(data_sub) <- c("eid","hdl_res","ptv_gene","prs_total","prs_std")
c1 <- sum(data_sub$ptv_gene>0)
model1 <- rlm(hdl_res~ptv_gene,data = data_sub)
model2 <- rlm(hdl_res~prs_std,data = data_sub)
model3 <- rlm(hdl_res~ptv_gene+prs_std,data = data_sub)
model4 <- rlm(hdl_res~ptv_gene+prs_std+ptv_gene:prs_std,data = data_sub)
lof_prs_results_hdl[i,1] <- genes[i]
lof_prs_results_hdl[i,2] <- BIC(model1)
lof_prs_results_hdl[i,3] <- BIC(model2)
lof_prs_results_hdl[i,4] <- BIC(model3)
lof_prs_results_hdl[i,5] <- BIC(model4)
lof_prs_results_hdl[i,6] <- summary(model4)$coeff[2,1]
lof_prs_results_hdl[i,7] <- summary(model4)$coeff[3,1]
lof_prs_results_hdl[i,8] <- summary(model4)$coeff[4,1]
lof_prs_results_hdl[i,9] <- summary(model4)$coeff[2,3]
lof_prs_results_hdl[i,10] <- summary(model4)$coeff[3,3]
lof_prs_results_hdl[i,11] <- summary(model4)$coeff[4,3]
lof_prs_results_hdl[i,12] <- c1
if(ceiling(i/10)==(i/10)){print(i)}
}
colnames(lof_prs_results_hdl) <- c("PTV_gene","model1_BIC","model2_BIC","model3_BIC","model4_BIC",
"PTV_E","PRS_E","PTV_PRS_E","PTV_t","PRS_t","PTV_PRS_t","N_carrier_PTV")
#LDL
genes <- colnames(lipids_pheno)[8:37]
lof_prs_results_ldl <- as.data.frame(matrix(NA,30,12))
for(i in 1:length(genes)){
PRS <- fread(paste("/camhpc/home/cchen11/lipid_interaction/ldl.pst_eff_a1_b05_phiauto_prscs.no_",genes[i],".totalscore",sep = ""),sep = " ")
data_sub <- merge(select(lipids_pheno,eid = eid.x,ldl_res,starts_with(genes[i])),select(PRS,eid = IID,totalscore = TOTALSCORE_AVG),by = "eid")
data_sub$prs_std <- scale(data_sub$totalscore)[,1]
data_sub <- na.omit(data_sub)
colnames(data_sub) <- c("eid","ldl_res","ptv_gene","prs_total","prs_std")
c1 <- sum(data_sub$ptv_gene>0)
model1 <- rlm(ldl_res~ptv_gene,data = data_sub)
model2 <- rlm(ldl_res~prs_std,data = data_sub)
model3 <- rlm(ldl_res~ptv_gene+prs_std,data = data_sub)
model4 <- rlm(ldl_res~ptv_gene+prs_std+ptv_gene:prs_std,data = data_sub)
lof_prs_results_ldl[i,1] <- genes[i]
lof_prs_results_ldl[i,2] <- BIC(model1)
lof_prs_results_ldl[i,3] <- BIC(model2)
lof_prs_results_ldl[i,4] <- BIC(model3)
lof_prs_results_ldl[i,5] <- BIC(model4)
lof_prs_results_ldl[i,6] <- summary(model4)$coeff[2,1]
lof_prs_results_ldl[i,7] <- summary(model4)$coeff[3,1]
lof_prs_results_ldl[i,8] <- summary(model4)$coeff[4,1]
lof_prs_results_ldl[i,9] <- summary(model4)$coeff[2,3]
lof_prs_results_ldl[i,10] <- summary(model4)$coeff[3,3]
lof_prs_results_ldl[i,11] <- summary(model4)$coeff[4,3]
lof_prs_results_ldl[i,12] <- c1
if(ceiling(i/10)==(i/10)){print(i)}
}
colnames(lof_prs_results_ldl) <- c("PTV_gene","model1_BIC","model2_BIC","model3_BIC","model4_BIC",
"PTV_E","PRS_E","PTV_PRS_E","PTV_t","PRS_t","PTV_PRS_t","N_carrier_PTV")
#TC
genes <- colnames(lipids_pheno)[8:37]
lof_prs_results_tc <- as.data.frame(matrix(NA,30,12))
for(i in 1:length(genes)){
PRS <- fread(paste("/camhpc/home/cchen11/lipid_interaction/tch.pst_eff_a1_b05_phiauto_prscs.no_",genes[i],".totalscore",sep = ""),sep = " ")
data_sub <- merge(select(lipids_pheno,eid = eid.x,tc_res,starts_with(genes[i])),select(PRS,eid = IID,totalscore = TOTALSCORE_AVG),by = "eid")
data_sub$prs_std <- scale(data_sub$totalscore)[,1]
data_sub <- na.omit(data_sub)
colnames(data_sub) <- c("eid","tc_res","ptv_gene","prs_total","prs_std")
c1 <- sum(data_sub$ptv_gene>0)
model1 <- rlm(tc_res~ptv_gene,data = data_sub)
model2 <- rlm(tc_res~prs_std,data = data_sub)
model3 <- rlm(tc_res~ptv_gene+prs_std,data = data_sub)
model4 <- rlm(tc_res~ptv_gene+prs_std+ptv_gene:prs_std,data = data_sub)
lof_prs_results_tc[i,1] <- genes[i]
lof_prs_results_tc[i,2] <- BIC(model1)
lof_prs_results_tc[i,3] <- BIC(model2)
lof_prs_results_tc[i,4] <- BIC(model3)
lof_prs_results_tc[i,5] <- BIC(model4)
lof_prs_results_tc[i,6] <- summary(model4)$coeff[2,1]
lof_prs_results_tc[i,7] <- summary(model4)$coeff[3,1]
lof_prs_results_tc[i,8] <- summary(model4)$coeff[4,1]
lof_prs_results_tc[i,9] <- summary(model4)$coeff[2,3]
lof_prs_results_tc[i,10] <- summary(model4)$coeff[3,3]
lof_prs_results_tc[i,11] <- summary(model4)$coeff[4,3]
lof_prs_results_tc[i,12] <- c1
if(ceiling(i/10)==(i/10)){print(i)}
}
colnames(lof_prs_results_tc) <- c("PTV_gene","model1_BIC","model2_BIC","model3_BIC","model4_BIC",
"PTV_E","PRS_E","PTV_PRS_E","PTV_t","PRS_t","PTV_PRS_t","N_carrier_PTV")
#TG
genes <- colnames(lipids_pheno)[8:37]
lof_prs_results_tg <- as.data.frame(matrix(NA,30,12))
for(i in 1:length(genes)){
PRS <- fread(paste("/camhpc/home/cchen11/lipid_interaction/tri.pst_eff_a1_b05_phiauto_prscs.no_",genes[i],".totalscore",sep = ""),sep = " ")
data_sub <- merge(select(lipids_pheno,eid = eid.x,tg_res,starts_with(genes[i])),select(PRS,eid = IID,totalscore = TOTALSCORE_AVG),by = "eid")
data_sub$prs_std <- scale(data_sub$totalscore)[,1]
data_sub <- na.omit(data_sub)
colnames(data_sub) <- c("eid","tg_res","ptv_gene","prs_total","prs_std")
c1 <- sum(data_sub$ptv_gene>0)
model1 <- rlm(tg_res~ptv_gene,data = data_sub)
model2 <- rlm(tg_res~prs_std,data = data_sub)
model3 <- rlm(tg_res~ptv_gene+prs_std,data = data_sub)
model4 <- rlm(tg_res~ptv_gene+prs_std+ptv_gene:prs_std,data = data_sub)
lof_prs_results_tg[i,1] <- genes[i]
lof_prs_results_tg[i,2] <- BIC(model1)
lof_prs_results_tg[i,3] <- BIC(model2)
lof_prs_results_tg[i,4] <- BIC(model3)
lof_prs_results_tg[i,5] <- BIC(model4)
lof_prs_results_tg[i,6] <- summary(model4)$coeff[2,1]
lof_prs_results_tg[i,7] <- summary(model4)$coeff[3,1]
lof_prs_results_tg[i,8] <- summary(model4)$coeff[4,1]
lof_prs_results_tg[i,9] <- summary(model4)$coeff[2,3]
lof_prs_results_tg[i,10] <- summary(model4)$coeff[3,3]
lof_prs_results_tg[i,11] <- summary(model4)$coeff[4,3]
lof_prs_results_tg[i,12] <- c1
if(ceiling(i/10)==(i/10)){print(i)}
}
colnames(lof_prs_results_tg) <- c("PTV_gene","model1_BIC","model2_BIC","model3_BIC","model4_BIC",
"PTV_E","PRS_E","PTV_PRS_E","PTV_t","PRS_t","PTV_PRS_t","N_carrier_PTV")
# merge results for FDR-correction
lof_prs_results_hdl$trait <- "HDL"
lof_prs_results_ldl$trait <- "LDL"
lof_prs_results_tc$trait <- "TC"
lof_prs_results_tg$trait <- "TG"
lof_prs_results_hdl$PTV_P <- pnorm(abs(lof_prs_results_hdl$PTV_t),lower.tail=F)*2
lof_prs_results_hdl$PRS_P <- pnorm(abs(lof_prs_results_hdl$PRS_t),lower.tail=F)*2
lof_prs_results_hdl$PTV_PRS_P <- pnorm(abs(lof_prs_results_hdl$PTV_PRS_t),lower.tail=F)*2
lof_prs_results_ldl$PTV_P <- pnorm(abs(lof_prs_results_ldl$PTV_t),lower.tail=F)*2
lof_prs_results_ldl$PRS_P <- pnorm(abs(lof_prs_results_ldl$PRS_t),lower.tail=F)*2
lof_prs_results_ldl$PTV_PRS_P <- pnorm(abs(lof_prs_results_ldl$PTV_PRS_t),lower.tail=F)*2
lof_prs_results_tc$PTV_P <- pnorm(abs(lof_prs_results_tc$PTV_t),lower.tail=F)*2
lof_prs_results_tc$PRS_P <- pnorm(abs(lof_prs_results_tc$PRS_t),lower.tail=F)*2
lof_prs_results_tc$PTV_PRS_P <- pnorm(abs(lof_prs_results_tc$PTV_PRS_t),lower.tail=F)*2
lof_prs_results_tg$PTV_P <- pnorm(abs(lof_prs_results_tg$PTV_t),lower.tail=F)*2
lof_prs_results_tg$PRS_P <- pnorm(abs(lof_prs_results_tg$PRS_t),lower.tail=F)*2
lof_prs_results_tg$PTV_PRS_P <- pnorm(abs(lof_prs_results_tg$PTV_PRS_t),lower.tail=F)*2
all_results <- rbind(lof_prs_results_hdl,lof_prs_results_ldl,lof_prs_results_tc,lof_prs_results_tg)
all_results_long <- as.data.frame(matrix(NA,nrow(all_results)*3,4))
colnames(all_results_long) <- c("trait","PTV","name","P")
for(i in 1:nrow(all_results)){
all_results_long[(3*i-2):(3*i),1] <- all_results$trait[i]
all_results_long[(3*i-2):(3*i),2] <- all_results$PTV_gene[i]
all_results_long[(3*i-2):(3*i),3] <- c("PTV","PRS","PTV_PRS")
all_results_long[(3*i-2):(3*i),4] <- as.numeric(all_results[i,14:16])
if(ceiling(i/100)==(i/100)){print(i)}
}
all_results_long$fdr_q <- p.adjust(all_results_long$P,"fdr")
temp <- as.data.frame(matrix(NA,nrow(all_results),5))
colnames(temp) <- c("trait","PTV_gene","PTV_fdr_q","PRS_fdr_q","PTV_PRS_fdr_q")
for(i in 1:nrow(temp)){
temp[i,1] <- all_results_long$trait[(3*i-2)]
temp[i,2] <- all_results_long$PTV[(3*i-2)]
temp[i,3] <- all_results_long$fdr_q[(3*i-2)]
temp[i,4] <- all_results_long$fdr_q[(3*i-1)]
temp[i,5] <- all_results_long$fdr_q[(3*i)]
if(ceiling(i/100)==(i/100)){print(i)}
}
all_results_merge <- merge(all_results,temp,by = c("trait","PTV_gene"))