library(survival)
read.delim("THCA.clin.merged.txt",header=T,sep="\t",row.names=1)->clinical

clinical[c("patient.bcr_patient_barcode","patient.days_to_death","patient.vital_status","patient.days_to_last_followup"),]->clinical2


gsub("-",".",as.matrix(clinical2)[1,])->clinical_name
paste(clinical_name,".01",sep="")->clinical_name2
toupper(clinical_name2)->clinical_name2
colnames(clinical2)<-clinical_name2

live<-as.numeric(clinical2[3,]!="alive")
days<-clinical2[4,]
days[which(live==1)]<-clinical2[2,which(live==1)]   

rbind(as.vector(t(live)),as.vector(t(days)))->clinical_last
colnames(clinical_last)<-clinical_name2


read.table("../gdac.broadinstitute.org_THCA.Merge_rnaseqv2__illuminahiseq_rnaseqv2__unc_edu__Level_3__RSEM_genes_normalized__data.Level_3.2016012800.0.0/gdac.broadinstitute.org_THCA.mRNAseq_Preprocess.Level_3.2016012800.0.0/THCA.uncv2.mRNAseq_RSEM_normalized_log2.txt",header=T,row.names=1)->a

a[is.na(a)]=-9
rownames(a)->names
agrep("\\?",names)->ques
unlist(lapply(names,function(x) strsplit(x,"\\|")[[1]][1]))->filter.names
filter.names[ques]<-names[ques]
filter.names[which(filter.names=="SLC35E2")][2]="SLC35E2|9906"
rownames(a)<-filter.names
colnames(a)->col.name
unlist(lapply(col.name,function(x) strsplit(x,"\\.")[[1]][length(strsplit(col.name[1],"\\.")[[1]])]))->type
a[,which(type%in%c("01","02","03","04","05","06","07","08","09"))]->tumor
a[,which(type=="11")]->normal
colnames(tumor)->exp_name

sapply(strsplit(exp_name,"\\."),function(x) paste(x[1:3],collapse="."))-> exp_name2
paste(exp_name2,".01",sep="")-> exp_name3
colnames(tumor)<-exp_name3

intersect(colnames(clinical_last),colnames(tumor))-> inter_BRAF_mut_exp
tumor["NRARP",inter_BRAF_mut_exp]-> PRC2_exp_sum

PRC2_exp_sum-> PRC2_exp_sum2

rbind(clinical_last[, inter_BRAF_mut_exp],unlist(PRC2_exp_sum[inter_BRAF_mut_exp]))->matrix_BRAF


matrix_BRAF [2,which(matrix_BRAF[2,]<0)]=0
rownames(matrix_BRAF)<-c("event","days","PRC2")
matrix(0,nrow(matrix_BRAF),ncol(matrix_BRAF))-> matrix_BRAF2
as.numeric(matrix_BRAF[1,])-> matrix_BRAF2[1,]
as.numeric(matrix_BRAF[2,])-> matrix_BRAF2[2,]
as.numeric(matrix_BRAF[3,])-> matrix_BRAF2[3,]

rownames(matrix_BRAF2)<-rownames(matrix_BRAF)
summary(coxph(formula = Surv(days,event) ~PRC2, as.data.frame(t(matrix_BRAF2))))

fit<-survfit(Surv(days,event) ~1, as.data.frame(t(matrix_BRAF2)))

#######

unlist(PRC2_exp_sum)-> PRC2_exp_sum

names(which(PRC2_exp_sum<=quantile(PRC2_exp_sum,seq(0,1,1/3))[2]))->bot_sample

names(which(PRC2_exp_sum>quantile(PRC2_exp_sum,seq(0,1,1/3))[3]))->top_sample

rep(0,length(top_sample)+length(bot_sample))-> PRC2_exp_sum2
names(PRC2_exp_sum2)<-c(top_sample, bot_sample)
PRC2_exp_sum2[top_sample]=1

c(top_sample, bot_sample)->diff_sample
rbind(clinical_last[, diff_sample],unlist(PRC2_exp_sum2[diff_sample]))->matrix_BRAF

matrix_BRAF [2,which(matrix_BRAF[2,]<0)]=0
rownames(matrix_BRAF)<-c("event","days","PRC2")
matrix(0,nrow(matrix_BRAF),ncol(matrix_BRAF))-> matrix_BRAF2
as.numeric(matrix_BRAF[1,])-> matrix_BRAF2[1,]
as.numeric(matrix_BRAF[2,])-> matrix_BRAF2[2,]
as.numeric(matrix_BRAF[3,])-> matrix_BRAF2[3,]

rownames(matrix_BRAF2)<-rownames(matrix_BRAF)
colnames(matrix_BRAF2)<-colnames(matrix_BRAF)
summary(coxph(formula = Surv(days,event) ~PRC2, as.data.frame(t(matrix_BRAF2))))
plot(survfit(formula = Surv(days, event) ~ PRC2,data = data.frame(t(matrix_BRAF2))),col=c("blue","red"),main="")
legend("topright",legend=c("NRARP_Low1/3%","NRARP_Top1/3"),col=c("blue","red"),lty=1)




boxplot(unlist(log(exp[gene,7:24])),unlist(log(exp[gene,25:28])))