学徒数据挖掘之谁说生存分析一定要按照表达量中位值或者平均值分组呢？

学徒和学员已经陆续出师，是时候把生信技能树的舞台交给后辈了！

下面是《生信入门第6期》学员的分享

呃提笔来写这个总结，貌似是上周发生的事了，但是不太想回忆，学习新东西会抱着满满的热情，但是翻旧账就会让人心里不是很舒畅...特别是没能完成的旧事重提。唔但是也算是慢慢习惯了完成一件事会记下点个人感悟，nice~这是个好习惯（最开始当然也只是Jimmy大大的任务模式啦）。

这里感谢老师费心帮助我完成这个，还找了位前辈小蚂蚁山神（就是昨天完成任务的那个：学徒数据挖掘之IgG4-RD患者的WGCNA分析）来帮助我，非常感谢大家的帮助:happy:

写完这个，又发生了一件特别让我感动的事情，不管多小的人物类似于我，Jimmy大大都会认真的回答你的小小疑问，遇到一个好老师真的很重要，超开心能认识这么棒的老师和一些志同道合的小伙伴。虽然学习生信的初衷可能只是为了发一篇简单的小文章好找工作一点，但是不知不觉中每天学习一点生信技能树的公众号，论坛和各种这个团队优秀伙伴的简书笔记成为了一个习惯，习惯着等任务，没有任务时自己找事情做，B站不再是天天刷吃播视频哈哈哈，B站真的能学习，老师的那么多视频得刷到啥时候啊，还不断更新~收藏的东西也不会吃灰:six::six::six:。害，就是我太懒了，不愿意拓展知识。

首先接到完成文章中一个生存分析的图的任务：

数据

文章中说明使用了cbioportal上的TCGA的数据：

• 这是方法和描述里的文字信息。作者利用TCGA里的LIHC数据做了生存分析，其中，一共有370个肝癌患者的资料，总生存期（Overall Survival， OS）生存分析用到了370个患者资料，无病生存期(disease-free survival )生存分析用到了319个患者资料，患者资料不一致可能是删去了缺失值。根据作者的描述，是以LHPP的表达量z~score后，<= -1为界，如果按照百分比的话大概是20%-25%

原文中的图为这样：

• 第一次接触这种...没怎么看懂无病生存期和总生存期的信息，搜索了背景知识。刚开始压根没想到要用网页工具重复这个图，毕竟身为一个任务，不该这样子完成。

基本术语：

• Event（事件）：在癌症研究中，事件可以是Relapse，Progression以及Death

• Survival time（生存时间）：一般指某个事件的开始到终止这段事件，如癌症研究中的疾病确诊到缓解或者死亡，其中有几个比较重要的肿瘤临床试验终点：

• 中位生存期又称为半数生存期，表示有且只有50%的个体可以活过这个时间。

• 比如共1000人参加临床试验，将每个人的生存时间按从小到大排名，第501人的生存时间为18个月，即表明该临床试验的中位生存期为18个月。

• 如果是评估某个癌种的中位生存期，一般从发现该肿瘤开始计算；如果是评估某项临床试验的中位生存期，一般从给药或随机开始。

• 总生存期（Overall Survival，OS）：指从随机化开始到任意原因死亡的时间(非肿瘤因素引起的死亡也被统计在内，比如受试者在统计时间内车祸身亡，其生存期的数据也属于有效数据。)，我们一般见到的5年生存率、10年生存率都是基于OS的

• 无进展生存期（progression-free survival，PFS）：指从开始到肿瘤发生任意进展或者发生死亡的时间；受试者只要“肿瘤恶化”或“死亡”二者其一先发生，则到达研究终点。疾病进展是指肿瘤增长，或肿瘤原发病灶转移，或发现新的病灶）等。相比OS包含了恶化这个概念，可用于评估一些治疗的临床效益

• 疾病进展时间（time to progress，TTP）：从开始到肿瘤发生任意进展或者进展前死亡的时间；TTP相比PFS只包含了肿瘤的恶化，不包含死亡。故若受试者尚未发生“肿瘤恶化”就已经先“死亡”，则此为受试者再也观察不到“肿瘤恶化”

• 无病生存期(disease-free survival ) ：指从开始到肿瘤复发或者任何原因死亡的时间；常用于根治性手术治疗或放疗后的辅助治疗，如乳腺癌术后内分泌疗法等：

• event free survival（EFS，无事件生存期）：指从开始到发生任何事件的时间，这里的事件包括肿瘤进展，死亡，治疗方案的改变，致死副作用等（主要用于病程较长的恶性肿瘤、或该实验方案危险性高等情况下）

• 中位生存期（Median Survival Time，MST）

• 生存概率（Survival probability）是指某段时间开始时存活的个人至该时间结束时仍然存活的可能性大小

• Censoring（删失）：这经常会在临床资料中看到，生存分析中也有其对应的参数，一般指不是由死亡引起的的数据丢失，可能是失访，可能是非正常原因退出，可能是时间终止而事件未发等等，一般在展示时以‘+’号显示
• left censored（左删失）：只知道实际生存时间小于观察到的生存时间
• right censored（右删失）：只知道实际生存时间大于观察到的生存时间
• interval censored（区间删失）：只知道实际生存时间在某个时间区间范围内

用在线xena下载数据，直接下载临床信息，全部都是整理好的，分14个数据集的和19个数据集的，19的那个。

xena数据库中会有两个数据：

• GDC（HTseq count数据中，Ensemble ID）

• TCGA（gene symol）

 

• LIHC_survival.txt.gz

 

• LIHC_clinicalMatrix

 

有了这么多数据文件，接下来就是该秀出我在生信技能树学习班获得的R语言知识啦：

rm(list=ls())
options(stringsAsFactors = F)
##读取数据
mRNA_HiSeqV2 = read.table('HiSeqV2.gz',header = T,sep = '\t',check.names = F,row.names = 1)
dim(mRNA_HiSeqV2)
mRNA_HiSeqV2[1:4,1:4]

#查看NA的数据
dim(mRNA_HiSeqV2)
exp = mRNA_HiSeqV2
exp = exp[apply(exp, 1, function(x) sum(x > 1) > 9), ]
dim(exp)
exp[1:4,1:4]


##临床信息（在这里没有用到）
mRNA_clinical = read.table('LIHC_clinicalMatrix',header = T,sep = '\t',row.names = 1)
dim(mRNA_clinical)
mRNA_clinical[1:4,1:4]


##生存相关信息
mRNA_survival = read.table('LIHC_survival.txt.gz',header = T,sep = '\t',row.names = 1)
dim(mRNA_survival)
mRNA_survival <- mRNA_survival[, -1] 
mRNA_survival[1:4,1:4]


save(mRNA_survival,mRNA_clinical,exp,group_list,file='LIHC.Rdata')


#01–09是癌症，10–19是正常，20–29是癌旁
根据样本ID的第14-15位，给样本分组（tumor和normal）
library(stringr)
table(str_sub(colnames(exp),14,15))
group_list = ifelse(as.numeric(str_sub(colnames(exp),14,15)) < 10,'tumor','normal')
group_list = factor(group_list,levels = c("normal","tumor"))
table(group_list)


##xena下载的数据是经过log的count值，所以要还原才可以进行差异分析
exprSet = exp
exprSet <- 2^exprSet-1
dim(exprSet)
exprSet[1:4,1:4]

#并且还要取整数，现在就变成count值了
exprSet <- floor(exprSet)
exprSet[1:4,1:4]


save(mRNA_survival,mRNA_clinical,exp,group_list,exprSet,file='LIHC1.Rdata')



rm(list=ls())
options(stringsAsFactors = F)
load('LIHC1.Rdata')
library("stringr")


先看以count数据的结果 OS
tumor_sample <- colnames(exprSet)[substr( colnames(exprSet),14,15) < 10] 
mRNA_survival = mRNA_survival[,-1]
pheno = rownames(mRNA_survival)[substr(rownames(mRNA_survival),14,15) < 10] 
fin_tumor = intersect(tumor_sample,pheno)

# 提取表达矩阵和临床信息 
exprSet = exprSet[,fin_tumor] 
pheno_fina = mRNA_survival[fin_tumor,] 

### 这里又归一化了？
exprSet = log2(exprSet+1)

代码超级长，绝大部分都是生信技能树jimmy老师写的，我只是一个代码搬运工以及修修补补。

• 下面是生存分析

# 开始生存分析
library(ggrisk) 
library(survival) 
library(survminer)

# 先做无病生存期，首先把DFI不明（NA）的样本给去掉 
DFI =  pheno_fina[,3:6] 
dfi = DFI[!is.na(DFI$DFI),] 
dfi_count = exprSet[,rownames(dfi)]
dd = scale(as.numeric( dfi_count["LHPP",]))

# 绘图 
mySurv<-Surv(dfi$DFI.time, dfi$DFI) 
### 以LHPP的表达量z-score后，<= -1为界，哦z-score是这样用的啊
dfi_group<-ifelse(dd <= -1 ,"low","high")
fi <-survfit(mySurv~dfi_group,data = dfi) 
survminer::ggsurvplot(fi, data = dfi, risk.table = TRUE, 
                      pval = T,  )

 

# 接下来绘制os的 
os = pheno_fina[,c(1,2,7,8)] 
os = os[!is.na(os$OS),] 
os_count =  exprSet[,rownames(os)]

dd = scale(as.numeric( os_count["LHPP",]))

mySurv<-Surv(os$OS.time, os$OS) 
dfi_group<-ifelse(dd <= -1 ,"low","high")
fi <-survfit(mySurv~dfi_group,data = dfi) 
survminer::ggsurvplot(fi, data = dfi, risk.table = TRUE, 
                      pval = T,  )

 

• 可以看出和原文的图还是挺像的。

这大概算是完成了

• 要知道一些生存分析背景知识
• 明白OS和DFS的定义，这两者使用的数据是不一样的，处理标准也不一样
• 对于OS来说，如果时间这一栏位NA，应该把NA替换成大生存时间（OS.time那一列的大值），而对于DFS来说，要手动处理结局事件，如果结局事件为NA，可考虑剔除。一方面可以直接利用TCGA的临床数据，自行制作OS，DFS的信息（死亡时间 - 诊断时间= OS.time，对于DFS.time来说也是类似的，复发时间或者死亡时间 -诊断时间 = DFS.time，根据定义去计算相应 的时间），也可直接去XENA上面下载已经处理好的临床数据。
• 需要找准分组的cutoff。原文指出临界标准为：LHPP的表达量经Z~score后，≤ -1 ，则为低表达， 否则为正常表达（高表达）。常见的分组方式可以包括：中位数、上/下四分位数，上/ 下三分位数，或者和原文一样用Z~score，再取一侧或双侧mean +- 1/2/3 * sd。

原文没有指明，他用的表达数据是Count，标准化后的Count，FPKM还是TPM，也 没有说是否对数据进行Log2(dat + 1)处理，不过表达量本身仅仅是用来把病人分组，何种归一化并不会影响基因表达量高低这个排序。

写到最后

如果你时间比较充裕，而且确实想把生存分析学会，学好，我在生信技能树多次分享过生存分析的细节也值得你认真读完；

“ID转换和生存分析”群的钉钉群号： 35371384，如果你下载钉钉软件，搜索进群这个能力都没有，我还是建议你放弃学生信哈！

如果你也想开启自己的生物信息学数据处理生涯，加入生信技能树小圈子，还等什么呢，赶快行动起来吧！ 发邮件（[e-mail]）给生信技能树创始人jimmy就有惊喜哦！当然了，不能是辣鸡或者骚扰邮件啦，带上自己的简历和想学习交流的诚心吧！