Hello,
I made the plot in attach with this function:
qqunif = function(p, BH=T, MAIN = " ", SUB=" ")
{
nn = length(p)
xx = -log10((1:nn)/(nn+1))
plot( xx, -sort(log10(p)),
main = MAIN, sub= SUB, cex.sub=1.3,
xlab=expression(Expected~~-log[10](italic(p))),
ylab=expression(Observed~~-log[10](italic(p))),
cex.lab=1.0,mgp=c(2,1,0))
abline(0,1,col='red')
if(BH) ## BH = include Benjamini Hochberg FDR
{
abline(-log10(0.05),1, col='black',lty=1)
text(0.5,1.9 , "FDR=0.05", col = "gray60",srt=20, cex=1)
abline(-log10(0.10),1, col='black',lty=1)
text(0.5, 1.6, "FDR=0.10", col = "gray60",srt=20, cex=1)
abline(-log10(0.25),1, col='black',lty=1)
text(0.5, 1.2, "FDR=0.25", col = "gray60",srt=20, cex=1)
#legend('topleft', c("FDR = 0.05","FDR = 0.10","FDR = 0.25"),
#col=c('black','black','black'),lty=c(1,1,1), cex=0.8)
if (BF)
{
abline(h=-log10(0.05/nn), col='black') ## bonferroni
}
}
}
biob272=read.table("/Users/ams/Desktop/biobank272LD.txt")
qqunif(biob272$V2)
> head(biob272)
V1 V2
1 rs2089177 0.581204
2 rs4360974 0.418456
3 rs6502526 0.416670
4 rs8069906 0.568030
5 rs9895995 0.266746
6 rs9905280 0.510032
But the red, abline doesn't look like it is 1:1 line.
Can you please advise?
Thanks
Ana
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
