The data

For this example I will use a dataset that I will cut in two to pretent they are two different populations. This is just for illustration purposes.

setwd("C:/Users/marti/OneDrive/Documents/GWA")

load("data_meta.RData")
ls()
## [1] "G"     "map"   "pheno" "Z"
partition<-sample(1:2,nrow(G),replace = T)
pheno_1<-pheno[partition==1,]
pheno_2<-pheno[partition==2,]
G1<-G[rownames(pheno_1),rownames(pheno_1)]
G2<-G[rownames(pheno_2),rownames(pheno_2)]
Z1<-Z[rownames(pheno_1),]
Z2<-Z[rownames(pheno_2),]

Fit two gblup

Notice that the models can be completely different

blup1<-gblup(rsp = "wt_birth",design = c(y~sex,~slgdt_cd),data = pheno_1,G=G1)
blup2<-gblup(rsp = "wt_birth",design = c(y~sex),data=pheno_2,G=G2)

print(blup1)
## gblup analysis of trait: wt_birth 
## 
## fixed effects equation:
## y ~ sex
## 
## random effects equation:
## ~G + slgdt_cd + In
## 
## log-likelihood: 283.5474 converged in: 7 iterations 
## 
## estimated fixed effects:
##               Estimate   StdError   test.st    p.value
## (Intercept) 1.49561770 0.03018116 49.554683 0.00000000
## sexM        0.06122943 0.02969517  2.061932 0.03921417
## 
## estimated variance components:
##            Estimate    StdError  prop.var
## G        0.01985394 0.008298874 0.1776490
## slgdt_cd 0.01539242 0.005875582 0.1377282
## In       0.07651301 0.007589883 0.6846228
print(blup2)
## gblup analysis of trait: wt_birth 
## 
## fixed effects equation:
## y ~ sex
## 
## random effects equation:
## ~G + In
## 
## log-likelihood: 346.3695 converged in: 7 iterations 
## 
## estimated fixed effects:
##               Estimate   StdError   test.st   p.value
## (Intercept) 1.51883656 0.02010573 75.542463 0.0000000
## sexM        0.04248117 0.02716534  1.563801 0.1178644
## 
## estimated variance components:
##      Estimate    StdError  prop.var
## G  0.01967956 0.007228473 0.2049805
## In 0.07632744 0.006641037 0.7950195

##perform sepparate GWA

gw1<-gwas(gblup = blup1,x = t(Z1))
head(gw1)
##                     ghat     var_ghat
## MARC0044150 -0.002791660 1.976481e-05
## ASGA0000014 -0.001490496 2.629636e-05
## H3GA0000026  0.003919056 1.980873e-05
## ASGA0000021  0.002791660 1.976481e-05
## ALGA0000009  0.003919056 1.980873e-05
## ALGA0000014  0.003919056 1.980873e-05
p1<-getpvalue(gw1,log.p = F)
manhattan_plot(p1,map=map)

gw2<-gwas(gblup = blup2,x = t(Z2))
head(gw2)
##                     ghat     var_ghat
## MARC0044150  0.005875252 2.129344e-05
## ASGA0000014 -0.001150058 3.185832e-05
## H3GA0000026 -0.006067477 2.131129e-05
## ASGA0000021 -0.005875252 2.129344e-05
## ALGA0000009 -0.006067477 2.131129e-05
## ALGA0000014 -0.006067477 2.131129e-05
p2<-getpvalue(gw2,log.p = F)
manhattan_plot(p2,map=map)

perform the meta-analysis GWA

For this it’s necesary to have the varian

w1<-1/((blup1$sigma["G"]^2)/gw1[,2])
w2<-1/((blup2$sigma["G"]^2)/gw2[,2])

z1<-gw1[,1]/sqrt(gw1[,2])
names(z1)<-rownames(gw1)
z2<-gw2[,1]/sqrt(gw2[,2])
names(z2)<-rownames(gw2)

sum(names(z1)!=names(z2))
## [1] 0
#stop if this is not zero
zm<-(z1*sqrt(w1)+z2*sqrt(w2))/sqrt(w1+w2)
pm<-getpvalue(gwas = zm,is.z = T,log.p = F)

manhattan_plot(pvalues = pm,map=map)

plot(z1,zm)

plot(z2,zm)