--- a
+++ b/km_plot.R
@@ -0,0 +1,86 @@
+library("ipred")
+library("survival")
+library("survivalROC")
+library("glmnet")
+library('kernlab')
+library('caret')
+setwd('D:/pan-cancer/brca')
+par(mfrow = c(3, 4))
+########################LIHC#########################
+kx2=read.csv('brca_cox2.csv',row.names= 1)
+ind=5
+results_fin<-matrix(nrow = ind, ncol=6)
+j=6
+set.seed(j)
+k=10
+folds <- createFolds(as.data.frame(t(kx2)),k)
+results<-matrix(nrow = k, ncol=4)
+eva<-matrix(nrow = nrow(kx2), ncol=k)
+count=0
+for(i in 1:k){
+  testset <- kx2[folds[[i]],] 
+  trainset <- kx2[-folds[[i]],] 
+  x=trainset[,-c(1:2)]
+  x=as.matrix(x)
+  tc=trainset$time
+  tc[tc==0]=0.001
+  y=Surv(tc,trainset$status)
+  cv.fit<-cv.glmnet(x,y,family="cox",maxit=10000,alpha=0,nfold=5)
+  fit<-glmnet(x,y,family="cox",alpha=0)
+  tt=predict(fit,x,s=cv.fit$lambda.min)
+  x2=testset[,-c(1:2)]
+  x2=as.matrix(x2)
+  tc2=testset$time
+  tc2[tc2==0]=0.001
+  y2=Surv(tc2,testset$status)
+  tt2=predict(fit,x2,s=cv.fit$lambda.min)
+  ci_k2=survConcordance(formula = y2~ tt2)
+  
+  
+  results[i,1]<-ci_k2$concordance
+  eva[(count+1):(count+nrow(tt2)),1]<-tc2
+  eva[(count+1):(count+nrow(tt2)),2]<-testset$status
+  eva[(count+1):(count+nrow(tt2)),3]=tt2
+  count=count+nrow(tt2)
+  
+}
+
+tc3=eva[,1]/365
+st3=eva[,2]
+rk3=eva[,3]
+rk4=eva[,4]
+y3=Surv(tc3,st3)
+ci_k3=survConcordance(formula = y3~ rk3)
+ci_k3$concordance
+mm=median(rk3)
+for (p in 1:length(rk3)){
+  if (rk3[p]>mm){
+    rk3[p]=1 
+  } else  
+  {rk3[p]=-1}
+}
+q2=survdiff(y3~rk3)
+p.val <- 1 - pchisq(q2$chisq, length(q2$n) - 1)
+p.val
+sd11=survfit(y3~rk3)
+
+
+
+
+################################################
+kt2=kx2[,-(1:2)]
+cl <- kmeans(kt2, 2)
+riskk=cl$cluster
+yy=Surv((kx2$time)/365,kx2$status)
+sp=survdiff(yy~riskk)
+p.valk2 <- 1 - pchisq(sp$chisq, length(sp$n) - 1)
+p.valk2
+riskk[riskk==2]=0
+#riskk_bl=riskk
+sd12=survfit(yy~riskk_bl)
+
+
+plot(sd11,conf.int=FALSE,mark.time=T,col=c("green","red"),lty=1,main="BRCA",
+     cex.main=1.5,xlab="Year",ylab = "Survival Probablity",cex.lab=1.4,cex.axis=1.2,lwd=2)
+lines(sd12,mark.time=T,col=c("grey","grey"),lty=2,lwd=2)
+legend('topright',c('DCAP-low risk','DCAP-high risk'),lty=c(1,1),col=c(3,2),cex=1.2)
\ No newline at end of file