panp.R
In panp: Presence-Absence Calls from Negative Strand Matching Probesets

### R code from vignette source 'panp.Rnw'

###################################################
### code chunk number 1: LibraryPreload
###################################################
library(gcrma)
library(panp)


###################################################
### code chunk number 2: load.ExpressionSet
###################################################
data(gcrma.ExpressionSet)


###################################################
### code chunk number 3: usage.info
###################################################
pa.calls()


###################################################
### code chunk number 4: run.pa.calls
###################################################
PA <- pa.calls(gcrma.ExpressionSet)


###################################################
### code chunk number 5: extract.PAcalls
###################################################
PAcalls <- PA$Pcalls
Pvalues <- PA$Pvals

write.table(PAcalls, file="PAcalls_gcrma.csv", sep=",", col.names=NA)
write.table(Pvalues, file="Pvalues_gcrma.csv", sep=",", col.names=NA)


###################################################
### code chunk number 6: show.some.output
###################################################
head(PAcalls[,1])

head(Pvalues[,1])


###################################################
### code chunk number 7: extract.names
###################################################
P_list_1 <-rownames(PAcalls)[PAcalls[,1]=="P"]
M_list_1 <-rownames(PAcalls)[PAcalls[,1]=="M"]
A_list_1 <-rownames(PAcalls)[PAcalls[,1]=="A"]


###################################################
### code chunk number 8: panp.Rnw:136-188
###################################################
data(NSMPnames.hgu133a)	# read in the NSMP names list
# create x and y ranges for NSMP intensities
NSMP_x <- sort(exprs(gcrma.ExpressionSet)[NSMPnames.hgu133a,1], decreasing=TRUE)
NSMP_y <- seq(0,1,1/(length(NSMP_x)-1))

# Plot expression densities of all probesets, then just NSMPs
plot(density(exprs(gcrma.ExpressionSet)[,1]),
	col="blue",
	xlim = c(1,13),
	ylim = c(0,.8),
	main = "Expression density: NSMPs vs. all, and NSMP survivor curve",
	xlab = "Log2(Intensity)",
	ylab = "Probability density")
lines(density(NSMP_x), col=6)

# interpolate over the NSMP exprs to draw survivor function
interp <- approxfun(NSMP_x, NSMP_y, yleft=1, yright=0) 
x = NSMP_x
curve(interp(x),add=TRUE, lwd=2)	 #add it to the plot

# reverse interpolate to get intensity values at p-value cutoffs
revInterp <- approxfun(NSMP_y, NSMP_x, yleft=1, yright=0)
rev01=revInterp(0.01)
rev02=revInterp(0.02)
# Pinpoint the x-y locations
# points(rev01, .01, pch=21, cex=2, lwd=2,col=1)
# points(rev02, .02, pch=21, cex=2, lwd=2,col=1)

# Draw horiz. lines & labels for both Pval cutoffs: tightCutoff, looseCutoff:
abline(h = 0.01, col = 1, lty = 2)
abline(h = 0.02, col = 1, lty = 2)
text(2.4, 0.01, pos=3, offset=0.1, cex=.7, as.character(0.01))
text(1.7, 0.02, pos=3, offset=0.1, cex=.7, as.character(0.02))

# vertical lines & labels for interpolated intensities at cutoffs
revTight=revInterp(0.01)
revLoose=revInterp(0.02)
abline(v = revTight, col = 1, lty = 2)
abline(v = revLoose, col = 1, lty = 2)
text(revLoose, .35, pos=2, offset=0.1, cex=.7, format.pval(revLoose,digits=3))
text(revTight, .30, pos=2, offset=0.1, cex=.7, format.pval(revTight,digits=3))
text(revLoose, .35, pos=4, cex=.8, "Log(intensity)")
text(revTight, .30, pos=4, cex=.8, "Log(intensity)")

lines(density(exprs(gcrma.ExpressionSet)[,1][PAcalls[,1]=="P"], bw=.1, n=512),col=2, lty=2, lwd=1)
lines(density(exprs(gcrma.ExpressionSet)[,1][PAcalls[,1]=="A"], bw=.1, n=512),col=3, lty=2, lwd=1)

legend(7.5,.8, c("NSMP exprs, survivor fcn.","NSMP exprs, density","All probesets, density","'Present' probesets, density","'Absent' probesets, density"),
col = c(1,6,4,2,3), lty=c(1,1,1,2,2), lwd=c(2,1,1,1,1), cex=.75,
 	text.col= "darkgreen",
 	bg='gray90')