TCGA癌症数据挖掘之预后模型建立和评价 by 生信技能树

众所周知，我们【生信技能树】团队这些年一直在身体力行的推广马拉松学习理念，目前唯一可行的生物信息学入门策略！是60个小时的线上全网直播视频课程，连续4个星期，每个星期5天，每天的晚上8~11点的3个小时互动授课（周三、周日休息）

与十万人一起学生信，你值得拥有下面的学习班：

数据挖掘（GEO,TCGA,单细胞）2022年5月场，快速了解一些生物信息学应用图表
生信入门课-2022年5月场，你的生物信息学第一课

同期也安排了助教和实习生一起协助整理讲师的直播授课时间线和授课内容相关笔记，提炼内容梗概方便学员们查漏补缺，精准复现！

下面是学员的（数据挖掘）直播配套笔记

一. “生存分析前的数据整理”

1.读入数据

表达矩阵只需要tumor数据，不要normal，将其去掉，新表达矩阵数据命名为exprSet；

临床信息需要进一步整理，成为生存分析需要的格式，新临床信息数据命名为meta。

由于不同癌症的临床信息表格列名可能不同，这里的代码需要根据实际情况修改。

rm(list=ls())
proj = "TCGA-KIRC"
load(paste0(proj,".Rdata"))
library(stringr)

2.整理表达矩阵

不需要正常样本；使用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) # 对行进行计算，每行中至少有一半的样本的值大于0
exprSet = exprSet[k,]
nrow(exprSet)

3.整理生存信息和临床信息

xena将生存信息和临床信息分开了。后续构建模型需要纳入一些临床信息，所以要合并到一起。

library(dplyr)
meta = left_join(surv,clinical,by = c("sample"= "submitter_id.samples"))
# 去掉表达矩阵里没有的样本
library(stringr)
k = meta$sample %in% colnames(exprSet);table(k)
meta = meta[k,]

# 去掉生存信息不全或者生存时间小于30天的样本，样本纳排标准不唯一，且差别很大
k1 = meta$OS.time >= 30;table(k1)
k2 = !(is.na(meta$OS.time)|is.na(meta$OS));table(k2)
meta = meta[k1&k2,]

# 选择有用的列
tmp = data.frame(colnames(meta))
meta = meta[,c(
  'sample',
  'OS',
  'OS.time',
  'race.demographic',
  'age_at_initial_pathologic_diagnosis',
  'gender.demographic' ,
  'tumor_stage.diagnoses'
)]

dim(meta)
rownames(meta) <- meta$sample
meta[1:4,1:4]

#简化meta的列名
colnames(meta)=c('ID','event','time','race','age','gender','stage')

#空着的值、not reported改为NA
meta[meta==""|meta=="not reported"]=NA

3.实现表达矩阵与临床信息的匹配

有的病人会有两个或两个以上的肿瘤样本，就有重复。两种可行的办法：

（1）以病人为中心，对表达矩阵的列按照病人ID去重复，每个病人只保留一个样本。

exprSet = exprSet[,sort(colnames(exprSet))]
k = !duplicated(str_sub(colnames(exprSet),1,12));table(k)
exprSet = exprSet[,k]

（2）以样本为中心，如果每个病人有多个样本则全部保留。(删掉上面这一段代码即可)

#调整meta行名与exprSet列名一一对应
s = intersect(rownames(meta),colnames(exprSet))
exprSet = exprSet[,s]
meta = meta[s,]
identical(rownames(meta),colnames(exprSet))

4. 整理生存分析的输入数据

生存分析的输入数据里，要求结局事件必须用0和1表示，0表示活着，1表示死了;生存时间的单位（月）;

table(meta$event)
range(meta$time)
meta$time = meta$time/30
range(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];g
meta$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-rank test

KM的p值是log-rank test得出的，可以批量操作

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)

4.批量单因素cox

coxfile = paste0(proj,"_cox.Rdata")
if(!file.exists(coxfile)){
  cox_results <-apply(exprSet , 1 , function(gene){
  meta$gene = gene
  #可直接使用连续型变量
  m = coxph(Surv(time, event) ~ gene, data =  meta)
  #也可使用二分类变量
  #meta$group=ifelse(gene>median(gene),'high','low') 
  #meta$group = factor(meta$group,levels = c("low","high"))
  #m=coxph(Surv(time, event) ~ group, data =  meta)
  
  beta <- coef(m)
  se <- sqrt(diag(vcov(m)))
  HR <- exp(beta)
  HRse <- HR * se
  
  #summary(m)
  tmp <- round(cbind(coef = beta, 
                     se = se, z = beta/se, 
                     p = 1 - pchisq((beta/se)^2, 1),
                     HR = HR, HRse = HRse,
                     HRz = (HR - 1) / HRse, 
                     HRp = 1 - pchisq(((HR - 1)/HRse)^2, 1),
                     HRCILL = exp(beta - qnorm(.975, 0, 1) * se),
                     HRCIUL = exp(beta + qnorm(.975, 0, 1) * se)), 3)
  
  return(tmp['gene',]) 
  #return(tmp['grouphigh',])#二分类变量
})
  cox_results=as.data.frame(t(cox_results))
  save(cox_results,file = coxfile)
}
load(coxfile)
table(cox_results$p<0.01)
table(cox_results$p<0.05)

lr = names(log_rank_p)[log_rank_p<0.01];length(lr)
cox = rownames(cox_results)[cox_results$p<0.01];length(cox)
length(intersect(lr,cox))
save(lr,cox,file = paste0(proj,"_logrank_cox_gene.Rdata"))

5.lasso回归

1.准备输入数据

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$event
library(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)  # Area under the curve: 0.9953

g + 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("AUC of min = ",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) # Area under the curve: 0.9136
g <- 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("AUC of min = ",format(round(as.numeric(auc(m)),2),nsmall = 2)),color = "#2fa1dd")+
  annotate("text",x = .75, y = .15,
           label = paste("AUC of 1se = ",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数据集建模

和上面的建模方法一样。

#计算lambda
x = t(train)
y = train_meta$event
cv_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) #Area under the curve: 0.7752

g + 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("AUC of min = ",format(round(as.numeric(auc(m)),2),nsmall = 2)),color = "#2fa1dd")

计算AUC取值范围在0.5-1之间，越接近于1越好。可以根据预测结果绘制ROC曲线。

两个模型的曲线画在一起

m2 <- roc(test_meta$event, re$prob_1se)
auc(m2)  # Area under the curve: 0.7426
g <- 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("AUC of min = ",format(round(as.numeric(auc(m)),2),nsmall = 2)),color = "#2fa1dd")+
  annotate("text",x = .75, y = .15,
           label = paste("AUC of 1se = ",format(round(as.numeric(auc(m2)),2),nsmall = 2)),color = "#f87669")

6.cox-forest

1.准备输入数据

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)

4.模型预测

fp <- predict(model,newdata = dat2)
library(Hmisc)
options(scipen=200)
with(dat2,rcorr.cens(fp,Surv(time, event)))

C-index用于计算生存分析中的COX模型预测值与真实之间的区分度（discrimination），也称为Harrell’s concordanceindex。C-index在0.5-1之间。0.5为完全不一致,说明该模型没有预测作用,1为完全一致,说明该模型预测结果与实际完全一致。

5.切割数据构建模型并预测

5.1 切割数据

用R包caret切割数据，生成的结果是一组代表列数的数字，用这些数字来给表达矩阵和meta取子集即可。

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,]

5.2 切割后的train数据集建模

和上面的建模方法一样。

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 模型可视化

ggforest(model, data =dat2)

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)))

ixxmu / mp_duty

TCGA癌症数据挖掘之预后模型建立和评价 #2112

TCGA癌症数据挖掘之预后模型建立和评价 by 生信技能树

1.读入数据

2.整理表达矩阵

3.整理生存信息和临床信息

3.实现表达矩阵与临床信息的匹配

4. 整理生存分析的输入数据

二.生存分析

1.准备输入数据

2.KM-plot

3.log-rank test

4.批量单因素cox

5.lasso回归

1.准备输入数据

2.构建lasso回归模型

2.1挑选合适的λ值

2.2 用这两个λ值重新建模

3.模型预测和评估

ROC曲线

5.切割数据构建模型并预测

5.1 切割数据

5.2 切割后的train数据集建模

4.模型预测

6.cox-forest

1.准备输入数据

2.构建coxph模型

逐步回归法构建最优模型

3.模型可视化-森林图

4.模型预测

5.切割数据构建模型并预测

5.1 切割数据

5.2 切割后的train数据集建模

5.3 模型可视化

5.4 用切割后的数据test数据集验证模型