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