-
Notifications
You must be signed in to change notification settings - Fork 0
/
run-meta-analyses-within-ancestries.R
146 lines (116 loc) · 5.52 KB
/
run-meta-analyses-within-ancestries.R
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
setwd("/home/l.yengo/loic/GIANT2020/trans-ethnic-analyses")
M <- 1385132
cc <- c(rep("character",3),rep("numeric",5))
fe_meta_analysis <- function(studies,outputname){
nstudies <- length(studies)
Bmat <- matrix(NA,nrow=M,ncol=nstudies)
SEmat <- matrix(NA,nrow=M,ncol=nstudies)
Nmat <- matrix(NA,nrow=M,ncol=nstudies)
Fmat <- matrix(NA,nrow=M,ncol=nstudies)
var_y <- rep(NA,nstudies)
for(i in 1:nstudies){
cat(paste0("\tReading ",studies[i],"..."))
tmp <- read.table(studies[i],header=TRUE,stringsAsFactors=FALSE,colClasses=cc)
Bmat[,i] <- tmp[,"beta"]
SEmat[,i] <- tmp[,"se"]
Nmat[,i] <- tmp[,"N"]
Fmat[,i] <- tmp[,"freq"]
var_y[i] <- median( tmp[,"N"] * 2 * tmp[,"freq"] * (1-tmp[,"freq"]) * tmp[,"se"] * tmp[,"se"], na.rm=TRUE)
cat(paste0("E[var(Y)] = ",round(var_y[i],4),".\n"))
}
cat("\tCorr of frequencies\n")
print(cor(Fmat,use="pairwise"))
cat("\tCorr of BETA's\n")
print(cor(Bmat,use="pairwise"))
slopeMat <- matrix(NA,nrow=nstudies,ncol=nstudies)
for(i in 1:nstudies){
for(j in 1:nstudies){
slopeMat[i,j] <- ifelse(i==j,1,coef(lm(Bmat[,i]~Bmat[,j]))[2])
}
}
cat("\tSlope of BETA's\n")
print(slopeMat)
wgtr <- 1/(SEmat^2)
wgts <- wgtr / rowSums(wgtr,na.rm=TRUE)
BETA <- rowSums(Bmat * wgts,na.rm=TRUE)
SE <- 1/sqrt(rowSums(wgtr,na.rm=TRUE))
Z <- BETA/SE
CHISQ <- Z*Z
P <- pchisq(q=CHISQ,df=1,lower.tail=FALSE)
N <- rowSums(Nmat,na.rm=TRUE)
FREQ <- rowSums(Nmat*Fmat,na.rm=TRUE)/N
## filters
f <- which(SE==0 | SE==Inf | abs(BETA)==Inf | N==0 | N==Inf | is.na(N) | is.na(BETA) | is.na(SE) | is.na(P) | is.na(FREQ) | FREQ==0 | FREQ==1)
if(length(f)>0){
print("Issue here!")
BETA[f] <- SE[f] <- P[f] <- N[f] <- FREQ[f] <- NA
}
tmp[,"freq"] <- FREQ
tmp[,"beta"] <- BETA
tmp[,"se"] <- SE
tmp[,"p"] <- P
tmp[,"N"] <- N
write.table(tmp,outputname,row.names=FALSE,quote=FALSE)
}
sourceDir <- "/gpfs/gpfs01/polaris/Q0286/loic/GIANT2020/HM3/"
eur <- paste0(sourceDir,c("23andMe_EUR_int_height_female.HM3_v2.txt",
"23andMe_EUR_int_height_male.HM3_v2.txt",
"HEIGHT_EUR_ALL_FINAL_allchr_HM3_v2.txt"))
eas <- paste0(sourceDir,c("23andMe_EAS_int_height_female.HM3_v2.txt",
"23andMe_EAS_int_height_male.HM3_v2.txt",
"HEIGHT_EAS_1KGP3_withOut_UKB_FINAL_allchr_HM3_v2.txt"))
sas <- paste0(sourceDir,c("23andMe_SAS_int_height_female.HM3_v2.txt",
"23andMe_SAS_int_height_male.HM3_v2.txt",
"HEIGHT_SAS_1KGP3_withOUTUKB_FINAL_allchr_HM3_v2.txt"))
his <- paste0(sourceDir,c("23andMe_HA_int_height_female.HM3_v2.txt",
"23andMe_HA_int_height_male.HM3_v2.txt",
"HEIGHT_HIS_1KGP3_FINAL_allchr_HM3_v2.txt"))
aa <- paste0(sourceDir,c("23andMe_AA_int_height_female.HM3_v2.txt",
"23andMe_AA_int_height_male.HM3_v2.txt",
"HEIGHT_AA_1KGP3_FINAL_allchr_HM3_v2.txt"))
eur2 <- paste0(sourceDir,c("23andMe_EUR_int_height_female.HM3_v2.txt",
"23andMe_EUR_int_height_male.HM3_v2.txt",
"HEIGHT_EUR_ALL_withoutUKB_FINAL_allchr_HM3_v2.txt",
"HEIGHT_UKB_v3EUR_maf0001_nosibs_LMM_inf_HM3_v2.txt"))
giant2020 <- paste0(sourceDir,c("HEIGHT_EUR_ALL_withoutUKB_FINAL_allchr_HM3_v2.txt",
"HEIGHT_UKB_v3EUR_maf0001_withsibs_LMM_inf_HM3_v2.txt"))
print("Analysis 1")
#fe_meta_analysis(eur,"meta-analyses/EUR.ma")
print("Analysis 2")
#fe_meta_analysis(eas,"meta-analyses/EAS.ma")
print("Analysis 3")
#fe_meta_analysis(sas,"meta-analyses/SAS.ma")
print("Analysis 4")
#fe_meta_analysis(his,"meta-analyses/HIS.ma")
print("Analysis 5")
#fe_meta_analysis(aa,"meta-analyses/AA.ma")
#print("Analysis 6 -- naive")
#fe_meta_analysis(paste0("meta-analyses/",dir("meta-analyses")),"meta-analyses/ALL.ma")
print("Analysis 7 -- no UKB sibs")
#fe_meta_analysis(eur2,"meta-analyses/EUR_noUKBsibs.ma")
print("Analysis 8 - nosibs + UKB")
#fe_meta_analysis(giant2020,"meta-analyses/GIANT2020_inclUKB_v2.ma")
### New analyses >23APR2020
eas <- c(paste0(sourceDir,c("23andMe_EAS_int_height_female.HM3_v2.txt",
"23andMe_EAS_int_height_male.HM3_v2.txt")),
"../download/formatted/GIANT2020noUKB_EAS.ma")
afr <- c(paste0(sourceDir,c("23andMe_AA_int_height_female.HM3_v2.txt",
"23andMe_AA_int_height_male.HM3_v2.txt")),
"../download/formatted/GIANT2020noUKB_AFR.ma")
sas <- c(paste0(sourceDir,c("23andMe_SAS_int_height_female.HM3_v2.txt",
"23andMe_SAS_int_height_male.HM3_v2.txt")),
"../download/formatted/GIANT2020noUKB_SAS.ma")
his <- c(paste0(sourceDir,c("23andMe_HA_int_height_female.HM3_v2.txt",
"23andMe_HA_int_height_male.HM3_v2.txt")),
"../download/formatted/GIANT2020noUKB_HIS.ma")
eur <- c(paste0(sourceDir,c("23andMe_EUR_int_height_female.HM3_v2.txt",
"23andMe_EUR_int_height_male.HM3_v2.txt",
"HEIGHT_UKB_v3EUR_maf0001_nosibs_v2_LMM_inf_HM3_v2.txt")),
"../download/formatted/GIANT2020noUKB_EUR.ma")
all <- c(eas,sas,afr,his,eur)
fe_meta_analysis(eas,"meta-analyses/EAS.ma")
fe_meta_analysis(sas,"meta-analyses/SAS.ma")
fe_meta_analysis(afr,"meta-analyses/AFR.ma")
fe_meta_analysis(his,"meta-analyses/HIS.ma")
#fe_meta_analysis(eur,"meta-analyses/EUR.ma")
#fe_meta_analysis(all,"meta-analyses/ALL.ma")