|
114
|
1 options("warn"=-1)
|
|
|
2
|
|
|
3 #from http://selection.med.yale.edu/baseline/Archive/Baseline%20Version%201.3/Baseline_Functions_Version1.3.r
|
|
|
4 # Compute p-value of two distributions
|
|
|
5 compareTwoDistsFaster <-function(sigma_S=seq(-20,20,length.out=4001), N=10000, dens1=runif(4001,0,1), dens2=runif(4001,0,1)){
|
|
|
6 #print(c(length(dens1),length(dens2)))
|
|
|
7 if(length(dens1)>1 & length(dens2)>1 ){
|
|
|
8 dens1<-dens1/sum(dens1)
|
|
|
9 dens2<-dens2/sum(dens2)
|
|
|
10 cum2 <- cumsum(dens2)-dens2/2
|
|
|
11 tmp<- sum(sapply(1:length(dens1),function(i)return(dens1[i]*cum2[i])))
|
|
|
12 #print(tmp)
|
|
|
13 if(tmp>0.5)tmp<-tmp-1
|
|
|
14 return( tmp )
|
|
|
15 }
|
|
|
16 else {
|
|
|
17 return(NA)
|
|
|
18 }
|
|
|
19 #return (sum(sapply(1:N,function(i)(sample(sigma_S,1,prob=dens1)>sample(sigma_S,1,prob=dens2))))/N)
|
|
|
20 }
|
|
|
21
|
|
|
22
|
|
|
23 require("grid")
|
|
|
24 arg <- commandArgs(TRUE)
|
|
|
25 #arg <- c("300143","4","5")
|
|
|
26 arg[!arg=="clonal"]
|
|
|
27 input <- arg[1]
|
|
|
28 output <- arg[2]
|
|
|
29 rowIDs <- as.numeric( sapply(arg[3:(max(3,length(arg)))],function(x){ gsub("chkbx","",x) } ) )
|
|
|
30
|
|
|
31 numbSeqs = length(rowIDs)
|
|
|
32
|
|
|
33 if ( is.na(rowIDs[1]) | numbSeqs>10 ) {
|
|
|
34 stop( paste("Error: Please select between one and 10 seqeunces to compare.") )
|
|
|
35 }
|
|
|
36
|
|
|
37 #load( paste("output/",sessionID,".RData",sep="") )
|
|
|
38 load( input )
|
|
|
39 #input
|
|
|
40
|
|
|
41 xMarks = seq(-20,20,length.out=4001)
|
|
|
42
|
|
|
43 plot_grid_s<-function(pdf1,pdf2,Sample=100,cex=1,xlim=NULL,xMarks = seq(-20,20,length.out=4001)){
|
|
|
44 yMax = max(c(abs(as.numeric(unlist(listPDFs[pdf1]))),abs(as.numeric(unlist(listPDFs[pdf2]))),0),na.rm=T) * 1.1
|
|
|
45
|
|
|
46 if(length(xlim==2)){
|
|
|
47 xMin=xlim[1]
|
|
|
48 xMax=xlim[2]
|
|
|
49 } else {
|
|
|
50 xMin_CDR = xMarks[listPDFs[pdf1][[1]][["CDR"]]>0.001][1]
|
|
|
51 xMin_FWR = xMarks[listPDFs[pdf1][[1]][["FWR"]]>0.001][1]
|
|
|
52 xMax_CDR = xMarks[listPDFs[pdf1][[1]][["CDR"]]>0.001][length(xMarks[listPDFs[pdf1][[1]][["CDR"]]>0.001])]
|
|
|
53 xMax_FWR = xMarks[listPDFs[pdf1][[1]][["FWR"]]>0.001][length(xMarks[listPDFs[pdf1][[1]][["FWR"]]>0.001])]
|
|
|
54
|
|
|
55 xMin_CDR2 = xMarks[listPDFs[pdf2][[1]][["CDR"]]>0.001][1]
|
|
|
56 xMin_FWR2 = xMarks[listPDFs[pdf2][[1]][["FWR"]]>0.001][1]
|
|
|
57 xMax_CDR2 = xMarks[listPDFs[pdf2][[1]][["CDR"]]>0.001][length(xMarks[listPDFs[pdf2][[1]][["CDR"]]>0.001])]
|
|
|
58 xMax_FWR2 = xMarks[listPDFs[pdf2][[1]][["FWR"]]>0.001][length(xMarks[listPDFs[pdf2][[1]][["FWR"]]>0.001])]
|
|
|
59
|
|
|
60 xMin=min(c(xMin_CDR,xMin_FWR,xMin_CDR2,xMin_FWR2,0),na.rm=TRUE)
|
|
|
61 xMax=max(c(xMax_CDR,xMax_FWR,xMax_CDR2,xMax_FWR2,0),na.rm=TRUE)
|
|
|
62 }
|
|
|
63
|
|
|
64 sigma<-approx(xMarks,xout=seq(xMin,xMax,length.out=Sample))$x
|
|
|
65 grid.rect(gp = gpar(col=gray(0.6),fill="white",cex=cex))
|
|
|
66 x <- sigma
|
|
|
67 pushViewport(viewport(x=0.175,y=0.175,width=0.825,height=0.825,just=c("left","bottom"),default.units="npc"))
|
|
|
68 #pushViewport(plotViewport(c(1.8, 1.8, 0.25, 0.25)*cex))
|
|
|
69 pushViewport(dataViewport(x, c(yMax,-yMax),gp = gpar(cex=cex),extension=c(0.05)))
|
|
|
70 grid.polygon(c(0,0,1,1),c(0,0.5,0.5,0),gp=gpar(col=grey(0.95),fill=grey(0.95)),default.units="npc")
|
|
|
71 grid.polygon(c(0,0,1,1),c(1,0.5,0.5,1),gp=gpar(col=grey(0.9),fill=grey(0.9)),default.units="npc")
|
|
|
72 grid.rect()
|
|
|
73 grid.xaxis(gp = gpar(cex=cex/1.1))
|
|
|
74 yticks = pretty(c(-yMax,yMax),8)
|
|
|
75 yticks = yticks[yticks>(-yMax) & yticks<(yMax)]
|
|
|
76 grid.yaxis(at=yticks,label=abs(yticks),gp = gpar(cex=cex/1.1))
|
|
|
77 if(length(listPDFs[pdf1][[1]][["CDR"]])>1){
|
|
|
78 ycdr<-approx(xMarks,listPDFs[pdf1][[1]][["CDR"]],xout=seq(xMin,xMax,length.out=Sample),yleft=0,yright=0)$y
|
|
|
79 grid.lines(unit(x,"native"), unit(ycdr,"native"),gp=gpar(col=2,lwd=2))
|
|
|
80 }
|
|
|
81 if(length(listPDFs[pdf1][[1]][["FWR"]])>1){
|
|
|
82 yfwr<-approx(xMarks,listPDFs[pdf1][[1]][["FWR"]],xout=seq(xMin,xMax,length.out=Sample),yleft=0,yright=0)$y
|
|
|
83 grid.lines(unit(x,"native"), unit(-yfwr,"native"),gp=gpar(col=4,lwd=2))
|
|
|
84 }
|
|
|
85
|
|
|
86 if(length(listPDFs[pdf2][[1]][["CDR"]])>1){
|
|
|
87 ycdr2<-approx(xMarks,listPDFs[pdf2][[1]][["CDR"]],xout=seq(xMin,xMax,length.out=Sample),yleft=0,yright=0)$y
|
|
|
88 grid.lines(unit(x,"native"), unit(ycdr2,"native"),gp=gpar(col=2,lwd=2,lty=2))
|
|
|
89 }
|
|
|
90 if(length(listPDFs[pdf2][[1]][["FWR"]])>1){
|
|
|
91 yfwr2<-approx(xMarks,listPDFs[pdf2][[1]][["FWR"]],xout=seq(xMin,xMax,length.out=Sample),yleft=0,yright=0)$y
|
|
|
92 grid.lines(unit(x,"native"), unit(-yfwr2,"native"),gp=gpar(col=4,lwd=2,lty=2))
|
|
|
93 }
|
|
|
94
|
|
|
95 grid.lines(unit(c(0,1),"npc"), unit(c(0.5,0.5),"npc"),gp=gpar(col=1))
|
|
|
96 grid.lines(unit(c(0,0),"native"), unit(c(0,1),"npc"),gp=gpar(col=1,lwd=1,lty=3))
|
|
|
97
|
|
|
98 grid.text("Density", x = unit(-2.5, "lines"), rot = 90,gp = gpar(cex=cex))
|
|
|
99 grid.text( expression(paste("Selection Strength (", Sigma, ")", sep="")) , y = unit(-2.5, "lines"),gp = gpar(cex=cex))
|
|
|
100
|
|
|
101 if(pdf1==pdf2 & length(listPDFs[pdf2][[1]][["FWR"]])>1 & length(listPDFs[pdf2][[1]][["CDR"]])>1 ){
|
|
|
102 pCDRFWR = compareTwoDistsFaster(sigma_S=xMarks, N=10000, dens1=listPDFs[[pdf1]][["CDR"]], dens2=listPDFs[[pdf1]][["FWR"]])
|
|
|
103 pval = formatC(as.numeric(pCDRFWR),digits=3)
|
|
|
104 grid.text( substitute(expression(paste(P[CDR/FWR], "=", x, sep="")),list(x=pval))[[2]] , x = unit(0.02, "npc"),y = unit(0.98, "npc"),just=c("left", "top"),gp = gpar(cex=cex*1.2))
|
|
|
105 }
|
|
|
106 grid.text(paste("CDR"), x = unit(0.98, "npc"),y = unit(0.98, "npc"),just=c("right", "top"),gp = gpar(cex=cex*1.5))
|
|
|
107 grid.text(paste("FWR"), x = unit(0.98, "npc"),y = unit(0.02, "npc"),just=c("right", "bottom"),gp = gpar(cex=cex*1.5))
|
|
|
108 popViewport(2)
|
|
|
109 }
|
|
|
110 #plot_grid_s(1)
|
|
|
111
|
|
|
112
|
|
|
113 p2col<-function(p=0.01){
|
|
|
114 breaks=c(-.51,-0.1,-.05,-0.01,-0.005,0,0.005,0.01,0.05,0.1,0.51)
|
|
|
115 i<-findInterval(p,breaks)
|
|
|
116 cols = c( rgb(0.8,1,0.8), rgb(0.6,1,0.6), rgb(0.4,1,0.4), rgb(0.2,1,0.2) , rgb(0,1,0),
|
|
|
117 rgb(1,0,0), rgb(1,.2,.2), rgb(1,.4,.4), rgb(1,.6,.6) , rgb(1,.8,.8) )
|
|
|
118 return(cols[i])
|
|
|
119 }
|
|
|
120
|
|
|
121
|
|
|
122 plot_pvals<-function(pdf1,pdf2,cex=1,upper=TRUE){
|
|
|
123 if(upper){
|
|
|
124 pCDR1FWR2 = compareTwoDistsFaster(sigma_S=xMarks, N=10000, dens1=listPDFs[[pdf1]][["CDR"]], dens2=listPDFs[[pdf2]][["FWR"]])
|
|
|
125 pFWR1FWR2 = compareTwoDistsFaster(sigma_S=xMarks, N=10000, dens1=listPDFs[[pdf1]][["FWR"]], dens2=listPDFs[[pdf2]][["FWR"]])
|
|
|
126 pFWR1CDR2 = compareTwoDistsFaster(sigma_S=xMarks, N=10000, dens2=listPDFs[[pdf2]][["CDR"]], dens1=listPDFs[[pdf1]][["FWR"]])
|
|
|
127 pCDR1CDR2 = compareTwoDistsFaster(sigma_S=xMarks, N=10000, dens2=listPDFs[[pdf2]][["CDR"]], dens1=listPDFs[[pdf1]][["CDR"]])
|
|
|
128 grid.polygon(c(0.5,0.5,1,1),c(0,0.5,0.5,0),gp=gpar(col=p2col(pFWR1FWR2),fill=p2col(pFWR1FWR2)),default.units="npc")
|
|
|
129 grid.polygon(c(0.5,0.5,1,1),c(1,0.5,0.5,1),gp=gpar(col=p2col(pCDR1FWR2),fill=p2col(pCDR1FWR2)),default.units="npc")
|
|
|
130 grid.polygon(c(0.5,0.5,0,0),c(1,0.5,0.5,1),gp=gpar(col=p2col(pCDR1CDR2),fill=p2col(pCDR1CDR2)),default.units="npc")
|
|
|
131 grid.polygon(c(0.5,0.5,0,0),c(0,0.5,0.5,0),gp=gpar(col=p2col(pFWR1CDR2),fill=p2col(pFWR1CDR2)),default.units="npc")
|
|
|
132
|
|
|
133 grid.lines(c(0,1),0.5,gp=gpar(lty=2,col=gray(0.925)))
|
|
|
134 grid.lines(0.5,c(0,1),gp=gpar(lty=2,col=gray(0.925)))
|
|
|
135
|
|
|
136 grid.text(formatC(as.numeric(pFWR1FWR2),digits=3), x = unit(0.75, "npc"),y = unit(0.25, "npc"),just=c("center", "center"),gp = gpar(cex=cex))
|
|
|
137 grid.text(formatC(as.numeric(pCDR1FWR2),digits=3), x = unit(0.75, "npc"),y = unit(0.75, "npc"),just=c("center", "center"),gp = gpar(cex=cex))
|
|
|
138 grid.text(formatC(as.numeric(pCDR1CDR2),digits=3), x = unit(0.25, "npc"),y = unit(0.75, "npc"),just=c("center", "center"),gp = gpar(cex=cex))
|
|
|
139 grid.text(formatC(as.numeric(pFWR1CDR2),digits=3), x = unit(0.25, "npc"),y = unit(0.25, "npc"),just=c("center", "center"),gp = gpar(cex=cex))
|
|
|
140
|
|
|
141
|
|
|
142 # grid.text(paste("P = ",formatC(pCDRFWR,digits=3)), x = unit(0.5, "npc"),y = unit(0.98, "npc"),just=c("center", "top"),gp = gpar(cex=cex))
|
|
|
143 # grid.text(paste("P = ",formatC(pFWRFWR,digits=3)), x = unit(0.5, "npc"),y = unit(0.02, "npc"),just=c("center", "bottom"),gp = gpar(cex=cex))
|
|
|
144 }
|
|
|
145 else{
|
|
|
146 }
|
|
|
147 }
|
|
|
148
|
|
|
149
|
|
|
150 ##################################################################################
|
|
|
151 ################## The whole OCD's matrix ########################################
|
|
|
152 ##################################################################################
|
|
|
153
|
|
|
154 #pdf(width=4*numbSeqs+1/3,height=4*numbSeqs+1/3)
|
|
|
155 pdf( output ,width=4*numbSeqs+1/3,height=4*numbSeqs+1/3)
|
|
|
156
|
|
|
157 pushViewport(viewport(x=0.02,y=0.02,just = c("left", "bottom"),w =0.96,height=0.96,layout = grid.layout(numbSeqs+1,numbSeqs+1,widths=unit.c(unit(rep(1,numbSeqs),"null"),unit(4,"lines")),heights=unit.c(unit(4,"lines"),unit(rep(1,numbSeqs),"null")))))
|
|
|
158
|
|
|
159 for( seqOne in 1:numbSeqs+1){
|
|
|
160 pushViewport(viewport(layout.pos.col = seqOne-1, layout.pos.row = 1))
|
|
|
161 if(seqOne>2){
|
|
|
162 grid.polygon(c(0,0,0.5,0.5),c(0,0.5,0.5,0),gp=gpar(col=grey(0.5),fill=grey(0.9)),default.units="npc")
|
|
|
163 grid.polygon(c(1,1,0.5,0.5),c(0,0.5,0.5,0),gp=gpar(col=grey(0.5),fill=grey(0.95)),default.units="npc")
|
|
|
164 grid.polygon(c(0,0,1,1),c(1,0.5,0.5,1),gp=gpar(col=grey(0.5)),default.units="npc")
|
|
|
165
|
|
|
166 grid.text(y=.25,x=0.75,"FWR",gp = gpar(cex=1.5),just="center")
|
|
|
167 grid.text(y=.25,x=0.25,"CDR",gp = gpar(cex=1.5),just="center")
|
|
|
168 }
|
|
|
169 grid.rect(gp = gpar(col=grey(0.9)))
|
|
|
170 grid.text(y=.75,substr(paste(names(listPDFs)[rowIDs[seqOne-1]]),1,16),gp = gpar(cex=2),just="center")
|
|
|
171 popViewport(1)
|
|
|
172 }
|
|
|
173
|
|
|
174 for( seqOne in 1:numbSeqs+1){
|
|
|
175 pushViewport(viewport(layout.pos.row = seqOne, layout.pos.col = numbSeqs+1))
|
|
|
176 if(seqOne<=numbSeqs){
|
|
|
177 grid.polygon(c(0,0.5,0.5,0),c(0,0,0.5,0.5),gp=gpar(col=grey(0.5),fill=grey(0.95)),default.units="npc")
|
|
|
178 grid.polygon(c(0,0.5,0.5,0),c(1,1,0.5,0.5),gp=gpar(col=grey(0.5),fill=grey(0.9)),default.units="npc")
|
|
|
179 grid.polygon(c(1,0.5,0.5,1),c(0,0,1,1),gp=gpar(col=grey(0.5)),default.units="npc")
|
|
|
180 grid.text(x=.25,y=0.75,"CDR",gp = gpar(cex=1.5),just="center",rot=270)
|
|
|
181 grid.text(x=.25,y=0.25,"FWR",gp = gpar(cex=1.5),just="center",rot=270)
|
|
|
182 }
|
|
|
183 grid.rect(gp = gpar(col=grey(0.9)))
|
|
|
184 grid.text(x=0.75,substr(paste(names(listPDFs)[rowIDs[seqOne-1]]),1,16),gp = gpar(cex=2),rot=270,just="center")
|
|
|
185 popViewport(1)
|
|
|
186 }
|
|
|
187
|
|
|
188 for( seqOne in 1:numbSeqs+1){
|
|
|
189 for(seqTwo in 1:numbSeqs+1){
|
|
|
190 pushViewport(viewport(layout.pos.col = seqTwo-1, layout.pos.row = seqOne))
|
|
|
191 if(seqTwo>seqOne){
|
|
|
192 plot_pvals(rowIDs[seqOne-1],rowIDs[seqTwo-1],cex=2)
|
|
|
193 grid.rect()
|
|
|
194 }
|
|
|
195 popViewport(1)
|
|
|
196 }
|
|
|
197 }
|
|
|
198
|
|
|
199
|
|
|
200 xMin=0
|
|
|
201 xMax=0.01
|
|
|
202 for(pdf1 in rowIDs){
|
|
|
203 xMin_CDR = xMarks[listPDFs[pdf1][[1]][["CDR"]]>0.001][1]
|
|
|
204 xMin_FWR = xMarks[listPDFs[pdf1][[1]][["FWR"]]>0.001][1]
|
|
|
205 xMax_CDR = xMarks[listPDFs[pdf1][[1]][["CDR"]]>0.001][length(xMarks[listPDFs[pdf1][[1]][["CDR"]]>0.001])]
|
|
|
206 xMax_FWR = xMarks[listPDFs[pdf1][[1]][["FWR"]]>0.001][length(xMarks[listPDFs[pdf1][[1]][["FWR"]]>0.001])]
|
|
|
207 xMin=min(c(xMin_CDR,xMin_FWR,xMin),na.rm=TRUE)
|
|
|
208 xMax=max(c(xMax_CDR,xMax_FWR,xMax),na.rm=TRUE)
|
|
|
209 }
|
|
|
210
|
|
|
211
|
|
|
212
|
|
|
213 for(i in 1:numbSeqs+1){
|
|
|
214 for(j in (i-1):numbSeqs){
|
|
|
215 pushViewport(viewport(layout.pos.col = i-1, layout.pos.row = j+1))
|
|
|
216 grid.rect()
|
|
|
217 plot_grid_s(rowIDs[i-1],rowIDs[j],cex=1)
|
|
|
218 popViewport(1)
|
|
|
219 }
|
|
|
220 }
|
|
|
221
|
|
|
222 dev.off()
|
|
|
223
|
|
|
224 cat("Success", paste(rowIDs,collapse="_"),sep=":")
|
|
|
225
|
