表达矩阵只需要tumor数据，不要normal，将其去掉，新表达矩阵数据命名为exprSet；
临床信息需要进一步整理，成为生存分析需要的格式，新临床信息数据命名为meta。
由于不同癌症的临床信息表格列名可能不同，这里的代码需要根据实际情况修改。
rm(list=ls())proj="TCGA-KIRC"load(paste0(proj,".Rdata"))library(stringr)
1.整理表达矩阵
不需要正常样本；使用logCPM或logTPM数据
exprSet=log2(edgeR::cpm(exp[,Group=='tumor'])+1)##可以仿照这个将RNA_seq测序的count数据转换成cpm数据，即表达矩阵，这个矩阵可用来画热图ncol(exprSet)
因前面的差异分析过滤标准有宽有严，保险起见，这里可以再次进行基因过滤，至少要在50%的样本里表达量大于0。
k=apply(exprSet,1,function(x){sum(x>0)>0.5*ncol(exprSet)});table(k)#对行进行计算，每行中至少有一半的样本的值大于0exprSet=exprSet[k,]nrow(exprSet)
1.整理生存信息和临床信息
xena将生存信息和临床信息分开了。后续构建模型需要纳入一些临床信息，所以要合并到一起。
1.实现表达矩阵与临床信息的匹配
有的病人会有两个或两个以上的肿瘤样本，就有重复。两种可行的办法：
（1）以病人为中心，对表达矩阵的列按照病人ID去重复，每个病人只保留一个样本。
exprSet=exprSet[,sort(colnames(exprSet))]k=!duplicated(str_sub(colnames(exprSet),1,12));table(k)exprSet=exprSet[,k]
#调整meta行名与exprSet列名一一对应s=intersect(rownames(meta),colnames(exprSet))exprSet=exprSet[,s]meta=meta[s,]identical(rownames(meta),colnames(exprSet))
1.整理生存分析的输入数据
table(meta$event)range(meta$time)meta$time=meta$time/30range(meta$time)
抹除stage里的重复信息
head(meta$stage)meta$stage=meta$stage%>%str_remove("stage")%>%str_to_upper()table(meta$stage,useNA="always")#不需要ABC可以去掉，需要的话就保留，不运行下面这句meta$stage=str_remove(meta$stage,"A|B|C")head(meta)save(meta,exprSet,proj,file=paste0(proj,"_sur_model.Rdata"))
二.生存分析
1.准备输入数据
rm(list=ls())proj="TCGA-KIRC"load(paste0(proj,"_sur_model.Rdata"))ls()exprSet[1:4,1:4]meta[1:4,1:4]
2.KM-plot
简单版本和进阶版本
library(survival)library(survminer)sfit<-survfit(Surv(time,event)~gender,data=meta)ggsurvplot(sfit,pval=TRUE)ggsurvplot(sfit,palette="jco",risk.table=TRUE,pval=TRUE,conf.int=TRUE)
连续型信息怎么作KM分析？例如年龄,基因？
连续型数据的离散化
年龄
group=ifelse(meta$age>median(meta$age,na.rm=T),"older","younger")table(group)sfit=survfit(Surv(time,event)~group,data=meta)ggsurvplot(sfit,pval=TRUE,data=meta,risk.table=TRUE)
基因
g=rownames(exprSet)[1];gmeta$gene=ifelse(exprSet[g,]>median(exprSet[g,]),'high','low')sfit=survfit(Surv(time,event)~gene,data=meta)ggsurvplot(sfit,pval=TRUE,data=meta,risk.table=TRUE)
3.log-ranktest
KM的p值是log-ranktest得出的，可以批量操作
logrankfile=paste0(proj,"_log_rank_p.Rdata")if(!file.exists(logrankfile)){log_rank_p<-apply(exprSet,1,function(gene){meta$group=ifelse(gene>median(gene),'high','low')data.survdiff=survdiff(Surv(time,event)~group,data=meta)p.val=1-pchisq(data.survdiff$chisq,length(data.survdiff$n)-1)return(p.val)})log_rank_p=sort(log_rank_p)save(log_rank_p,file=logrankfile)}load(logrankfile)table(log_rank_p<0.01)table(log_rank_p<0.05)
1.批量单因素cox
5.lasso回归
rm(list=ls())proj="TCGA-KIRC"load(paste0(proj,"_sur_model.Rdata"))ls()exprSet[1:4,1:4]meta[1:4,1:4]load(paste0(proj,"_logrank_cox_gene.Rdata"))exprSet=exprSet[cox,]
2.构建lasso回归模型
输入数据是表达矩阵(仅含tumor样本)和每个病人对应的生死（顺序必须一致）。
x=t(exprSet)#x行名为样本，列名为基因y=meta$eventlibrary(glmnet)
2.1挑选合适的λ值
Lambda是构建模型的重要参数。他的大小关系着模型选择的基因个数
#调优参数set.seed(1006)#选取不同的数，画出来的效果不同cv_fit<-cv.glmnet(x=x,y=y)plot(cv_fit)#系数图fit<-glmnet(x=x,y=y)plot(fit,xvar="lambda")
两条虚线分别指示了两个特殊的λ值,一个是lambda.min,一个是lambda.1se,这两个值之间的lambda都认为是合适的。lambda.1se构建的模型最简单，即使用的基因数量少，而lambda.min则准确率更高一点，使用的基因数量更多一点。
2.2用这两个λ值重新建模
model_lasso_min<-glmnet(x=x,y=y,lambda=cv_fit$lambda.min)model_lasso_1se<-glmnet(x=x,y=y,lambda=cv_fit$lambda.1se)
选中的基因与系数存放于模型的子集beta中，用到的基因有一个s0值，没用的基因只记录了“.”，所以可以用下面代码挑出用到的基因。
head(model_lasso_min$beta,20)choose_gene_min=rownames(model_lasso_min$beta)[as.numeric(model_lasso_min$beta)!=0]choose_gene_1se=rownames(model_lasso_1se$beta)[as.numeric(model_lasso_1se$beta)!=0]length(choose_gene_min)length(choose_gene_1se)save(choose_gene_min,file=paste0(proj,"_lasso_choose_gene_min.Rdata"))save(choose_gene_1se,file=paste0(proj,"_lasso_choose_gene_1se.Rdata"))
3.模型预测和评估
newx参数是预测对象。输出结果lasso.prob是一个矩阵，第一列是min的预测结果，第二列是1se的预测结果，预测结果是概率，或者说百分比，不是绝对的0和1。
将每个样本的生死和预测结果放在一起，直接cbind即可。
lasso.prob<-predict(cv_fit,newx=x,s=c(cv_fit$lambda.min,cv_fit$lambda.1se))re=cbind(y,lasso.prob)head(re)re=as.data.frame(re)colnames(re)=c('event','prob_min','prob_1se')re$event=as.factor(re$event)
ROC曲线
library(pROC)library(ggplot2)m<-roc(meta$event,re$prob_min)g<-ggroc(m,legacy.axes=T,size=1,color="#2fa1dd")auc(m)#Areaunderthecurve:0.9953g+theme_minimal()+geom_segment(aes(x=0,xend=1,y=0,yend=1),colour="grey",linetype="dashed")+annotate("text",x=.75,y=.25,label=paste("AUCofmin=",format(round(as.numeric(auc(m)),2),nsmall=2)),color="#2fa1dd")
计算AUC取值范围在0.5-1之间，越接近于1越好。可以根据预测结果绘制ROC曲线。
两个模型的曲线画在一起
m2<-roc(meta$event,re$prob_1se)auc(m2)#Areaunderthecurve:0.9136g<-ggroc(list(min=m,se=m2),legacy.axes=T,size=1)g+theme_minimal()+scale_color_manual(values=c("#2fa1dd","#f87669"))+geom_segment(aes(x=0,xend=1,y=0,yend=1),colour="grey",linetype="dashed")+annotate("text",x=.75,y=.25,label=paste("AUCofmin=",format(round(as.numeric(auc(m)),2),nsmall=2)),color="#2fa1dd")+annotate("text",x=.75,y=.15,label=paste("AUCof1se=",format(round(as.numeric(auc(m2)),2),nsmall=2)),color="#f87669")
5.切割数据构建模型并预测
5.1切割数据
用R包caret切割数据，生成的结果是一组代表列数的数字，用这些数字来给表达矩阵和meta取子集即可。
library(caret)set.seed(12345679)sam<-createDataPartition(meta$event,p=.5,list=FALSE)head(sam)
可查看两组一些临床参数切割比例
train<-exprSet[,sam]test<-exprSet[,-sam]train_meta<-meta[sam,]test_meta<-meta[-sam,]prop.table(table(train_meta$stage))prop.table(table(test_meta$stage))prop.table(table(test_meta$race))prop.table(table(train_meta$race))
5.2切割后的train数据集建模
和上面的建模方法一样。
#计算lambdax=t(train)y=train_meta$eventcv_fit<-cv.glmnet(x=x,y=y)plot(cv_fit)
#构建模型model_lasso_min<-glmnet(x=x,y=y,lambda=cv_fit$lambda.min)model_lasso_1se<-glmnet(x=x,y=y,lambda=cv_fit$lambda.1se)#挑出基因head(model_lasso_min$beta)choose_gene_min=rownames(model_lasso_min$beta)[as.numeric(model_lasso_min$beta)!=0]choose_gene_1se=rownames(model_lasso_1se$beta)[as.numeric(model_lasso_1se$beta)!=0]length(choose_gene_min)length(choose_gene_1se)
4.模型预测
用训练集构建模型，预测测试集的生死，注意newx参数变了。
lasso.prob<-predict(cv_fit,newx=t(test),s=c(cv_fit$lambda.min,cv_fit$lambda.1se))re=cbind(event=test_meta$event,lasso.prob)re=as.data.frame(re)colnames(re)=c('event','prob_min','prob_1se')re$event=as.factor(re$event)head(re)
再画ROC曲线
library(pROC)library(ggplot2)m<-roc(test_meta$event,re$prob_min)g<-ggroc(m,legacy.axes=T,size=1,color="#2fa1dd")auc(m)#Areaunderthecurve:0.7752g+theme_minimal()+geom_segment(aes(x=0,xend=1,y=0,yend=1),colour="grey",linetype="dashed")+annotate("text",x=.75,y=.25,label=paste("AUCofmin=",format(round(as.numeric(auc(m)),2),nsmall=2)),color="#2fa1dd")
m2<-roc(test_meta$event,re$prob_1se)auc(m2)#Areaunderthecurve:0.7426g<-ggroc(list(min=m,se=m2),legacy.axes=T,size=1)g+theme_minimal()+scale_color_manual(values=c("#2fa1dd","#f87669"))+geom_segment(aes(x=0,xend=1,y=0,yend=1),colour="grey",linetype="dashed")+annotate("text",x=.75,y=.25,label=paste("AUCofmin=",format(round(as.numeric(auc(m)),2),nsmall=2)),color="#2fa1dd")+annotate("text",x=.75,y=.15,label=paste("AUCof1se=",format(round(as.numeric(auc(m2)),2),nsmall=2)),color="#f87669")
6.cox-forest
rm(list=ls())proj="TCGA-KIRC"if(!require(My.stepwise))install.packages("My.stepwise")load(paste0(proj,"_sur_model.Rdata"))load(paste0(proj,"_lasso_choose_gene_1se.Rdata"))g=choose_gene_1se
2.构建coxph模型
将用于建模的基因（例如lasso回归选中的基因）从表达矩阵中取出来，，可作为列添加在meta表噶的后面,组成的数据框赋值给dat。
library(stringr)e=t(exprSet[g,])colnames(e)=str_replace_all(colnames(e),"-","_")dat=cbind(meta,e)dat$gender=as.numeric(factor(dat$gender))dat$stage=as.numeric(factor(dat$stage))colnames(dat)
逐步回归法构建最优模型
输出结果行数太多，所以我注释掉了
library(survival)library(survminer)#不能允许缺失值dat2=na.omit(dat)library(My.stepwise)vl<-colnames(dat)[c(5:ncol(dat))]#My.stepwise.coxph(Time="time",#Status="event",#variable.list=vl,#data=dat2)
使用输出结果里的最后一个模型
model=coxph(formula=Surv(time,event)~stage+age+AL357140.2+C1DP1+HCCAT5+AC131097.2+LINC01522+AC011497.2+PROX1+AC021171.1+INAFM2+GREB1L+CCL22+SLAMF9+LINC01675+AP001893.3+AC092296.1+ZNF320+MZT1P2+CDC42BPG+AL157832.1+AC040934.1+AC018659.8+CHI3L2,data=dat2)
3.模型可视化-森林图
ggforest(model,data=dat2)
fp<-predict(model,newdata=dat2)library(Hmisc)options(scipen=200)with(dat2,rcorr.cens(fp,Surv(time,event)))
C-index用于计算生存分析中的COX模型预测值与真实之间的区分度（discrimination），也称为Harrell’sconcordanceindex。C-index在0.5-1之间。0.5为完全不一致,说明该模型没有预测作用,1为完全一致,说明该模型预测结果与实际完全一致。
library(caret)set.seed(12345679)sam<-createDataPartition(meta$event,p=.5,list=FALSE)train<-exprSet[,sam]test<-exprSet[,-sam]train_meta<-meta[sam,]test_meta<-meta[-sam,]
e=t(train[g,])colnames(e)=str_replace_all(colnames(e),"-","_")dat=cbind(train_meta,e)dat$gender=as.numeric(factor(dat$gender))dat$stage=as.numeric(factor(dat$stage))colnames(dat)library(My.stepwise)dat2=na.omit(dat)vl<-colnames(dat2)[c(5:ncol(dat2))]#My.stepwise.coxph(Time="time",#Status="event",#variable.list=vl,#data=dat2)model=coxph(formula=Surv(time,event)~stage+AC092651.1+MZT1P2+NOC2LP2+CCL22+AC021171.1+INAFM2+LINC01522+AC018630.2+STK19B+ZNF320+GREB1L+NARF+SEMA3A+COL18A1_AS1+HCCAT5+C1DP1+AF230666.2+LRFN1+TGM3+AC092296.1+CDC42BPG+RHNO1+AC107982.3+AL157832.1+AC002070.1,data=dat2)
5.3模型可视化
5.4用切割后的数据test数据集验证模型
e=t(test[g,])colnames(e)=str_replace_all(colnames(e),"-","_")test_dat=cbind(test_meta,e)test_dat$gender=as.numeric(factor(test_dat$gender))test_dat$stage=as.numeric(factor(test_dat$stage))fp<-predict(model,newdata=test_dat)library(Hmisc)with(test_dat,rcorr.cens(fp,Surv(time,event)))
转载自《生信技能树》，如有侵权，请联系删除。
在填写之前，请确认您项目的核心团队至少有一名成员是清华的在校生、校友及教师