突变位点生存分析

options(download.file.method='libcurl')
options(url.method='libcurl')
BiocManager::install("maftools")
rm(list = ls())
library(stringr)
options(stringsAsFactors = F)
maf<-"C:\\Users\\12162\\Desktop\\微信公众号\\学徒作业\\test_1\\maf\\gdc_download_20200427_034400.797422.tar\\gdc_download_20200427_034400\\6c93f518-1956-4435-9806-37185266d248\\TCGA.BRCA.varscan.6c93f518-1956-4435-9806-37185266d248.DR-10.0.somatic.maf.gz"
laml<-read.maf(maf=maf)
data<-laml@data
#样本水平
getSampleSummary(laml)
#基因水平
getGeneSummary(laml)
#将summary结果输出
write.mafSummary(maf=laml,basename='laml')
#对summary结果进行可视化
pdf("sum.pdf")
plotmafSummary(maf=laml)
dev.off()

#oncoplot可视化
pdf("top10.pdf")
oncoplot(maf=laml,top=10)
dev.off()

#提取含有TP53和KMT2C的样本
tumor<-rbind(data[data$Hugo_Symbol=="TP53",],data[data$Hugo_Symbol=="KMT2C",])
#load("gdc.Rdata")
#整理临床信息----
library(XML)
result <- xmlParse("./clinical/00049989-fa21-48fb-8dda-710c0dd5932e/nationwidechildrens.org_clinical.TCGA-A2-A0CT.xml")
rootnode <- xmlRoot(result)
rootsize <- xmlSize(rootnode)
print(rootnode[2])
xmldataframe <- xmlToDataFrame(rootnode[2])
head(t(xmlToDataFrame(rootnode[2])))

xmls = dir("clinical/",pattern = "*.xml$",recursive = T)
td = function(x){
  result <- xmlParse(file.path("clinical/",x))
  rootnode <- xmlRoot(result)
  xmldataframe <- xmlToDataFrame(rootnode[2])
  return(t(xmldataframe))
}
cl = lapply(xmls,td)
cl_df <- t(do.call(cbind,cl))
tmp = data.frame(colnames(cl_df))
clinical = cl_df[,c(
  'bcr_patient_barcode',
  'vital_status',
  'days_to_death',
  'days_to_last_followup',
  'race_list',
  'days_to_birth',
  'gender' ,
  'stage_event'
)]
clinical = data.frame(clinical)
rownames(clinical) <- NULL
rownames(clinical) <- stringr::str_to_upper(clinical$bcr_patient_barcode)
dim(clinical)
meta = clinical
#简化列名
colnames(meta)=c('ID','event','death','last_followup','race','age','gender','stage')
#根据选取基因进行分组----
tmp<-tumor$Tumor_Sample_Barcode
tmp<-substr(tmp,1,12)
tmp<-match(tmp,as.character(meta$ID))
group_list<-c(1:nrow(meta))
group_list[sort(as.numeric(na.omit(tmp)))]<-"TP53/KMT2C+"
group_list<-ifelse(group_list=="TP53/KMT2C+","TP53/KMT2C+","TP53/KMT2C-")
table(group_list)

#整理生存分析的输入数据----
#1.由随访时间和死亡时间计算生存时间(月)
is.empty.chr = function(x){
  ifelse(stringr::str_length(x)==0,T,F)
}
is.empty.chr(meta[2,3])
meta[,3][is.empty.chr(meta[,3])]=0
meta[,4][is.empty.chr(meta[,4])]=0
meta$time=(as.numeric(meta[,3])+as.numeric(meta[,4]))/30
#2.根据生死定义event，活着是0，死的是1
meta$event=ifelse(meta$event=='Alive',0,1)
#生存分析----
library(survival)
library(survminer)
sfit <- survfit(Surv(time, event)~group_list, data=meta)
ggsurvplot(sfit, conf.int=F, pval=TRUE,surv.median.line = "hv")
pdf("result.pdf")
ggsurvplot(sfit,palette = c("#E7B800", "#2E9FDF"),
           risk.table =TRUE,pval =TRUE,
           conf.int =TRUE,xlab ="Time in months", 
           ggtheme =theme_light())
dev.off()