Mercurial > repos > dereeper > mlmm

comparison source_library/plot_MLMM.Ago.r @ 1:380b364980f9 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 475f4d7d8442a0d75e103af326ae5881c4d2a4ac
author dereeper
date Mon, 16 Apr 2018 08:50:05 -0400
parents
children
comparison
equal deleted inserted replaced
0:6b7107812931 1:380b364980f9
1 ## fonction ###########################################################################
2 plot_MLMM.Ago<-function(x) {
3 output1 <-x$out[order(x$out$Pos),]
4 output_ok<-output1[order(output1$Chr),]
5
6 maxpos<-c(0,cumsum(as.numeric(aggregate(output_ok$Pos,list(output_ok$Chr),max)$x+max(cumsum(as.numeric(aggregate(output_ok$Pos,list(output_ok$Chr),max)$x)))/200)))
7 plot_col<-rep(c('gray10','gray60'),ceiling(max(unique(output_ok$Chr))/2))
8 # plot_col<-c('blue','darkgreen','red','cyan','purple')
9 size<-aggregate(output_ok$Pos,list(output_ok$Chr),length)$x
10 a<-rep(maxpos[1],size[1])
11 b<-rep(plot_col[1],size[1])
12 for (i in 2:max(unique(output_ok$Chr))){
13 a<-c(a,rep(maxpos[i],size[i]))
14 b<-c(b,rep(plot_col[i],size[i]))}
15
16 output_ok$xpos<-output_ok$Pos+a
17 output_ok$col<-b
18 output_ok$col[output_ok$mk=='qtl']<-'cyan'
19 output_ok$col[output_ok$SNP %in% x$cof]<-'red'
20
21 d<-(aggregate(output_ok$xpos,list(output_ok$Chr),min)$x+aggregate(output_ok$xpos,list(output_ok$Chr),max)$x)/2
22
23 plot(output_ok$xpos,-log10(output_ok$pval),col=output_ok$col,pch=20,ylab='-log10(pval)',xaxt='n',xlab='chromosome')
24 axis(1,tick=FALSE,at=d,labels=c(1:max(unique(output_ok$Chr))))
25 xline(output_ok$xpos[output_ok$mk=='qtl'], col='cyan', lwd=0.1)
26 abline(h=x$bonf_thresh,lty=3,col='black')
27
28
29 if (length(output_ok$pval[-log10(output_ok$pval) > x$bonf_thresh]) > 0) {
30 text(output_ok$xpos[-log10(output_ok$pval) > x$bonf_thresh], -log10(output_ok$pval[-log10(output_ok$pval) > x$bonf_thresh]), output_ok$SNP[-log10(output_ok$pval) > x$bonf_thresh], pos=3, cex=0.7)
31 legend("topright", lty=3, paste("bonf thresh :", x$bonf_thresh ,sep=" "))
32 } else {
33 legend("topright", lty=3, paste("bonf thresh :", x$bonf_thresh ,sep=" "))
34 }
35
36 }