NULK

Data from Science 2017, 358:365-8 Recent natural selection causes adaptive evolution of an avian polygenic trait

source("manhattan.R")
plink='/Users/gc5k/bin/plink_mac/plink'

assoc=paste(plink, " --linear --bfile NLUK_auto --pheno NLUK_auto.eigenvec --out NLUK_auto --allow-no-sex --autosome-num 34")
system(assoc)

layout(matrix(c(1,2,2,1,3,3),2,3, byrow=T))
evec=read.table("NLUK_auto.eigenvec", as.is = T)
evec$geo=ifelse(evec$V1 == "Wytham_UK", "Wytham", ifelse(evec$V1 == "Oosterhout", "Oosterhot", "Veluwe forests"))
evec$col=ifelse(evec$V1 == "Wytham_UK", "grey", ifelse(evec$V1=="Oosterhout", "yellow", "green"))
plot(evec$V3, evec$V4, xlab="PC 1", ylab="PC 2", col=evec$col, pch=16, bty='l')
legend("topleft", legend = unique(evec$geo), col=unique(evec$col), bty='n', pch=16)

dat=read.table("NLUK_auto.assoc.linear", as.is = T, header = T)

manhattan(dat, limitchromosomes = 1:max(dat$CHR), pch=16, cex=0.5, bty='l')

gc=qchisq(median(dat$P, na.rm = T), 1, lower.tail = F)/qchisq(0.5, 1, lower.tail = F) #LambdaGC

dat1=dat
dat1$P=pchisq(dat$STAT^2/gc, 1, lower.tail = F)
manhattan(dat1, limitchromosomes = 1:max(dat1$CHR), pch=16, cex=0.5, bty='l')