1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
|
# read filename from the command line
args <- commandArgs(trailingOnly=TRUE)
filename <- args[1]
outfile <- args[2]
pdf(file=outfile,height=10,width=7)
# read file
LOD <- as.matrix(read.table(filename,header=FALSE))
N <- ncol(LOD)
L <- nrow(LOD)
LOD <- matrix(LOD,ncol=N)
# sort LOD by last column
O <- order(LOD[,N])
LOD <- matrix(LOD[O,],ncol=N)
# create PP table
PP <- (10**LOD)/(1+10**LOD)
# split last column out of tables
aveLOD <- as.vector(LOD[,N])
avePP <- as.vector( PP[,N])
LOD <- matrix(LOD[,1:N-1],ncol=N-1)
PP <- matrix(PP[,1:N-1],ncol=N-1)
# set up the two plotting surfaces
par(mfrow=c(2,1))
#-------------- Plot 1 -------------------
plot(avePP,xlab="Split",ylab="PP",type="n",ylim=c(0,1),xaxt="n")
axis(side=1,at=1:L,1:L)
# plot results from another analysis for comparison
# lines(avePP_other,col=hsv(0.63,0.4,1),lwd=2)
# Plot the bars representing the range for each split
minPP <- apply(PP, 1, min, na.rm=TRUE)
maxPP <- apply(PP, 1, max, na.rm=TRUE)
for(i in 1:L) {
lines(c(i,i),c(minPP[i],maxPP[i]),lwd=2)
}
# Plot the estimate for each split, and connect the dots
lines(avePP,col=hsv(1,1,1),lwd=2)
#-------------- Plot 2 -------------------
plot(aveLOD,xlab="Split",ylab="LOD10", type="n",xaxt="n")
axis(side=1,at=1:L,1:L)
# plot results from another analysis for comparison
# lines(LOD,col=hsv(0.63,0.4,1),lwd=2)
# Plot the bars representing the range for each split
minLOD <- apply(LOD, 1, min, na.rm=TRUE)
maxLOD <- apply(LOD, 1, max, na.rm=TRUE)
for(i in 1:L) {
lines(c(i,i),c(minLOD[i],maxLOD[i]),lwd=2)
}
# Plot the estimate for each split, and connect the dots
lines(aveLOD,col=hsv(1,1,1),lwd=2)
|
