下面是學員的（ 資料挖掘 ）直播配套筆記
一. “生存分析前的資料整理”
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)))
複制