突变位点生存分析
好奇怪，最近大家完成学徒作业的积极性很差，是习题太难了吗？一个简单突变位点做生存分析居然拖了一两个月才有人提交笔记！前面的题目见：学徒作业-两个基因突变联合看生存效应 （2020-04-26出题），下面看其中一个学徒的答案哦，同时也欢迎大家继续提交笔记给我哈，有机会认识我！加油哈，广大粉丝们
1.本次选用BRCA的maf数据和临床数据，主要使用其中的varscan数据
2.使用R包maftools读取maf文件，并可视化top10突变基因
3.选取两基因对BRCA临床样本进行分组
所选取两基因都未发生突变的样本为一组
剩余样本为一组
4.使用logrank进行生存分析
1.读取maf文件并对数据进行可视化
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()
2.挑选top10基因，进行可视化
#oncoplot可视化
pdf("top10.pdf")
oncoplot(maf=laml,top=10)
dev.off()
3.整理临床数据，挑选目的基因
#提取含有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)
4.针对基因突变与否，进行生存分析
#整理生存分析的输入数据----
#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()