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
-------------- next part --------------
A non-text attachment was scrubbed...
Name: qqpl.png
Type: image/png
Size: 89152 bytes
Desc: not available
URL:
<https://stat.ethz.ch/pipermail/r-help/attachments/20200228/2d4ffabf/attachment.png>
Hi Ana, Look carefully at that red line. It goes through (0,0) and scoots off the plot at (2.5,2.5). As you have specified that intercept and slope in your code, poor abline is doing the best it can. Do not punish it for doing what you request. Jim On Sat, Feb 29, 2020 at 6:10 AM Ana Marija <sokovic.anamarija at gmail.com> wrote:> > 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 at 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.
Hi Jim, I have in my code: abline(0,1,col='red') can you please tell me how to change my code to have it indeed running from 0 to 1? Thanks Ana On Fri, Feb 28, 2020 at 3:34 PM Jim Lemon <drjimlemon at gmail.com> wrote:> > Hi Ana, > Look carefully at that red line. It goes through (0,0) and scoots off > the plot at (2.5,2.5). As you have specified that intercept and slope > in your code, poor abline is doing the best it can. Do not punish it > for doing what you request. > > Jim > > On Sat, Feb 29, 2020 at 6:10 AM Ana Marija <sokovic.anamarija at gmail.com> wrote: > > > > 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 at 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.