HELP.R
In HELP: Tools for HELP data analysis

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

###################################################
### code chunk number 1: HELP.Rnw:45-46
###################################################
library(HELP)


###################################################
### code chunk number 2: HELP.Rnw:48-49
###################################################
data.path <- system.file("pipeline.data", package="HELP")


###################################################
### code chunk number 3: HELP.Rnw:81-95
###################################################
msp1.files <- dir(data.path, pattern="532[._]pair", full.names=TRUE)
hpa2.files <- dir(data.path, pattern="635[._]pair", full.names=TRUE)
N <- length(msp1.files)
for (i in 1:N) {
	if (i == 1) {
		pairs <- readPairs(msp1.files[i],hpa2.files[i])
	}
	else {
		pairs <- readPairs(msp1.files[i],hpa2.files[i],pairs)
	}
}
class(pairs)
dim(pairs)
sampleNames(pairs)


###################################################
### code chunk number 4: HELP.Rnw:105-111
###################################################
chips <- sub("[_.]*532[_.]*pair.*","",basename(msp1.files))
chips ## CHIP_IDs
samplekey <- file.path(data.path, "sample.key.txt")
chips <- readSampleKey(file=samplekey,chips=chips)
chips ## SAMPLEs (from supplied key)
sampleNames(pairs) <- chips


###################################################
### code chunk number 5: HELP.Rnw:121-126
###################################################
ndf.file <- file.path(data.path, "HELP.ndf.txt")
ngd.file <- file.path(data.path, "HELP.ngd.txt")
pairs <- readDesign(ndf.file, ngd.file, pairs)
pData(featureData(pairs))[1:10,c("CHR","START","STOP")]
getFeatures(pairs,"SEQUENCE")[1]


###################################################
### code chunk number 6: HELP.Rnw:137-139
###################################################
calcTm("ATCTAGGAAGATTTAGAGAGGCAATGTGTCATTTAGCATCTAATTTTACC")
calcTm(getFeatures(pairs,"SEQUENCE")[1:4])


###################################################
### code chunk number 7: HELP.Rnw:145-147
###################################################
calcGC("ATCTAGGAAGATTTAGAGAGGCAATGTGTCATTTAGCATCTAATTTTACC")
calcGC(getFeatures(pairs,"SEQUENCE")[1:4])


###################################################
### code chunk number 8: HELP.Rnw:163-165
###################################################
getSamples(pairs)[1:4,]
calcPrototype(pairs,center=FALSE)[1:4]


###################################################
### code chunk number 9: fig1
###################################################
plotChip(pairs, sample="Brain2")


###################################################
### code chunk number 10: fig2
###################################################
x <- rep(1:100,100)
y <- rep(1:100,each=100)
z <- x*(1001:11000/1000)
z <- z-mean(z)
z <- z*(sample(1:10000/10000)+1)
plotChip(x,y,z,main="Curved gradient",xlab="x",ylab="y")


###################################################
### code chunk number 11: HELP.Rnw:222-224
###################################################
exprs(pairs) <- log2(exprs(pairs))
exprs2(pairs) <- log2(exprs2(pairs))


###################################################
### code chunk number 12: fig3
###################################################
plotFeature(pairs[,"Brain2"],cex=0.5)


###################################################
### code chunk number 13: HELP.Rnw:247-255
###################################################
slopes <- (seq(from=1,to=10,length.out=20)-5.5)/4.5
x <- c()
for (i in 1:20) {
for (j in 1:10) {
x <- c(x,seq(from=0,to=slopes[i]*j,length.out=10))
}
}
y <- rep(1:20,each=100)


###################################################
### code chunk number 14: fig4
###################################################
plot(density(x),main="")
for (i in 1:20) {
d <- density(x[which(y==i)])
lines(d$x,0.3*d$y/max(d$y),col=rainbow(n=20,start=0,end=0.66)[21-i])
}
lines(density(x))


###################################################
### code chunk number 15: fig5
###################################################
x2 <- quantileNormalize(x,y,num.bins=20,num.steps=1,mode="discrete")
d2 <- density(x2)
plot(c(min(d2$x), max(d2$x)), c(0, 20), type="p", cex=0, ylab="Density (by bin)", xlab="", main="")
for (i in 1:20) {
d <- density(x2[which(y==i)])
lines(d$x, 0.5*i+0.25*(1+i/20)*20*d$y/max(d$y), col=rainbow(n=20,start=0,end=0.66)[21-i], lty="solid")
}


###################################################
### code chunk number 16: HELP.Rnw:304-307
###################################################
rand <- which(getFeatures(pairs, "TYPE")=="RAND")
msp.rand <- getSamples(pairs, element="exprs")[rand,]
hpa.rand <- getSamples(pairs, element="exprs2")[rand,]


###################################################
### code chunk number 17: HELP.Rnw:310-317
###################################################
msp.rand.med <- apply(msp.rand, 2, median)
msp.rand.mad <- apply(msp.rand, 2, mad)
hpa.rand.med <- apply(hpa.rand, 2, median)
hpa.rand.mad <- apply(hpa.rand, 2, mad)
msp.fail <- msp.rand.med + 2.5*msp.rand.mad
hpa.fail <- hpa.rand.med + 2.5*hpa.rand.mad
hpa.meth <- apply(hpa.rand, 2, quantile, 0.99)


###################################################
### code chunk number 18: HELP.Rnw:320-327
###################################################
norand <- which(getFeatures(pairs, "TYPE")=="DATA")
size <- as.numeric(getFeatures(pairs, "SIZE"))[norand]
msp <- getSamples(pairs, "Brain2", element="exprs")[norand]
hpa <- getSamples(pairs, "Brain2", element="exprs2")[norand]
nofail <- which(msp>msp.fail["Brain2"] | hpa>hpa.fail["Brain2"])
msp.norm <- msp
msp.norm[nofail] <- quantileNormalize(msp[nofail],size[nofail])


###################################################
### code chunk number 19: HELP.Rnw:330-335
###################################################
meth <- which(msp>msp.fail["Brain2"] & hpa<=hpa.meth["Brain2"])
hpa.norm <- hpa
hpa.norm[meth] <- quantileNormalize(hpa[meth],size[meth])
nometh <- which(hpa>hpa.meth["Brain2"])
hpa.norm[nometh] <- quantileNormalize(hpa[nometh],size[nometh])


###################################################
### code chunk number 20: HELP.Rnw:339-344
###################################################
pairs.norm <- pairs
exprs(pairs.norm)[norand, "Brain2"] <- msp.norm
exprs2(pairs.norm)[norand, "Brain2"] <- hpa.norm
getSamples(pairs, element="exprs")[1:5,]
getSamples(pairs.norm, element="exprs")[1:5,]


###################################################
### code chunk number 21: HELP.Rnw:363-367
###################################################
data <- getSamples(pairs,element="exprs2")
seqids <- getFeatures(pairs,'SEQ_ID')
weight <- getSamples(pairs,element="exprs")
combineData(data, seqids, weight, trim=0.2)[1:5,]


###################################################
### code chunk number 22: HELP.Rnw:371-372
###################################################
combineData(pairs, feature.group='SEQ_ID', trim=0.2)[1:5,]


###################################################
### code chunk number 23: fig6
###################################################
plotPairs(pairs,element="exprs2")