#YPS163 2011April 13, plot
rm=(list=ls())
setwd("~/github/LOH_H2O2_2012-master/figures")
debug = 0;
FileList = list.files( path="../data.H2O2-LOH/", pattern="YPS163"); FileList;
filename = 'YPS163,Apr13,2011.H2O2.LOH.csv';
fullFileName = paste('../data.H2O2-LOH/',filename, sep='');
mylabel = 'YPS163 20110413'
tb = read.csv(fullFileName, colClasses=c("character",NA, NA, "character", rep("numeric",8 ), NA));
names(tb) = c("Strain", "OD600", "Dilution","Date","H2O2stock", "White", "Black", "halfBlack", "quarterBlack", "ThreeQBlack", "QQBlack", "Other", "Notes")
######## set zeros
mycolumns = c("White","Black","halfBlack", "quarterBlack","ThreeQBlack", "QQBlack", "Other");
for( i in 1:length(tb[,1])) {
for ( j in mycolumns) {
if( is.na(tb[i,j]) ) { tb[i,j]= 0 }
}
}
#adjust for low-counts
tb$Black[tb$Black<=0] = 0.5
tb$halfBlack[tb$halfBlack<=0] = 0.5
tb$Black[tb$Black<0]=NA; #remove weird experimental data, such as low-lead concentration effect
tb$H2O2 = tb$H2O2stock/2
tb$tot = tb$White + tb$Black + tb$halfBlack + tb$quarterBlack + tb$ThreeQBlack + tb$QQBlack + tb$Other
tb.ori = tb;
tb = tb[ ! is.na(tb$White), ]
tb$Dilution = tb$Dilution / tb$Dilution[1]
######## normalize all data
mycolumns = c("White","Black","halfBlack", "quarterBlack","ThreeQBlack", "QQBlack", "Other","tot");
for ( j in mycolumns) {
tb[,j] = tb[,j] * tb$Dilution
}
####### find out means, sd
H2O2 = sort(unique( tb$H2O2))
#s = H2O2
tbm = data.frame(cbind(H2O2))
for ( i in 1:length(H2O2)) {
c = H2O2[i]
tmp = tb[ tb$H2O2==c, ]
tbm$tot[i] = mean(tmp$tot, na.rm=T)
tbm$tot.sd[i] = sd(tmp$tot, na.rm=T)
tbm$White[i] = mean(tmp$White, na.rm=T)
tbm$Black[i] = mean(tmp$Black, na.rm=T)
tbm$Black.sd[i] = sd(tmp$Black, na.rm=T)
tbm$halfBlack[i] = mean(tmp$halfBlack, na.rm=T)
tbm$halfBlack.sd[i] = sd(tmp$halfBlack, na.rm=T)
tbm$quarterBlack[i] = mean(tmp$quarterBlack, na.rm=T)
tbm$ThreeQBlack[i] = mean(tmp$ThreeQBlack, na.rm=T)
tbm$QQBlack[i] = mean(tmp$QQBlack, na.rm=T);
}
tbm = tbm[tbm$tot>2, ] #remove plates with <2 colonies. At least 1 white and 1 black is needed.
###### some manual curations here
# ...
###### calculate fractions
tbf = tbm;
tbf$s = tbf$tot / max(tbf$tot)
for ( j in c("Black","halfBlack", "quarterBlack","ThreeQBlack","QQBlack")) {
tbf[, j] = tbf[,j] / tbf$tot
}
for ( j in c("tot.sd", "Black.sd","halfBlack.sd")) {
tbf[, j] = tbf[,j] / max(tbf$tot); #bug fixed, 2013 Dec 10, 18:50
tbf[, j] = ifelse(tbf[, j]<=1/700, 1/700, tbf[, j] ) #at least 0.2% error
}
head(tbf)
plot(tbf$tot.sd ~ tbf$H2O2)
#It is not clear whether color figures would be charged more or not by Aging Cell.
#### define functions
logistical.viability <- function( T, w, t ) { ret <- 1 /( 1 + ( t / T )^ w ); }
logistical.black <- function(b.max, b.min, T, w, t ) { ret <- b.max - (b.max - b.min) /( 1 + ( t / T )^ w ); }
# genome.integrity <- function(b.max, b.min, T, w, t) { 2 * (b.max - b.min) / (1 + (t/T)^w) + (1 - 2 * b.max); }
######## calclulate Cv, Cb
require(nlme)
t= tbf$H2O2; s= tbf$s; b = tbf$Black
model1 = function(v,w) { 1/( 1 + ( t / v )^ w ) }
fm.s <- gnls( s ~ model1(v,w) , start = list( v = 0.02, w=2) );
t2 = seq(0,1,by=0.001)
fit.s = logistical.viability ( fm.s$coefficients[1], fm.s$coefficients[2], t2 );
Cv = fm.s$coefficients[1]
#estimate starting values
b.tmp = sort( b[ b>0 ] );
b.min = sum(b.tmp[1])/1
b.tmp = rev(b.tmp);
b.max = sum(b.tmp[1:3])/3
modelBlack = function(v, w) { b.max - (b.max-b.min)/(1 + (t/v)^w) }
modelBlack1 = function(v) { b.max - (b.max-b.min)/(1 + (t/v)) }
#fm.b <- gnls( b ~ modelBlack1(v) , start = list( v = 0.05) );
fm.b <- gnls( b ~ modelBlack(v,w) , start = list( v = 0.005, w= 4 )); fm.b
fit.b = logistical.black(b.max, b.min, fm.b$coef[1], fm.b$coef[2], t2 )
Cb = fm.b$coef[1]
####################################################################
# do the plot
####################################################################
#### specificy the layout
#pdf(paste("_20131210", filename, "pdf", sep="."), width=6,height=9);
pdf(paste("_20131210", filename, "pdf", sep="."), width=5,height=8);
#tiff(paste("_20131210", filename, "tif", sep="."), width=240,height=500);
mat = matrix( seq(1,3), nrow=3, ncol= 1 );
layout(mat, heights= c( 1.15, 1, 1.2) );
### viability
par(mar=c(0,5,5,1))
plot( tbf$s ~ tbf$H2O2, col="blue", axes=F, xlab='', ylab='Viability', ylim=c(-0.1, 1.2), type='p', pch=19);
axis( 1, at = pretty(tbf$H2O2), labels=F, tcl=0.3);
axis( 2, at=pretty(tbf$s), tcl=0.2, las=2 )
lines( fit.s ~ t2, col="blue",lty=1);
#legend ( max(H2O2)*0.5, 0.5, c("viability","black", "half black"), col=c("blue","black"), lty=c(2,1), pch=c(1,16) )
arrows( tbf$H2O2, (tbf$s - tbf$tot.sd), tbf$H2O2, (tbf$s + tbf$tot.sd), length=0.1, angle=90,code=3, lty=2 );
#legend ( max(H2O2)*0.5, 0.5, c("viability"), col=c("blue"), lty=c(1), pch=c(16) )
points( Cv, 0.5, pch=15, col="red", cex=1.2 );
arrows( Cv, 0.5, Cv, -1, lty=2, col="red");
mtext( "Cv",side=1,at=c(Cv*0.85), line=-1, cex=0.8 );
box()
title(mylabel)
############ full blacks
par(mar=c(0,5,0,1))
plot( tbf$Black ~ tbf$H2O2, pch=16, xlab='',ylab="Black", log='y', ylim=c(0.8E-3, 0.15), axes=F)
axis( 2, at = c(0.001,0.01,0.05,0.1), labels=c(0.001,0.01,0.05,0.1), tcl=0.2, las=2);
axis( 1, at = pretty(tbf$H2O2), labels=F, tcl=0.3);
arrows( tbf$H2O2, (tbf$Black - tbf$Black.sd), tbf$H2O2, (tbf$Black + tbf$Black.sd), length=0.1, angle=90,code=3, lty=2 );
lines( fit.b ~ t2, lty=1, col='black');
points ( Cb, ( b.max/2 + b.min/2), pch=15, col="red", cex=1.2);
arrows( Cb, (b.max/2 + b.min/2), Cb, 1E-8, lty=2, col="red", length=0.1);
mtext( "Cb",side=1,at=c(Cb*1.2), line=-1, cex=0.8 );
box()
####################half black plots
#estimate starting values
b0.5 = tbf$halfBlack
# Do I need to adjust the dilution bias in [0.02,0.03]? No. Plot as is.
# b0.5[6:8] = b0.5[6:8] + 0.01
b0.5.tmp = sort( b0.5[ b0.5>0 ] );
b0.5.min = sum(b0.5.tmp[1])*1.3
b0.5.tmp = rev(b0.5.tmp);
b0.5.max = sum(b0.5.tmp[1])*1.05 #account for increasing trend
modelBlack0.5B = function(v) { b0.5.max - (b0.5.max-b0.5.min)/(1 + (t/v)^2) }
#fm.b0.5 <- gnls( b0.5 ~ modelBlack0.5B(v) , start = list( v = 0.05) ); fm.b0.5
#fit.b0.5 = logistical.black(b0.5.max, b0.5.min, fm.b0.5$coef[1], 2, t2 )
#logistical.black <- function(b.max, b.min, v, w, t ) { ret <- b.max - (b.max - b.min) /( 1 + ( t / T )^ w ); }
modelBlack0.5 = function(v, w) { b0.5.max - (b0.5.max-b0.5.min)/(1 + (t/v)^w) }
fm.b0.5 <- gnls( b0.5 ~ modelBlack0.5(v,w) , start = list( v = 0.008, w=1.5) ); fm.b0.5
fit.b0.5 = logistical.black(b0.5.max, b0.5.min, fm.b0.5$coef[1], fm.b0.5$coef[2], t2 )
#logistical.blackC <- function(b.max, b.min, v, w, u, t ) { ret <- b.max - (b.max - b.min) /( u + ( t / v )^ w ); }
#modelBlack0.5C = function(v, w, u) { b0.5.max - (b0.5.max-b0.5.min)/(u + (t/v)^w) }
#fm.b0.5 <- gnls( b0.5 ~ modelBlack0.5C(v,w,u) , start = list( v = 0.008, w=1.2, u=1.1) ); fm.b0.5
#fit.b0.5 = logistical.blackC(b0.5.max, b0.5.min, fm.b0.5$coef[1], fm.b0.5$coef[2], fm.b0.5$coef[3], t2 )
#plot(b0.5 ~ t )
#lines(fit.b0.5 ~ t2)
Cb0.5 = fm.b0.5$coef[1]
#par(new=T)
par(mar=c(5,5,0,1))
#plot( tbf$halfBlack ~ tbf$H2O2, pch=16, xlab='H2O2',ylab="half-black", col='red')
plot( b0.5 ~ t, pch=16, xlab='H2O2',ylab="Half-black", col='black', log='y', ylim=c(0.8E-3,0.15), axes=F)
axis( 1, at = pretty(t), labels=T, tcl=0.3);
#axis( 2, at = c(0.001,0.01,0.05,0.1), labels=T, tcl=0.2, las=2);
axis( 2, at = c(0.001,0.01,0.05,0.1), labels=c(0.001,0.01,0.05,0.1), tcl=0.2, las=2);
box()
lines(fit.b0.5 ~ t2, col='black')
arrows( tbf$H2O2, (tbf$halfBlack - tbf$halfBlack.sd), tbf$H2O2, (tbf$halfBlack + tbf$halfBlack.sd),
length=0.1, angle=90,code=3, lty=2, lwd=1 );
points ( Cb0.5, ( b0.5.max/2 + b0.5.min/2), pch=15, col="red", cex=1.2);
arrows( Cb0.5, (b0.5.max/2 + b0.5.min/2), Cb0.5, 1E-8, lty=2, col="red", length=0.1);
mtext( "Cb0.5",side=1,at=c(Cb*1.2),line=-1, cex=0.8 );
dev.off()
