Mercurial > repos > davidvanzessen > shm_csr

comparison baseline/comparePDFs.r @ 0:c33d93683a09 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Thu, 13 Oct 2016 10:52:24 -0400
parents
children a24f8c93583a
comparison
equal deleted inserted replaced
-1:000000000000 0:c33d93683a09
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