Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
siRNA screen scatterplots
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.