# output basal loh and max loh

rm(list=ls());

#postscript("05176.M34.ps",width=8, height=8);

file = "050306.M34.1.tab";

tb= read.table( file, header=T, sep="\t");

row.num = c(1:4, 5,5, 6,6, 7,7, 8:9, 10,10, 11,11, 12,12, 13,13, 14,14 );

tb2 = tb[row.num, ];

labels = names( tb );

#normalize tb2

for( j in 5:13 ) {

for( i in 1: length(tb2[,1]) ) {

tb2[i,j] = tb2[i,j] * tb2[i,2] * tb2[i,3];

}

}

#generate tb.m

n=length(tb2$t) / 2

col.end = 11;

tb.m = data.frame( matrix( nrow=n, ncol= col.end ) ) # mean values

names( tb.m ) = c( labels[ c(1, 5:13) ], "total" );

for( i in 1: length(tb.m[,1]) ) {

tb.m[i,1] = tb2$t[2*i -1]

#tb.sd[i,1] = tb.m[i, 1]

for( j in 2: 10 ) {

tb.m[i, j] = mean( tb2[ c(2*i-1, 2*i), j+3] )

#tb.sd[i,j] = sd( tb2[ c(2*i-1, 2*i), j+3] )

}

tb.m$total[i] = sum( tb.m[i, 2:10], na.rm=T );

}

# output to out,

# columns in out

header = c("t","half.over.black","Pb","R0.5", "R0.75", "s" );

out = data.frame( matrix( nrow= length(tb.m[,1]) , ncol= length(header) ) );

names( out ) = header;

out$t = tb.m$t; # "t"

#calculate the ratios and s

for( i in length(out[,1]):1 ) { #row

out$s[i] = tb.m$total[i] / tb.m$total[1]; #s

out[i,2] = tb.m$B0.5[i] / tb.m$black[i];

}

#calculate the risks (in contrast to percentages)

for( i in length(out[,1]):1 ) { #row

out$Pb[i] = tb.m$black[i] / tb.m$total[i]; # full blacks

out$R0.5[i] = tb.m$B0.5[i] / ( tb.m$total[i] - 2 * tb.m$black[i] ); # half blacks

}

###generate plots

tokens = unlist( strsplit(file, "/") );

shortfile = tokens[1];

#postscript( paste( "_", shortfile, "loh.042506.ps", sep=".") );

par(mar=c(4,4,4,4));

#plot the ratio

y.max = max( out[,2], na.rm=T ) * 1.2;

y.min = min( out[,2], na.rm=T ) ;

cols = c("blue","black","green");

if ( y.max == Inf ) { y.max = 7; };

plot(out[,2] ~ out$t,type='l', col=cols[1],

main= paste( shortfile, " ", "Half vs Full Blacks" ),

ylim=c(y.min, y.max), xlab="t (hours)",ylab="Half/Full Ratio");

#overlay viability plot.

par(new=T);

plot( out$s ~ out$t,type='l',xlab="",ylab="",axes=F,lty=2);

#plot the percentage and risks

y.max = max( out[,3:4], na.rm=T ) * 1.2;

y.min = min( out[,3:4], na.rm=T ) ;

par(new=T);

plot(out$Pb ~ out$t, type='l', col=cols[2],

ylim=c(y.min, y.max), xlab="",ylab="",axes=F);

lines( out[,4] ~ out$t, col=cols[3] );

axis(4, at=pretty(c(y.min, y.max)) );

mtext("Percentage", side=4, line=3);

#legend

labels = c("Half/Full", "P(full)", "R(1/2)", "viability");

ltypes = c(1,1,1,2);

legend( 0, y.max * 0.8, labels, col=c(cols,"black"), lty = ltypes);

#dev.off();

#} ###file loop####

#dev.off();

# q("no");