siRNA screen scatterplots

CompareCompanies.R

getReps = function(df,avedef){
  dfA = df[df$RepeatNumber%in%c(1,3,5),]
  dfB = df[df$RepeatNumber%in%c(2,4,6),]
  dfaveA = aggregate(dfA,by=list(dfA$Gene),FUN=avedef)
  dfaveA$Gene = dfaveA$Group.1
  dfaveB = aggregate(dfB,by=list(dfB$Gene),FUN=avedef)
  dfaveB$Gene = dfaveB$Group.1
  dfreps = merge(dfaveA,dfaveB,by="Gene")
  return(dfreps)
}

makeQ = function(df,method="fdr",test="wilcox.test"){
  if(test=="wilcox.test") {summfunc = median; test = wilcox.test}else{summfunc = mean; test = t.test}
  ctrl = aggregate(df[[var]],by=list(paste(df$Plate,df$RepeatNumber,sep="_")),FUN=summfunc)$x
  getp = function(x) test(x,ctrl)$p.value
  dfp = aggregate(df[[var]],by=list(df$Gene),FUN=getp)
  colnames(dfp) = c("Gene","p")
  dfp$q = p.adjust(dfp$p,method)
  q = dfp$q
  names(q)=dfp$Gene
  return(q)
}

dharmapath="IRRADIATION\\ANALYSISOUT\\IRRADIATION_RawData.txt"
sigmapath="IRRADIATION_SIGMA\\ANALYSISOUT\\IRRADIATION_SIGMA_RawData.txt"

dharma=read.delim(dharmapath,sep="\t",stringsAsFactors=FALSE)
sigma=read.delim(sigmapath,sep="\t",stringsAsFactors=FALSE)

dharma = dharma[dharma$Treatment=="IRR",]
sigma = sigma[sigma$Treatment=="IRR",]

common = intersect(unique(dharma$Gene),unique(sigma$Gene))

dharma = dharma[dharma$Gene%in%common,]
sigma = sigma[sigma$Gene%in%common,]

avedef = median
var = "NormalisedCount"
dharma[[var]] = dharma[[var]]/1000
sigma[[var]] = sigma[[var]]/1000
axrng = c(1,60)

dharmave = aggregate(dharma,by=list(dharma$Gene),FUN=avedef)
dharmave$Gene = dharmave$Group.1
dharmave = dharmave[order(dharmave[[var]]),]
dharmaq = makeQ(dharma,method="fdr")
dharmave$q = dharmaq[dharmave$Gene]

sigmave = aggregate(sigma,by=list(sigma$Gene),FUN=avedef)
sigmave$Gene = sigmave$Group.1
sigmave = sigmave[order(sigmave[[var]]),]
sigmaq = makeQ(sigma,method="fdr")
sigmave$q = sigmaq[sigmave$Gene]

n=20
sigmahits = c(head(sigmave$Gene[sigmave$q<0.05],n=n),tail(sigmave$Gene[sigmave$q<0.05],n=n))
dharmahits = c(head(dharmave$Gene[dharmave$q<0.05],n=n),tail(dharmave$Gene[dharmave$q<0.05],n=n))
print(intersect(sigmahits,dharmahits))

ss = sigma[sigma$Gene%in%c("EMPTY",sigmahits),]
dd = dharma[dharma$Gene%in%c("EMPTY",dharmahits),]
op=par(mfrow=c(2,1))
stripchart(log(NormalisedCount)~Gene,data=ss,method="jitter",vertical=TRUE,las=2,main="Sigma-Aldrich")
stripchart(log(NormalisedCount)~Gene,data=dd,method="jitter",vertical=TRUE,las=2,main="Dharmacon")
par(op)

#sigmahits = sigmave$Gene[sigmave$q<0.05]
#dharmahits = dharmave$Gene[dharmave$q<0.05]

both = merge(dharmave,sigmave,by="Gene",suffixes=c(".dharma",".sigma"))

dharmareps = getReps(dharma,avedef)
sigmareps = getReps(sigma,avedef)

png("Correlations.png",width=3000,height=1000,pointsize=60)
op=par(mfrow=c(1,3),mar=c(5, 5, 2, 1) + 0.1)

cstr = formatC(cor(log(dharmareps[[paste(var,"x",sep=".")]]),log(dharmareps[[paste(var,"y",sep=".")]])),2)
plot(NULL,main=paste("Correlation",cstr),xlab="Median fluorescence (replicates 1, 3 & 5)\nDharmacon",ylab="Median fluorescence (replicates 2, 4 & 6)\nDharmacon",xlim=axrng,ylim=axrng,log="xy")
abline(a=0,b=1,lty=2,lwd=2,col="black")
points(dharmareps[[paste(var,"x",sep=".")]],dharmareps[[paste(var,"y",sep=".")]],col=rgb(0,0,0,0.2),pch=16)
hitpts = dharmareps[dharmareps$Gene%in%dharmahits,]
neut = dharmareps[dharmareps$Gene=="EMPTY",]
points(hitpts[[paste(var,"x",sep=".")]],hitpts[[paste(var,"y",sep=".")]],col="red",pch=1,cex=1,lwd=2)
points(neut[[paste(var,"x",sep=".")]],neut[[paste(var,"y",sep=".")]],col="yellow",pch=16,cex=1)

cstr = formatC(cor(log(sigmareps[[paste(var,"x",sep=".")]]),log(sigmareps[[paste(var,"y",sep=".")]])),2)
plot(NULL,main=paste("Correlation",cstr),xlab="Median fluorescence (replicates 1, 3 & 5)\nSigma-Aldrich",ylab="Median fluorescence (replicates 2, 4 & 6)\nSigma-Aldrich",xlim=axrng,ylim=axrng,log="xy")
abline(a=0,b=1,lty=2,lwd=2,col="black")
points(sigmareps[[paste(var,"x",sep=".")]],sigmareps[[paste(var,"y",sep=".")]],col=rgb(0,0,0,0.2),pch=16)
hitpts = sigmareps[sigmareps$Gene%in%sigmahits,]
neut = sigmareps[sigmareps$Gene=="EMPTY",]
points(hitpts[[paste(var,"x",sep=".")]],hitpts[[paste(var,"y",sep=".")]],col="blue",pch=16,cex=0.5)
points(neut[[paste(var,"x",sep=".")]],neut[[paste(var,"y",sep=".")]],col="yellow",pch=16,cex=1.0)

sig=both[[paste(var,"sigma",sep=".")]]
dha=both[[paste(var,"dharma",sep=".")]]
cstr = formatC(cor(log(both[[paste(var,"sigma",sep=".")]]),log(both[[paste(var,"dharma",sep=".")]])),2)
plot(NULL,main=paste("Correlation",cstr),xlab="Median fluorescence (all replicates)\nSigma-Aldrich",ylab="Median fluorescence (all replicates)\nDharmacon",xlim=axrng,ylim=axrng,log="xy")
abline(v=sig[both$Gene=="EMPTY"],h=dha[both$Gene=="EMPTY"],col=c("red","blue"),lty=3,lwd=3)
abline(a=0,b=1,lty=2,lwd=2,col="black")
points(sig,dha,col=rgb(0,0,0,0.2),pch=16)
points(sig[both$Gene%in%sigmahits],dha[both$Gene%in%sigmahits],col="blue",pch=16,cex=0.5)
points(sig[both$Gene%in%dharmahits],dha[both$Gene%in%dharmahits],col="red",pch=1,cex=1,lwd=2)
points(sig[both$Gene=="EMPTY"],dha[both$Gene=="EMPTY"],col="yellow",pch=16,cex=1.0)

par(op)
dev.off()

#boxplot(log(Count)~Gene,data=dharma)