生存分析有个难点是删失(cersored)数据处理，删失数据是指整个数据收集过程都没发生事件的数据。说的有点拗口，举例说我们做某癌症生存分析，那么事件就是因此癌症导致病人死亡。如果有个病人随访期间不幸被车撞死了，那么这个病人记录到此为止，但是我们的事件并没有发生；或者病人突然搬到国外去了，无法继续随访记录，那这个人数据收集结束了，但事件也没有发生。这例子是右删失数据，还有左删失等感兴趣朋友可自行了解，我们一般都是碰到右删失数据。
Kaplan-Meier生存曲线的生存率公式如下，n是存活总人数，d是事件发生数目。
R语言用 survival 和 survminer 包可以很容易实现生存曲线分析与可视化，下面展示一个TCGA数据分析例子。要提醒一下有些教程让人把数据处理成0，1，其实应该把数据处理成1和2，用1表示删失，2表示事件发生(死亡)。
首先准备生存数据，TCGA数据可以直接下载tsv格式临床数据，但是不建议使用，下载的这个数据有不少地方列错位。我曾经下载json版数据用python自己提取，这是个方法。更好方法是去cBioPortal(cBioPortal for Cancer Genomics)下载其整理好的，下载后用excel打开检查一下，还是可能部分条目有问题的。刚刚说把数据转换为1和2，下面就先这样处理一下。
> library(survival, quietly = TRUE)
> library(survminer, quietly = TRUE)
> library(tidyverse, quietly = TRUE)
# 定义转换函数，LIVING(删失)是1，死亡是2. 符号 > 和 + 是因为交互模式，写代码记得去掉。
> statusNum <- function(x){
+ if(x == "LIVING"){
+ return(1)}
+ else{
+ return(2)}
+ }
# 临床数据列非常多，选择我们需要的几列，同时改改列名方便R操作
> clinicaL <- read_tsv("20190211Survival/cBioPortal-blca_tcga_pub_2017_clinical_data.tsv") %>% dplyr::select(`Patient ID`, `Overall Survival (Months)`, `Overall Survival Status`) %>% filter(!is.na(`Overall Survival (Months)`) & !is.na(`Overall Survival Status`)) %>% rename(time=`Overall Survival (Months)`, patient_id=`Patient ID`) %>% mutate(status=map_dbl(`Overall Survival Status`, statusNum)) %>% dplyr::select(patient_id, time, status)
# 数据如下
> head(clinicaL)
# A tibble: 6 x 3
patient_id time status
<chr> <dbl> <dbl>
1 TCGA-2F-A9KP 12.0 2
2 TCGA-2F-A9KQ 94.8 1
3 TCGA-2F-A9KR 105. 2
4 TCGA-2F-A9KT 109. 1
5 TCGA-2F-A9KW 8.34 2
6 TCGA-4Z-AA7M 16.3 1
然后读取基因表达数据，选取有生存数据的病人，根据自己需要基因表达量分组。
> geneExpr <- read_tsv("20190211Survival/ALL_FPKM_UQ_ENTREZID.tsv")
> colnames(geneExpr) %>% head()
[1] "ENSEMBL" "entrezgene_id" "hgnc_symbol" "TCGA-FD-A5BX" "TCGA-K4-A3WS" "TCGA-E7-A8O7"
# 筛选同时有生存数据和表达数据的病人
> patientList <- intersect(clinicaL$patient_id, colnames(geneExpr))
> head(patientList)
[1] "TCGA-2F-A9KP" "TCGA-2F-A9KQ" "TCGA-2F-A9KR" "TCGA-2F-A9KT" "TCGA-2F-A9KW" "TCGA-4Z-AA7M"
> length(patientList)
[1] 399
# 为了后面方便，直接把表达数据样本按照临床数据样本顺序排列
> clinicaL2 <- dplyr::filter(clinicaL, patient_id %in% patientList) %>% dplyr::distinct(patient_id, .keep_all = TRUE)
> geneExpr2 <- dplyr::select(geneExpr, hgnc_symbol, clinicaL2$patient_id)
> dim(clinicaL2)
[1] 399 3
> dim(geneExpr2)
[1] 20084 400
# 我这里就随便选个基因
> head(genE)
hgnc_symbol TCGA-2F-A9KP TCGA-2F-A9KQ TCGA-2F-A9KR TCGA-2F-A9KT TCGA-2F-A9KW
"SCFD2" "17.6915874147098" "17.1654098740173" "17.7186146961568" "16.8297194432778" "16.3051962379739"
# 因为之前表达数据样本排列根据临床样本顺序来的，我就直接把表达数据添加到临床数据表格了
> clinicaL3 <- mutate(clinicaL2, SCFD2=genE[-1]) %>% dplyr::arrange(desc(SCFD2)) %>% mutate(SCFD2_Expr=c(rep("High", 200), rep("Low", 199)))
> head(clinicaL3, n=3)
# A tibble: 3 x 5
patient_id time status SCFD2 SCFD2_Expr
<chr> <dbl> <dbl> <chr> <chr>
1 TCGA-BT-A20X 8.25 2 18.5989977626721 High
2 TCGA-E7-A85H 12.9 1 18.4226732766918 High
3 TCGA-FD-A3SJ 24.3 2 18.4119859441587 High
> tail(clinicaL3, n=3)
# A tibble: 3 x 5
patient_id time status SCFD2 SCFD2_Expr
<chr> <dbl> <dbl> <chr> <chr>
1 TCGA-LT-A5Z6 15.6 1 15.6753402657604 Low
2 TCGA-FJ-A3ZF 17.2 1 15.5461411355296 Low
3 TCGA-ZF-A9R7 21.8 1 14.6207572703697 Low
现在把样品分组信息添加到了临床信息表，下面用 survival 和 survminer 分析按SCFD2表达分组生存是否差异。
fit <- survfit(Surv(time, status) ~ SCFD2_Expr, data = clinicaL3)
# 画图
ggsurvplot(fit, pval = TRUE)
2.png
从这结果看按照SCFD2基因表达分2组生存率无显著差异。对于图片美化， ggsurvplot 产生图片是ggplot2对象，可以直接用ggplot2进行任意修改。
> p <- ggsurvplot(fit, pval = TRUE)
# 换个主题把坐标翻转
> p$plot <- p$plot + theme_dark() + coord_flip()
> p$plot
3.png
参考资料
Survival Analysis in R
survminer R package: Survival Data Analysis and Visualization - Easy Guides - Wiki - STHDA
网友评论