##########################################################################
# The R script reads a data file of Ks (two columns, Ks and SE)          #  
# and conducts a bootstrap K-S test based on estimated lambda (rate of   #
# gene death                                                             #
##########################################################################
# example of a test
# testdata<-readdata('1.yn00.txt'); 
# lambda<-lambda_est(testdata);
# simulation<-simt(testdata,lambda); # a dataset for simulated Ks 
# pvalue<-test_d(testdata,simulation); 
# pvalue   #<- pvalue output 
###################################################################
library(stats);
library(Matching);
err<-function(ks){
   s.lo<-loess(SE~Ks,ks, control = loess.control(surface = "direct"))
   s.lo
}
trunc<-function(ks,low,high){
	kst<-c();
	j<-0;
	for (i in 1:length(ks$Ks)) {
		if (ks$Ks[i]>low && ks$Ks[i]<high)
		{
		j<-j+1;
		kst$Ks[j]<-ks$Ks[i];
		kst$SE[j]<-ks$SE[i];
		}
	}
	kst
}
readdata<-function(file) {
    #read data
    ks<-read.table(file=file,header=TRUE);
    #remove ks =0.0000 or ks>2
    low<-0.0001;
    high<-2.0;
    kst<-trunc(ks,low,high);
    kst
}
lambda_est<-function(kst) {
	#exponential rate
	n<-length(kst$Ks);
	lambda<-sum(kst$Ks+2*kst$SE)/(2*n);
	lambda
}
simt<-function(kst,lambda){
    suppressWarnings(1);
    size<-length(kst$Ks);
    #lambda<-mean(kst$Ks);
    high<-2.0;
    low<-0.0001;
    #simulate Ks at Log scale
    U<-runif(1000); #generate a much larger sample
    T<--lambda*log(U); # T ~ exponential(lambda)  
    s.Ks<-c();
    j<-0;
    k<-0;
    t<-c();
    while (j<size) {
        k<-k+1;
     	if (T[k]<high && T[k]> low) {
           t<-c(t,T[k]);
           j<-j+1;
        };
   };
   s.T<-predict(err(kst),data.frame(t),se=F);
   z<-rnorm(size);
   s.Ks<-t+s.T*z;
   abs(s.Ks);
}   
test_d<-function(kst,s.Ks){
	   suppressWarnings(TRUE);
	   size<-min(length(kst$Ks), length(s.Ks));
        #do a ks test and get the p-value
        test<-ks.boot(kst$Ks[1:size],s.Ks[1:size],nboots=1000);
        test$ks.boot.pvalue
}