Mercurial > repos > davidvanzessen > mutation_analysis
comparison baseline/comparePDFs.r @ 0:8a5a2abbb870 draft default tip
Uploaded
author | davidvanzessen |
---|---|
date | Mon, 29 Aug 2016 05:36:10 -0400 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:8a5a2abbb870 |
---|---|
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 |