# test all genes for survival diff
library(survival);
#library(multtest);

Data = read.delim("Data.txt", header=TRUE, sep="\t");

Filter=which(apply(Data[1,],2,function(x)all(is.na(x))));

if(length(Filter)!=0)
	Data = Data[,-Filter];

GeneCnt = dim(Data)[1]-2;
SampleCnt = dim(Data)[2]-1;

Time = as.numeric(Data[1,]);
Time = Time[-1];

Status = as.numeric(Data[2,]);
Status = Status[-1];

Status = replace(Status, Status==2, 0);
Status = replace(Status, Status==3, 0);

SurvObj=Surv(Time, Status);

Results=NULL;
Results$GeneNames = as.character(Data[,1]);
Results$GeneNames = Results$GeneNames [-1:-2];
Results$Median=numeric(GeneCnt);
Results$PValues=numeric(GeneCnt);

Groups=matrix(0,GeneCnt, SampleCnt);

for(i in 1:GeneCnt)
	{
		Results$Median[i] = median(as.numeric(Data[2+i,2:(SampleCnt+1)]));
		
		for(j in 1:SampleCnt)
			if(Data[2+i,1+j]>=Results$Median[i])
				Groups[i, j]=1;
		
		sdf = survdiff(SurvObj ~ Groups[i,]);
		Results$PValues[i] = 1 - pchisq(sdf$chisq, length(sdf$n) - 1);
	}

#estimate FDR
#Results$FDRp=fwer2fdr(Results$PValues, method="conservative")$adjp
Results$FDRp=numeric(GeneCnt);

Sorted=sort(Results$PValues, index=TRUE)

Results$PValues=Results$PValues[Sorted$ix]
Results$GeneNames=Results$GeneNames[Sorted$ix]
Results$Median=Results$Median[Sorted$ix]

Groups=Groups[Sorted$ix,]

for(i in GeneCnt:1)
	{
		if(i==GeneCnt)
			{
				Results$FDRp[i]=Results$PValues[i]
			}
		else
			{
				Tmp=Results$PValues[i]*GeneCnt/i;
				
				if(Tmp<Results$FDRp[i+1])
					Results$FDRp[i]=Tmp
				else
					Results$FDRp[i]=Results$FDRp[i+1]
			}
	}


#write text results
write.table(Results, file="results.txt", sep="\t", row.names=FALSE)


#output pdf-s for significant genes
system("rm -rf KM-Plots; mkdir KM-Plots");

setwd("KM-Plots/");

#output pdf
for(i in 1:GeneCnt)
	if(Results$PValues[i]<0.05)
		{
			pdf(paste(Results$GeneNames[i],".pdf", sep=""));
			plot(survfit(SurvObj ~ Groups[i,]), col=c("green", "red"), legend.text=c("Under-expressed", "Over-expressed"))
			title(paste(Results$GeneNames[i]," (p=", round(Results$PValues[i],5), " / FDR-corrected p=", round(Results$FDRp[i],5), ")", sep=""))
			dev.off();
		}

setwd("..");