#!/usr/bin/Rscript

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

#  fullannot  # colnames  rownames
#  TSS1500Ind # TSS1500Ind$SID  names(TSS1500Ind$SID) TSS1500Ind$SID[[1]]  TSS1500Ind$SID[["ZYX"]]
#  TSS200Ind  
#  UTR5Ind    
#  EXON1Ind  
#  GENEBODYInd
#  UTR3Ind
   
#  ISLANDInd # ISLANDInd$SID
#$`chrY:9930123-9930391`
#[1] "cg14210405"

#  NSHOREInd  
#  SSHOREInd  
#  NSHELFInd  
#  SSHELFInd

sf_df <- matrix(nrow=2000, ncol=32)
sf_df <- as.data.frame(sf_df)
colnames(sf_df) <- c("BLCA_N" , "BLCA_T",
"BRCA_N" , "BRCA_T",
"CHOL_N" , "CHOL_T",
"COAD_N" , "COAD_T",
"ESCA_N" , "ESCA_T",
"HNSC_N" , "HNSC_T",
"KIRC_N" , "KIRC_T",
"KIRP_N" , "KIRP_T",
"LIHC_N" , "LIHC_T",
"LUAD_N" , "LUAD_T",
"LUSC_N" , "LUSC_T",
"PRAD_N" , "PRAD_T",
"READ_N" , "READ_T",
"STES_N" , "STES_T",
"THCA_N" , "THCA_T",
"UCEC_N" , "UCEC_T")

tumors <- scan("tumor_list",what="")
for(tumor in tumors) {

#sf file names
   sf_site <- paste(tumor,".site.ratio",sep="")
   sf_ele  <- paste(tumor,".ele.ratio",sep="")
   sf_island <-paste(tumor,".islands",sep="")
   
#read in files
   print(paste("dealing with",tumor,sep=" "))
   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)
   exp_file <- paste("/home/zhanghl/resource/raw_counts/",tumor,".paired.expression",sep="")
   exp <- read.table(exp_file,header=TRUE,row.names=1,nrow=5,sep="\t",stringsAsFactors=FALSE)
   
#common samples exp methy
   methy_sub_names <- substr(colnames(methy),1,12)
   exp_sub_names <- substr(colnames(exp),1,12)
   com_sub_names <- unique(intersect(methy_sub_names,exp_sub_names))
   t_names <- paste(com_sub_names,".01",sep="")
   n_names <- paste(com_sub_names,".11",sep="")
   
#idx x and aa probes    
   anno_x_idx <- fullannot[,"CHR"]=="X"
   anno_aa_idx <- fullannot[,"CHR"]!="X" & fullannot[,"CHR"]!="Y" & !is.na(fullannot[,"CHR"])
   anno_x_probes <- rownames(fullannot)[anno_x_idx]
   anno_aa_probes <- rownames(fullannot)[anno_aa_idx]
   common_x_probes <- intersect(anno_x_probes,rownames(methy))
   common_aa_probes<- intersect(anno_aa_probes,rownames(methy))
   
#idx tumor and normal samples
   t_idx <- t_names
   n_idx <- n_names
   x_t_median <- apply(methy[common_x_probes,t_idx],2,median,na.rm=TRUE)
   x_n_median <- apply(methy[common_x_probes,n_idx],2,median,na.rm=TRUE)
   aa_t_median<- apply(methy[common_aa_probes,t_idx],2,median,na.rm=TRUE)
   aa_n_median<- apply(methy[common_aa_probes,n_idx],2,median,na.rm=TRUE)
   methy_ratio_t <- x_t_median/aa_t_median
   methy_ratio_n <- x_n_median/aa_n_median
   
   idx <- which(tumors%in%tumor)
   sf_df[,(idx*2-1)] <- c(methy_ratio_n,rep(NA,(2000-length(methy_ratio_n))))
   sf_df[,idx*2]  <- c(methy_ratio_t,rep(NA,(2000-length(methy_ratio_t))))
   
   }
   write.table(sf_df,"sf_change_per_probe",row.names=FALSE,col.names=TRUE,quote=FALSE,sep="\t")
   
   
   
   setwd("F:/paper1 related files/Altered expression ratios of X chromosome to autosomes in tumors(v3)/version3/fig4.methy_ratio")
   data <- read.table("sf_change_per_probe",header=TRUE,sep="\t",stringsAsFactors = FALSE)
   png("methy_ratio.png")
   boxplot(data,col=c("yellow","green"),las=2 ,ylab="Methylation ratio of X to AA",cex.lab=1.2,cex.axis=1.2,notch = FALSE)
   abline(v=2.5)
   abline(v=4.5)
   abline(v=6.5)
   abline(v=8.5)
   abline(v=10.5)
   abline(v=12.5)
   abline(v=14.5)
   abline(v=16.5)
   abline(v=18.5)
   abline(v=20.5)
   abline(v=22.5)
   abline(v=24.5)
   abline(v=26.5)
   abline(v=28.5)
   abline(v=30.5)
   
   dev.off()
   
   names <- colnames(data)    
   i=1
   while(i<32){
   #test
   result=t.test(data[,i],data[,i+1],paired = TRUE)
   print(paste(names[i],result$p.value,length(na.omit(data[,i])),sep="    "))
   i=i+2
   }