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

###################################################
### code chunk number 1: coRNAi.Rnw:46-57
###################################################
makefig = function(expr, name, type="pdf", ppi=180, width=4.5, height=4.5){
  fn = paste(name, type, sep=".")
  if(!file.exists(fn)){
    switch(type,
       pdf = pdf(fn, width=width, height=height),
       jpeg = jpeg(fn, width=width*ppi, height=height*ppi, pointsize=24),
       stop(sprintf("Invalid type '%s'", type)))
    expr
    dev.off()
    }
  }


###################################################
### code chunk number 2: load libraries
###################################################
library("coRNAi")


###################################################
### code chunk number 3: load data
###################################################
data(screen1_raw)
data(screen2_raw)
#data(num2name)


###################################################
### code chunk number 4: make into data frame
###################################################
dfR1 = cellHTS2df(screen1_raw,neutral=c("Fluc"))
dfR2 = cellHTS2df(screen2_raw,neutral=c("Fluc"))
dfR1$value = log2(dfR1$value)
dfR2$value = log2(dfR2$value)


###################################################
### code chunk number 5: coRNAi.Rnw:121-122
###################################################
dfR2[1:2,]


###################################################
### code chunk number 6: coRNAi.Rnw:132-144
###################################################
mypars = list(cex.lab = 1,cex.main=1)
par(mfrow=c(1,2))

BoxPlotShorth(value~replicate,dfR1[dfR1$controlStatus=="sample",],
              las=2,main="",xlab="plates",boxfill="lightblue",
              outline=FALSE, ylab=expression(log[2](intensity)),pars=mypars)

BoxPlotShorth(value~plate,dfR2[dfR2$controlStatus=="sample" & 
                               dfR2$replicate==1,], las=2,
              main="",xlab="plates",
              boxfill="lightblue",outline=FALSE, 
              ylab=expression(log[2](intensity)),pars=mypars)


###################################################
### code chunk number 7: normaliztion methods
###################################################
dfScl1 = cellHTS2df(normalizePlates(screen1_raw,method="shorth",
  scale="multiplicative",log=TRUE),neutral="Fluc")
dfScl2 = cellHTS2df(normalizePlates(screen2_raw,method="shorth",
  scale="multiplicative",log=TRUE),neutral="Fluc")


###################################################
### code chunk number 8: coRNAi.Rnw:162-170
###################################################
par(mfrow=c(1,2))
BoxPlotShorth(value~replicate,dfScl1[dfScl1$controlStatus=="sample",],
              las=2, main="",xlab="plates",boxfill="lightpink",
              outline=FALSE,ylab=expression(log[2](intensity)))
BoxPlotShorth(value~plate,dfScl2[dfScl2$controlStatus=="sample"&
                                 dfScl2$replicate==1,],
              boxfill="lightpink",las=2,main= "",
              xlab="plates",outline=FALSE,ylab=expression(log[2](intensity)))


###################################################
### code chunk number 9: coRNAi.Rnw:184-187
###################################################
plotScreen(split(dfScl1$value,dfScl1$replicate),fill = 
           c("red","white","blue"),zrange=c(-4,4),ncol=5,
            main="screen 1",legend.label=NULL)


###################################################
### code chunk number 10: coRNAi.Rnw:198-202
###################################################
plotScreen(split(dfScl2$value[dfScl2$replicate==1],
                 dfScl2$plate[dfScl2$replicate==1]),fill = 
           c("red","white","blue"),zrange=c(-9,9),
           main="screen 2",legend.label=NULL)


###################################################
### code chunk number 11: coRNAi.Rnw:213-225
###################################################

par(mfrow=c(1,2))
WithinScreenPlot(dfScl1,what="value",main="",smooth=T,pch=".")
WithinScreenPlot(dfScl2,what="value",main="",smooth=T,pch=".")
points(dfScl2$value[dfScl2$Pair=="P2 P1" & dfScl2$replicate==2 & 
                    dfScl2$Direction==1],dfScl2$value[dfScl2$Pair=="P2 P1" & 
                      dfScl2$replicate==2 & dfScl2$Direction==2],
       col="red",pch=19)
points(dfScl2$value[dfScl2$Pair=="P57 P1" & dfScl2$replicate==2 &
                    dfScl2$Direction==1],dfScl2$value[dfScl2$Pair=="P57 P1" & 
                      dfScl2$replicate==2 & dfScl2$Direction==2],
       col="red",pch=19)


###################################################
### code chunk number 12: show data
###################################################
dfScl2[dfScl2$Pair =="P2 P1",c("Pair","Direction","replicate","value")]
dfScl2$value[dfScl2$Pair =="P2 P1" & dfScl2$Direction==1 & 
             dfScl2$Replicate==2]=NA


###################################################
### code chunk number 13: coRNAi.Rnw:242-252
###################################################
data(faultyscreen)
par(mfrow=c(1,2))
fdf = cellHTS2df(normalizePlates(faultyscreen,method="shorth",
  scale="multiplicative",log=TRUE),neutral="Fluc")
WithinScreenPlot(fdf[fdf$replicate==2,],main="with plate column 13"
                 ,pch=".",smooth=FALSE)
systerr = which(fdf$Pair%in%(unique(fdf$Pair[grep(13,fdf$well)])))
fdf$value[systerr]=NA
WithinScreenPlot(fdf[fdf$replicate==2,],main="without plate column 13"
                 ,pch=".",smooth=FALSE)


###################################################
### code chunk number 14: coRNAi.Rnw:261-262
###################################################
BetweenScreenPlot(dfScl1,smooth=FALSE)


###################################################
### code chunk number 15: add batch info
###################################################
dfScl1$batch[dfScl1$replicate%in%1:5]=1
dfScl1$batch[dfScl1$replicate%in%6:10]=2


###################################################
### code chunk number 16: add weights
###################################################
dfScl1 = weightDf(dfScl1,exclude = c("controlN2", "controlP2","controlP1N1"
                           ,"controlP1","double","controlN1"))

dfScl2 = weightDf(dfScl2,exclude = c("controlN2", "controlP2","controlP1N1"
                           ,"controlP1","double","controlN1"))


###################################################
### code chunk number 17: main effects estimate screen1
###################################################
lmres1 = lmmain(dfScl1,per = "batch") 


###################################################
### code chunk number 18: main effects estimate screen2
###################################################
lmres2D= lmmain(dfScl2,per = c("replicate","Direction"))
lmres2 = lmmain(dfScl2,per = c("replicate"))


###################################################
### code chunk number 19: coRNAi.Rnw:310-314
###################################################

barplot(t(sapply(lmres1,function(x) x$coefficient)),las=2,beside=T,
        col=c("skyblue","steelblue"))
legend("topleft",legend=c("batch1","batch2"),fill=c("skyblue","steelblue"))


###################################################
### code chunk number 20: update
###################################################
dfScl1 = updateDf(dfScl1,lmres1,per = "batch")
dfScl2 = updateDf(dfScl2,lmres2,per=c("replicate"))
head(dfScl1)


###################################################
### code chunk number 21: coRNAi.Rnw:335-341
###################################################

par(mfrow=c(1,2))
MainFitPlot(lmres1,xlab=expression(hat(y)[0*k]+hat(m)[i*k]+hat(m)[j*k]),
            ylab=expression(epsilon[i*j*k]),pch=".",main="screen 1")
MainFitPlot(lmres2,xlab=expression(hat(y)[0*k]+hat(m)[i*k]+hat(m)[j*k]),
            ylab=expression(epsilon[i*j*k]),pch=".",main="screen 2",sd.fit=FALSE)


###################################################
### code chunk number 22: coRNAi.Rnw:353-359
###################################################
pairs=sample(unique(dfScl1$Pair[dfScl1$Type=="comb"]),50)
sub = dfScl1[dfScl1$Pair%in%pairs,]
boxplot(residuals~Pair,sub,las=2,axis.cex=0.5,col ="palevioletred",
        ylab=expression(epsilon[i*j*k]),xlab="RNAi combinations",
        names = NA,xaxt="n")
abline(h=0,lwd=2)


###################################################
### code chunk number 23: ebayes
###################################################
ebfit = df2lmFit(dfScl1)
eb = eBayes(ebfit)  
tt1 =interactiontable(eb,ord.t=TRUE) 
ebfit = df2lmFit(dfScl2)
eb = eBayes(ebfit)  
tt2 =  interactiontable(eb)
head(tt2)


###################################################
### code chunk number 24: coRNAi.Rnw:388-389
###################################################
Pplot(tt2$P,maint = "Screen 2")


###################################################
### code chunk number 25: read key
###################################################
data(key)
colkey = ifelse(key$cellCycle==0,"grey","orange")
names(colkey)=as.character(key$GeneID)
key = key$cellCycle
names(key) = names(colkey)


###################################################
### code chunk number 26: threshold
###################################################
thrs1 = 0.001


###################################################
### code chunk number 27: coRNAi.Rnw:415-423
###################################################
tt1t = tt1
thrs1 = 0.1

g1= data2graph(tt1,thres=thrs1,thresBy="ord.p.adj",nodecolor=colkey,
  sizethres=0.3, writedot=TRUE, filename = 'interactionGraph1.dot', 
  shape='ellipse',scaleFactor=10,fixedsize=FALSE,fontsize=100,
  penwidth=20,width=2,gamma.col=0)
system('neato -Tpdf -o interactionGraph1.pdf interactionGraph1.dot')


###################################################
### code chunk number 28: coRNAi.Rnw:429-433
###################################################
g2 = data2graph(tt2,thres=thrs1,writedot=TRUE, scaleFactor=10,
  filename = 'interactionGraph2.dot',fontsize=20,penwidth=10,
  width=2,fixedsize=TRUE,nodecolor="grey",sizethres=0.3,gamma.col=0)
system('neato -Tpdf -o interactionGraph2.pdf interactionGraph2.dot')


###################################################
### code chunk number 29: coRNAi.Rnw:438-447
###################################################
mat = tt2matrix(tt1,what="size")
thrs=0.9
cormat = cortestmatrices(mat,method="spearman")
cormat[[2]][abs(cormat[[1]])<0.8]=1
c1 = data2graph(list("size"=cormat[[1]],"pvalues"=cormat[[2]]),
  thres=thrs,writedot=TRUE, scaleFactor=2, filename = 'correlationGraph1.dot',
  shape='ellipse',fixedsize=FALSE,nodecolor=colkey,fontsize=30,penwidth=5,
  gamma.col=0)
system('neato -Tpdf -o correlationGraph1.pdf correlationGraph1.dot')


###################################################
### code chunk number 30: coRNAi.Rnw:453-459
###################################################
mat = tt2matrix(tt2,what="size")
cormat = cortestmatrices(mat,method="spearman")
c2 = data2graph(list("size"=cormat[[1]],"pvalues"=cormat[[2]]),thres=
  thrs1,writedot=TRUE, scaleFactor=20, filename = 'correlationGraph2.dot',
  width=2,fontsize=50,penwidth=4,fixedsize=TRUE,nodecolor="grey")
system('neato -Tpdf -o correlationGraph2.pdf correlationGraph2.dot')


###################################################
### code chunk number 31: coRNAi.Rnw:482-484
###################################################
InteractLevelPlot(tt1,key = key,thresh=thrs1,by="P.Value",zerolimit=0.1,
                  colorRampPalette(c("blue", "white", "red")))


###################################################
### code chunk number 32: coRNAi.Rnw:493-498
###################################################

PlotHeatmap(tt2,dendrogram="none",margins = c(2,2), 
            colorRampPalette(c("blue", "white","red")),
            lmat=rbind( c(0, 3), c(2,1), c(0,4) ), 
            lhei=c(0.1, 4, 0.1 ) ,lwid=c(0.1,2))


###################################################
### code chunk number 33: fitness scale
###################################################
N0=15000
sc2 = screen2_raw
Data(sc2) = log(Data(sc2)/N0)


###################################################
### code chunk number 34: reanalys
###################################################
ksdf2 = cellHTS2df(normalizePlates(sc2,method="shorth",
  scale="multiplicative",log=TRUE),neutral="Fluc")
ksdf2 = weightDf(ksdf2,exclude = c("controlN2", "controlP2",
                         "controlP1N1","controlP1","double"))
lmres2 = lmmain(ksdf2,per = c("replicate","Direction"))
ksdf2 = updateDf(ksdf2,lmres2,per = c("replicate","Direction"))


###################################################
### code chunk number 35: coRNAi.Rnw:546-549
###################################################
par(mfrow=c(1,2))
qqnorm(ksdf2$residuals,main="log(log(N/N0))",pch=".")
qqnorm(dfScl2$residuals,main="log(N)",pch=".")


###################################################
### code chunk number 36: coRNAi.Rnw:559-560
###################################################
ROCstat=coRNAi:::ROCstat


###################################################
### code chunk number 37: coRNAi.Rnw:579-654
###################################################
## First find bench mark results
ResTable=tt1
par(mfrow=c(1,2))
t = 0.001 # p-value threshold = 0.001

bench.neg.int = ResTable[ResTable$ord.p<t & ResTable$size<0,]
bench.pos.int = ResTable[ResTable$ord.p<t & ResTable$size>0,]

## function for extracting p values from moderated t-test and ordinary t-test on subset of data

GetStats = function(df2){
  df = length(unique(df2$replicate))*2-1
  lm2 = lmmain(df2)
  df2=updateDf(df2,lm2)
  ebfit2 = df2lmFit(df2)
  ebfit2 = eBayes(ebfit2)
  resTable = interactiontable(ebfit2,ord.t=TRUE)
  resTable
}

ordstat=array(NA,c(3,101,10))
modstat= array(NA,c(3,101,10))
logstat = array(NA,c(3,101,10))
resTT = list()
for (i in 1:10){
  dfsub = dfScl1[dfScl1$replicate%in%i,]
  resTT = GetStats(dfsub)
  thr = seq(0,1,0.01)
  ordstat[,,i]=sapply(thr,function(x) ROCstat(thrs=x,stats="ord.p",testSet=resTT, trueNegSet = bench.neg.int, truePosSet = bench.pos.int))
  modstat[,,i]=sapply(thr,function(x) ROCstat(thrs=x,stats="mod.p",testSet=resTT, trueNegSet = bench.neg.int, truePosSet = bench.pos.int))
  thr = seq(1,0,-0.01)
  logstat[,,i]=sapply(thr,function(x) ROCstat(thrs=x,stats="size",testSet=resTT, trueNegSet = bench.neg.int, truePosSet = bench.pos.int))
}

plotROC = function(ordstat,modstat,logstat,main="",...){

  meanord = apply(ordstat,1:2,mean,na.rm=T)
  meanmod = apply(modstat,1:2,mean,na.rm=T)
  meanlog = apply(logstat,1:2,mean,na.rm=T)

  col1 = hsv(0.8,0.7,0.7,0.4)
  col2 = hsv(0.3,.7,.5,0.4)
  col3 = hsv(0.6,.7,.7)
  
  plot(NULL,ylim=c(0,1),xlim=c(0,0.3),ylab = "true positive rate", xlab = "false positive rate",main=main,...)
  
  lines(meanord[2,],meanord[1,],col=col1,lwd=4)
  lines(meanmod[2,],meanmod[1,],col=col2,lwd=4)
  lines(meanlog[2,],meanlog[1,],col=col3,lwd=4)
  
  
  legend("bottomright", legend = c("moderated t","normal t","effect size"), col= c(col2,col1,col3),lwd=4)
}

plotROC(ordstat,modstat,logstat,main="a")

ordstat=array(NA,c(3,101,25))
modstat= array(NA,c(3,101,25))
logstat = array(NA,c(3,101,25))
cc = combn(1:10,2)
r1= which(cc[1,]<6)
r2 =which(cc[2,]>5)
combs = cc[,intersect(r1,r2)]

for ( i in 1:ncol(combs)){
  dfsub = dfScl1[dfScl1$replicate%in%combs[,i],]
  
  resTT = GetStats(dfsub)
  thr = seq(0,1,0.01)
  ordstat[,,i]=sapply(thr,function(x) ROCstat(thrs=x,stats="ord.p",testSet=resTT, trueNegSet = bench.neg.int, truePosSet = bench.pos.int))
  modstat[,,i]=sapply(thr,function(x) ROCstat(thrs=x,stats="mod.p",testSet=resTT, trueNegSet = bench.neg.int, truePosSet = bench.pos.int))
  thr = seq(1,0,-0.01)
  logstat[,,i]=sapply(thr,function(x) ROCstat(thrs=x,stats="size",testSet=resTT, trueNegSet = bench.neg.int, truePosSet = bench.pos.int))
}
plotROC(ordstat,modstat,logstat,main="b")


###################################################
### code chunk number 38: coRNAi.Rnw:674-683
###################################################
lmm = lmmain(dfScl2)
mm = estmodel(dfScl2,estimate="median")
par(mfrow=c(1,2))
plot(mm$coefficient,lmm$coefficient,pch=".",ylab="OLS estimates",
     xlab="median estimates",main="main effects")
abline(0,1,col="red")
plot(mm$residuals,lmm$residuals,pch=".",ylab="OLS estimates",
     xlab="median estimates",main="residuals")
abline(0,1,col="red")


###################################################
### code chunk number 39: coRNAi.Rnw:693-696
###################################################
par(mfrow=c(1,2))
MainFitPlot(mm,main="fitted by estmodel",pch=".")
MainFitPlot(lmm,main = "fitted by lmmain",pch=".")


###################################################
### code chunk number 40: sessionInfo
###################################################
sessionInfo()