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

gompertz bootstrap fitting, using gnls in R

#092304Thu
 library(nlme)
 library(survival)

## functions
 gompfit <- function( data ) { #data is a vector
   tmp.fit <- survfit(Surv( data ) )
   s <- c( 1, tmp.fit$surv );
   t <- c( 0, tmp.fit$time );
   fm <- gnls( s ~ exp( (I/G)* ( 1 - exp(G * t) )  ) ,
                  start = list( I = 0.001, G = 0.2 ),
                  );
   I.fit <- fm$coefficients[1];
   G.fit <- fm$coefficients[2];
   c( I.fit, G.fit);
 }

 bootstrap.gompfit <- function ( in.vector ) {
   data <- in.vector[ ! is.na(in.vector) ];
   boot.num <- length(data);
   boot.out <- data.frame(matrix(ncol=3,nrow=boot.num));
   names(boot.out) <- c("I","G","avg.rls");
   for ( i in 1:boot.num ) {
     # data.tmp <- BY4741;
     data.tmp <- sample( data, size=length(data), replace=T);
     gf.out <- gompfit(data.tmp);
     boot.out[i, ] <- c( gf.out, mean(data.tmp) );
   }
   ret = boot.out;
 }

## the main section
 rls <- read.table( "rls.tab", sep="\t", header=T);
 strains <- names(rls);

 #buffers for results
 return.labels <- c("I", "se.I", "G", "se.G", "rls", "se.rls" );
 gomp.out <- data.frame( matrix(nrow=length(names(rls)), ncol=length(return.labels) ) );
 names(gomp.out) <- return.labels;
 row.names(gomp.out) <- names(rls);

 #go over each column anb call bootstrap.gompfit
 group.1 = c(1:20 );
 group.2 = c(21:30 );
 group.3 = c(31:length(strains) );
 
 for ( j in group.1 ) {
   boot.array = bootstrap.gompfit( rls[,j] );
   gomp.out[j,c(1,3,5)] = mean(boot.array);
   gomp.out[j,c(2,4,6)] = sd(boot.array);
 }

 #gomp.out;

 for ( j in group.2 ) {
   boot.array = bootstrap.gompfit( rls[,j] );
   gomp.out[j,c(1,3,5)] = mean(boot.array);
   gomp.out[j,c(2,4,6)] = sd(boot.array);
 }

 
 for ( j in group.3 ) {
   boot.array = bootstrap.gompfit( rls[,j] );
   gomp.out[j,c(1,3,5)] = mean(boot.array);
   gomp.out[j,c(2,4,6)] = sd(boot.array);
 }

 

 write.table(gomp.out, "092204.gomp.bootstrap.merged.lab-strains.csv", row.names=T, quote=T, sep="\t");


#q("no")