QQ plot

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

QQ plot

anikaM
Hello,

I was plotting QQ plot (in attach) via:

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=45, cex=1)
    abline(-log10(0.10),1, col='black',lty=1)
    text(0.5, 1.6, "FDR=0.10", col = "gray60",srt=45, cex=1)
    abline(-log10(0.25),1, col='black',lty=1)
    text(0.5, 1.2, "FDR=0.25", col = "gray60",srt=45, 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
    }
  }
}

db=fread("/Users/ams/Desktop/GWAS_all.txt", header=T)
dd=db[db$P<1e-3,]
qqunif(dd$P)

what do I need to change in my code so that FDR lines look like on the
2nd attachment?

______________________________________________
[hidden email] 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.

first.png (200K) Download Attachment
second.png (286K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: QQ plot

Jim Lemon-4
Hi Ana Marija,
It would help if we had the data, but I suspect that the data that we
don't have is not uniformly distributed:

library(plotrix)
ams_norm<-rnorm(100)
png("ams_norm.png")
qqunif(rescale(ams_norm,c(0,1)))
dev.off()
ams_unif<-runif(100)
png("ams_unif.png")
qqunif(ams_unif)
dev.off()

Jim

On Thu, Aug 15, 2019 at 3:27 AM Ana Marija <[hidden email]> wrote:

>
> Hello,
>
> I was plotting QQ plot (in attach) via:
>
> 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=45, cex=1)
>     abline(-log10(0.10),1, col='black',lty=1)
>     text(0.5, 1.6, "FDR=0.10", col = "gray60",srt=45, cex=1)
>     abline(-log10(0.25),1, col='black',lty=1)
>     text(0.5, 1.2, "FDR=0.25", col = "gray60",srt=45, 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
>     }
>   }
> }
>
> db=fread("/Users/ams/Desktop/GWAS_all.txt", header=T)
> dd=db[db$P<1e-3,]
> qqunif(dd$P)
>
> what do I need to change in my code so that FDR lines look like on the
> 2nd attachment?
> ______________________________________________
> [hidden email] 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.

______________________________________________
[hidden email] 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.

ams_norm.png (3K) Download Attachment
ams_unif.png (3K) Download Attachment