本篇教程来自Biostars的一篇帖子：https://www.biostars.org/p/153013/，内容涉及TCGA RNASeq及Clinical Data数据整理和生存分析（总生存期和DFS）。
数据下载
原文中的数据来自FireBrowse，我已经下载好：http://58.58.190.9:12345/s/48lUeSCIrIbQiez，下载解压后两个文件夹Clinical和RNA, 解压到自己工程根目录。
临床数据整理
与大部分国内教程不同的是，作者仔细的考虑了生存时间、censoring 、无病生存期（DFS）的区别，并对这些事件进行了分开统计；而国内绝大部分教程都是总生存时间。我们知道，生存分析中复发应该是一个事件，某种意义上，复发≈死亡。DFS一般指术后到复发的时间间隔。
首先读取数据
# 读取临床数据
clinical <- read.table('Clinical/KIRC.merged_only_clinical_clin_format.txt',header=T, row.names=1, sep='\t')
clinical <- as.data.frame(t(clinical))
# 查看RNA-seq中存在的病例
clinical$IDs <- toupper(clinical$patient.bcr_patient_barcode)
sum(clinical$IDs %in% colnames(z_rna))
查看一下clinical数据有哪些列
> colnames(clinical)
[1] "admin.disease_code"
[2] "admin.file_uuid"
[3] "admin.month_of_dcc_upload"
[4] "admin.patient_withdrawal.withdrawn"
[5] "admin.project_code"
[6] "admin.year_of_dcc_upload"
[7] "patient.additional_studies"
[8] "patient.additional_studies.additional_study.disease_code"
[9] "patient.additional_studies.additional_study.project_code"
[10] "patient.age_at_initial_pathologic_diagnosis"
[11] "patient.bcr_patient_barcode"
[12] "patient.bcr_patient_uuid"
...
[515] "patient.year_of_tobacco_smoking_onset"
[516] "IDs"
这是一个巨大的数据框，一个一个指定所选列是不现实的。我们可以使用grep查找：
# 找出有days_to_new_tumor_event_after_initial_treatment这个关键字的列，这些列中记录着复发数据
# grep方法返回一个下标数组：17 265 291 317 343，即这些列含有这个关键字
ind_keep <- grep('days_to_new_tumor_event_after_initial_treatment',colnames(clinical))
# 将复发病人数据筛选出来
new_tum <- as.matrix(clinical[,ind_keep])
dim(new_tum)
# [1] 537 5
new_tum（由于列名太长所以显示不全）
image
首先将肿瘤复发事件整理到一个向量new_tum_collapsed
new_tum_collapsed <- c()
# 遍历new_tum每一行
for (i in 1:dim(new_tum)[1]){
# 如果这一行所有值不全是NA（dim(new_tum)[2]即列数（5））
if (sum(is.na(new_tum[i,])) < dim(new_tum)[2]){
# TCGA数据中days_to_last_followup数据往往有多个，我们应该选择最近的一个，即最大的一个
m <- max(new_tum[i,],na.rm=T)
# 将follow-up最大值加入 new_tum_collapsed 向量
new_tum_collapsed <- c(new_tum_collapsed,m)
} else {
# 否则加入NA，即无肿瘤复发事件
new_tum_collapsed <- c(new_tum_collapsed,'NA')
}
}
对死亡事件和最后随访事件如法炮制：
# do the same to death
ind_keep <- grep('days_to_death',colnames(clinical))
death <- as.matrix(clinical[,ind_keep])
death_collapsed <- c()
for (i in 1:dim(death)[1]){
if ( sum ( is.na(death[i,])) < dim(death)[2]){
m <- max(death[i,],na.rm=T)
death_collapsed <- c(death_collapsed,m)
} else {
death_collapsed <- c(death_collapsed,'NA')
}
}
# and days last follow up here we take the most recent which is the max number
ind_keep <- grep('days_to_last_followup',colnames(clinical))
fl <- as.matrix(clinical[,ind_keep])
fl_collapsed <- c()
for (i in 1:dim(fl)[1]){
if ( sum (is.na(fl[i,])) < dim(fl)[2]){
m <- max(fl[i,],na.rm=T)
fl_collapsed <- c(fl_collapsed,m)
} else {
fl_collapsed <- c(fl_collapsed,'NA')
}
}
然后将所有结果整合到一起，命名为'new_tumor_days', 'death_days', 'followUp_days'
# 将复发、死亡和follow_up数据整合一起
all_clin <- data.frame(new_tum_collapsed,death_collapsed,fl_collapsed)
colnames(all_clin) <- c('new_tumor_days', 'death_days', 'followUp_days')
# 如果某个病例有复发事件，使用new_tumor_days作为生存事件，否则使用followUp_days
all_clin$new_time <- c()
for (i in 1:length(as.numeric(as.character(all_clin$new_tumor_days)))){
all_clin$new_time[i] <- ifelse(
is.na(as.numeric(as.character(all_clin$new_tumor_days))[i]), as.numeric(as.character(all_clin$followUp_days))[i],
as.numeric(as.character(all_clin$new_tumor_days))[i])
}
# 如果某个病例有死亡事件，使用death_days作为生存事件，否则使用followUp_days
all_clin$new_death <- c()
for (i in 1:length(as.numeric(as.character(all_clin$death_days)))){
all_clin$new_death[i] <- ifelse (
is.na(as.numeric(as.character(all_clin$death_days))[i]),
as.numeric(as.character(all_clin$followUp_days[i],
as.numeric(as.character(all_clin$death_days))[i])
}
现在我们需要创建一个矢量来审查数据，这意味着我们需要告诉R哪些病人已经死亡或有新的肿瘤。在这种情况下，如果一个病人有一个“days_to_death”，它将被赋值为2，否则赋值为1（使用1，2要比0，1更好）
即使对于复发病人，我们也使用死亡事件进行censor。因为，虽然病人可能复发然后死亡，但如果病人死了，就没有复发事件。因此更准确的是审查死亡事件而不是复发。
因此，我们统一使用死亡作为审查事件。
# create vector for death censoring
table(clinical$patient.vital_status)
# alive dead
# 372 161
# 将0，1 修改为1，2
all_clin$death_event <- ifelse(clinical$patient.vital_status == 'alive', 1,2)
#finally add row.names to clinical
rownames(all_clin) <- clinical$IDs
all_clin:
image
下面的生存分析，我们要用到的是最后3列数据，如果你想看DFS，那就用new_time作为生存时间，如果你想看总体生存，那就用new_death。
现在是时候整理RNA-seq数据了
RNA Seq 数据整理
library(survival)
library(limma)
# 读取转录组数据
rna <- read.table('RNA/KIRC.rnaseqv2__illuminahiseq_rnaseqv2__unc_edu__Level_3__RSEM_genes_normalized__data.data.txt',nrows=20533, header=T,row.names=1,sep='\t')
# 删除第一行
rna <- rna[-1,]
# 移除低表达数据(>50%的样品中表达量为0)
rem <- function(x){
x <- as.matrix(x)
x <- t(apply(x,1,as.numeric))
r <- as.numeric(apply(x,1,function(i) sum(i == 0)))
remove <- which(r > dim(x)[2]*0.5)
return(remove)
}
remove <- rem(rna)
rna <- rna[-remove,]
# 显示正常和肿瘤样本
table(substr(colnames(rna),14,14))
# get the index of the normal/control samples
n_index <- which(substr(colnames(rna),14,14) == '1')
t_index <- which(substr(colnames(rna),14,14) == '0')
# 使用limma包的voom方法对数据进行标准化
vm <- function(x){
cond <- factor(ifelse(seq(1,dim(x)[2],1) %in% t_index, 1, 0))
d <- model.matrix(~1+cond)
x <- t(apply(x,1,as.numeric))
ex <- voom(x,d,plot=F)
return(ex$E)
}
rna_vm <- vm(rna)
# 只保留barcode的1-12位
colnames(rna_vm) <- gsub('\\.','-',substr(colnames(rna),1,12))
# 查看样本分布
hist(rna_vm)
# 删除rna数据
rm(rna)
与大部分国内教程不同的是，这里作者使用了Z分数（Z-Score）而不是Fold Change（FC）来表示不同样本之间的基因差异表达。因为一旦使用FC，就意味着在整个样本水平上平均表达数据，这样会丢失样本的“特质性(heterogeity)”。z分数可以回答这样一个问题："一个给定分数距离平均数多少个标准差?"。在平均数之上的分数会得到一个正的标准分数，在平均数之下的分数会得到一个负的标准分数。 z分数是一种可以看出某分数在分布中相对位置的方法。
式中, ：原始数据；
：平均数；s： 标准差。
z = [(value gene X in tumor Y)-(mean gene X in normal)]/(standard deviation X in normal)
# calculate z-scores
scal <- function(x,y){
mean_n <- rowMeans(y) # mean of normal
sd_n <- apply(y,1,sd) # SD of normal
# z score as (value - mean normal)/SD normal
res <- matrix(nrow=nrow(x), ncol=ncol(x))
colnames(res) <- colnames(x)
rownames(res) <- rownames(x)
for(i in 1:dim(x)[1]){
for(j in 1:dim(x)[2]){
res[i,j] <- (x[i,j]-mean_n[i])/sd_n[i]
}
}
return(res)
}
z_rna <- scal(rna_vm[,t_index],rna_vm[,n_index])
# set the rownames keeping only gene name
rownames(z_rna) <- sapply(rownames(z_rna), function(x) unlist(strsplit(x,'\\|'))[[1]])
rm(rna_vm) #we don't need it anymore
生存分析
首先，我们根据上面获得的Z-score 对病人分组， Z > +/- 1.96 大致相当于 p=0.05 或 两倍 SD差距，这些基因我们视为差异基因（比较LFC>1有何不同？）。
# create event vector for RNASeq data
event_rna <- t(apply(z_rna, 1, function(x) ifelse(abs(x) > 1.96,1,0)))
# since we need the same number of patients in both clinical and RNASeq data take the indices for the matching samples
ind_tum <- which(unique(colnames(z_rna)) %in% rownames(all_clin))
ind_clin <- which(rownames(all_clin) %in% colnames(z_rna))
# pick your gene of interest
ind_gene <- which(rownames(z_rna) == 'TP53')
# check how many altered samples we have
table(event_rna[ind_gene,])
# 生存分析用的是time=new_death,status = death_event
# run survival analysis
s <- survfit(Surv(as.numeric(as.character(all_clin$new_death))[ind_clin],all_clin$death_event[ind_clin])~event_rna[ind_gene,ind_tum])
s1 <- tryCatch(survdiff(Surv(as.numeric(as.character(all_clin$new_death))[ind_clin],all_clin$death_event[ind_clin])~event_rna[ind_gene,ind_tum]), error = function(e) return(NA))
# extraect the p.value
pv <- ifelse ( is.na(s1),next,(round(1 - pchisq(s1$chisq, length(s1$n) - 1),3)))[[1]]
# plot the data
plot(survfit(
Surv(as.numeric(as.character(all_clin$new_death))[ind_clin],all_clin$death_event[ind_clin])~event_rna[ind_gene,ind_tum]),
col=c(1:3),
frame=F,
lwd=2,
main=paste('KIRK',rownames(z_rna)[ind_gene],sep='\n'))
# add legend
legend(1800,0.995,legend=paste('p.value = ',pv[[1]],sep=''),bty='n',cex=1.4)
legend(max(as.numeric(as.character(all_clin$death_days)[ind_clin]),na.rm = T)*0.7,0.94,
legend=c(paste('NotAltered=',x1),paste('Altered=',x2)),bty='n',cex=1.3,lwd=3,col=c('black','red'))
image
如果你觉得图不好看可以用survminer，这里重点不是画图。
问题
-
多个基因的生存分析怎么做？
lapply或者for
-
多因素Cox怎么做？
survfit(Surv(as.numeric(as.character(all_clin$new_death)) ind_clin],all_clin$death_event[ind_clin])~event_rna[ind_gene,ind_tum])这里，~ 后面添加不同的自变量即可
可以看到，我们习以为常的生存分析，原来有这么多问题！平时就是照着别人的代码敲，没仔细思考为什么。这篇文章告诉我们，没有一成不变的方法。
