Sunday, May 21, 2006

Over full and half LOH plot in R

# 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");