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

2020-09-04 14:33:03 浏览数 (1)

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

数据

文章中说明使用了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）是指某段时间开始时存活的个人至该时间结束时仍然存活的可能性大小
生存概率 = 1 — frac{某人群活过某段时间例数}{该人群同时段期观察例数}
  • 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语言知识啦：

代码语言：javascript复制
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老师写的，我只是一个代码搬运工以及修修补补。

  • 下面是生存分析
代码语言：javascript复制
# 开始生存分析
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,  )
代码语言：javascript复制
# 接下来绘制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)处理，不过表达量本身仅仅是用来把病人分组，何种归一化并不会影响基因表达量高低这个排序。

dfs free na time

0 人点赞