#!/usr/bin/Rscript

#calculate how many down and up -regulated genes are hyper and hypo at promoter
#dealing with down-regulated genes

load("/home/zhanghl/workshop/paper1/methylation/anno/fullannotInd.rda")

#  fullannot  # colnames  rownames
#  TSS1500Ind # TSS1500Ind$SID  names(TSS1500Ind$SID) TSS1500Ind$SID[[1]]  TSS1500Ind$SID[["ZYX"]]

com_down_genes <- scan("/home/zhanghl/workshop/paper1/expression/analysis_common_up_down_genes/survival_gene_list_down",what="")
tumors <- c("BRCA","KIRC","PRAD","THCA","UCEC")
#readin methy data
for (tumor in tumors){
   spc_file <- paste("/home/zhanghl/workshop/paper1/expression/analysis_common_up_down_genes/",tumor,".sig.down.genes",sep="")    
   spc_genes <- scan(spc_file,what="")
   com_genes <- intersect(spc_genes,com_down_genes)
   com_genes_symbol <- sub("\|[a-zA-Z0-9?]*","",com_genes,fixed=FALSE)
   inter_genes_anno <- intersect(names(TSS1500Ind$SID),com_genes_symbol)
       
   methy_file <- paste("/home/zhanghl/workshop/paper1/methylation/methylaiton_data_v3/",tumor,".paired.methylation2",sep="")
   methy <- read.table(methy_file,header=TRUE,row.names=1,sep="t",stringsAsFactors=FALSE)
   methy_patients <- unique(substr(colnames(methy),1,12))
   t_names <- paste(methy_patients,".01",sep="")
   n_names <- paste(methy_patients,".11",sep="")
   


   sf <- matrix(ncol=3,nrow=length(inter_genes_anno))
   sf <- as.data.frame(sf)
   colnames(sf) <- c("normal","Tumor","hypo_hyper")
   rownames(sf) <- inter_genes_anno
   
   for (i in inter_genes_anno){
       probes <- TSS1500Ind$SID[[i]]
       methy_ele_t_media <- apply(methy[probes,t_names],2,median,na.rm=TRUE)
       methy_ele_n_median <- apply(methy[probes,n_names],2,median,na.rm=TRUE)
       
       sf[i,1] <- median(methy_ele_n_median,na.rm=TRUE)
       sf[i,2] <- median(methy_ele_t_media,na.rm=TRUE)
       
       test <- wilcox.test(methy_ele_t_media,methy_ele_n_median,paired=TRUE,exact = NULL, correct = TRUE)
       sig <- test$p.value
       sf[i,3] <- ifelse(sig>0.05,"NULL",ifelse(sf[i,1]>sf[i,2],"hypo","hyper"))
       
   }    
   write.table(sf,paste(tumor,".down.genes",sep=""),row.names=TRUE,col.names=TRUE,sep="t",quote=FALSE)
   hyper_num <- sum(sf[,3]=="hyper")
   hypo_num  <- sum(sf[,3]=="hypo")
   null_num  <- sum(sf[,3]=="NULL")
   print (paste("hyper_num:",hyper_num,"hypo_num:",hypo_num,"null_null:",null_num,sep=" "))
   }