首先是在自己服务器里安装R，这一步可以直接借用conda来完成，推荐用jupyter notebook来运行R简单好用~，具体设置的可以参考“jupyter notebook中添加 R内核"
测试数据和代码见：https://github.com/leiwaaping/survival-analysis-with-TCGA-dataset
以下是文字介绍和部分内容我自己的理解说明：
做这个分析的时候，生存分析相关的一个包RTCGA刚好升级到只兼容3.6版本以上的R，然而我的R还是以前安装的旧版，所以一直报错，中间折腾了好久也只是报告说某个包无法安装完成，特定安装了被报有问题的包也依然无法解决问题，于是回到bioconda的文件里看每个包的限制条件，才发现是版本问题，所以吸取教训，看报错和找到真正的问题点很重要。
简单来说生存分析的就是已知一群人的生存情况（有可能是只活了一段时间，有可能是一直活着的），现在得知了这群人的另外一些数据，依着这些数据对这群人进行分组，然后每组的人分别进行生存率作图（横坐标是生存的时间长度，纵坐标是这段时间的这一组人的生存率，如一开始有100人，30天后10人死亡，则此时生存率从100%下降到90%）如果中途有人的数据丢失，病人不再回到医院复诊无法确认是否仍然生存，则记一个“+”在图上标出。一般来说最好有200个以上的样本才比较好做生存分析（考虑数据缺补情况）。生存分析的P值是指当前结果是随机出现的概率，P值越小，结果越可信
生存分析概念
在病人的临床信息中心，会有“临床试验终点"(end point)，根据实验设计不一样，又可以分为不同的研究指标（OS，PFS，DFS，TTP），根据研究目的不同可以选择不同的终点。
Overall survival（OS，总生存）: 0=存活；1=任何原因的死亡；t=随机化开始至死亡的时间。 这个是肿瘤临床试验中的最佳指标，因为好记录。临床试验中常用到5年OS生存率，肿瘤患者经过各种中和治疗后，生存五年以上的患者比例。也有用10年生存率来表示疗效的。
Progression-free survival（PFS，无进展生存期）：0=存活；1=肿瘤发生任何方面的进展（复发转移）或者任何原因的死亡；t=随机化开始至肿瘤进展/死亡的时间。常用的有3年PFS生存分析，t>3的患者，统一t=3，p=0。PFS相比OS包含了恶化这个概念，可用于评估一些治疗的临床效益。
Time to progress（TTP，疾病进展时间）：0=存活；1=肿瘤发生任何方面的进展（复发转移）。质保函了肿瘤恶化，不包含死亡。t=随机化开始至肿瘤进展的时间
disease-free survival（DFS，无病生存期）：0=存活；1=疾病复发或者任何原因的死亡；t=随机化开始至疾病复发/死亡的时间。
Censoring（删失）：这经常会在临床资料中看到，生存分析中也有其对应的参数，一般指不是由死亡引起的的数据丢失，可能是失访，可能是非正常原因退出，可能是时间终止而事件未发等等，一般在展示时以‘+’号显示，可分为：
- left censored（左删失）：只知道实际生存时间小于观察到的生存时间
- right censored（右删失）：只知道实际生存时间大于观察到的生存时间
- interval censored（区间删失）：只知道实际生存时间在某个时间区间范围内
基础版生存分析画图
已知样本的生存时间，当前状态（overall survival time/ PFS（按肿瘤进展来划分，考虑了死亡/复发等情况）），和分组数据，有一些可能要根据某数值大小来分组，可选用下文更详细的分组代码。有一些数据直接就分好组了，可以直接整合并画图
library(survival)
library(ggplot2)
library(ggpubr)
library(magrittr)
#install.packages("survminer")
library(survminer)
library(dplyr)
groupedclin<-as.data.frame(info)
colnames(groupedclin)
colnames(groupedclin)=c("overall_survival_Event","Overall_Survival_Time","PFS","PFS_T","group")
head(groupedclin)
# How many samples of each type?
table(groupedclin$group)
# Tabulate by outcome
xtabs(~group+overall_survival_Event, data=groupedclin) %>% addmargins()
coxph(Surv(Overall_Survival_Time, overall_survival_Event)~group, data=groupedclin)
sfit <- survfit(Surv(Overall_Survival_Time, overall_survival_Event)~group, data=groupedclin)
ggsurvplot(sfit, conf.int=TRUE, pval=TRUE)
cox比例风险分析
- The Cox proportional-hazards model (Cox, 1972) Cox比例风险模型是常用与医学研究中的回归模型，通常用于研究预测变量与生存时间的关系
- Kaplan-Meier curves 与 logrank test tests属于单因素分析的例子，他们研究的是单一变量与生存的关系,并且Kaplan-Meier 与log-rank检验只适用于分类变量，却并不适用于数值型变量，比如我们常见的基因表达
- Cox比例风险回归分析的优势在于可以分析分类变量与数值变量，并且将生存分析的范围由单因素拓展到多因素的分析.
cox比例风险分析，同时评估几个因素对生存的影响。（如预后治疗方法、年龄、性别同时纳入考虑生存率）
#PFS=status PFS_T = time
#单因素Cox回归分析
res.cox <- coxph(Surv(PFS_T, PFS) ~ treatmentArm, data = info)
summary(res.cox)
#批量进行单因素Cox生存分析
covariates <- c("age", "gender", "TNstage", "SMART","group")
univ_formulas <- sapply(covariates,function(x) as.formula(paste('Surv(PFS_T, PFS)~', x)))
#univ_formulas
univ_models <- lapply( univ_formulas, function(x){coxph(x, data = info)})
#univ_results
#要整理输出，注意要一个因素一行才行，像treatmentArm这样子chr类型的，会被当成不同种类而有不同的p-value（见下文多因素的输出）
#chr类型会被堪称不同种类，所以要变成numeric（看age）
univ_results <- lapply(univ_models,
function(x){
x <- summary(x)
p.value<-signif(x$wald["pvalue"], digits=2)
wald.test<-signif(x$wald["test"], digits=2)
beta<-signif(x$coef[1], digits=2);#coeficient beta
HR <-signif(x$coef[2], digits=2);#exp(beta)
HR.confint.lower <- signif(x$conf.int[,"lower .95"],2)
HR.confint.upper <- signif(x$conf.int[,"upper .95"],2)
HR <- paste0(HR, " (",
HR.confint.lower, "-", HR.confint.upper, ")")
res<-c(beta, HR, wald.test, p.value)
names(res)<-c("beta", "HR (95% CI for HR)","wald.test", "p.value")
return(res)
#return(exp(cbind(coef(x),confint(x))))
})
res <- t(as.data.frame(univ_results, check.names = FALSE))
as.data.frame(res)
同时考虑多个因素，看这些因素如何共同影响生存时间
#多因素Cox生存分析
res.cox <- coxph(Surv(PFS_T, PFS) ~ x + treatment + age +gender + STAGE + batch, data = dat)
res.cox
如果需要每个基因逐一放入cox模型中，检查基因的风险度，同时增加协变量来矫正
library(survival)
#（100位数以内）取消科学计数法
options(scipen = 100)
coxR<-lapply(11:(dim(dat)[2]),function(x){
#coxR<-lapply(11:19,function(x){
n=dat[,x]
res.cox <- coxph(Surv(PFS_T, PFS) ~ n + treatment + age +gender + STAGE + batch, data = dat) # complex Coxph, ori expr
s <- summary(res.cox)
gene<-colnames(dat)[x]
#也可以单独提取信息，但千万要注意选定好p-value 和HR，不同特征的位置不一样
p.value<-signif(signif(s$coefficients)[1,5], digits=3)
HR <-exp(signif(s$coefficients[1], digits=3))
res<-c(gene,HR,p.value)
return(res)
})
res <- t(as.data.frame(coxR, check.names = FALSE))
row.names(res)=res[,1]
res=as.data.frame(res[,-1])
res[] <- lapply(res, function(x) {
if(is.factor(x)) as.numeric(as.character(x)) else x
})
colnames(res)=c("HR","pvalue")
res_p0=res[res[, "pvalue"]<0.05,]
res_p0=res_p0[order(res_p0[,"HR"],decreasing=TRUE),]
head(res_p0)
Cox回归分析代码参考：
R语言生存分析03-Cox比例风险模型
Cox分析summary结果含义：
Hazard ratios. The exponentiated coefficients (exp(coef) ), 瞬时“掉血”速度，一般用来衡量药物的治疗效果。约远离1表示干预效果越大（有利或有害两个方向）
Cox Proportional-Hazards Model
下面的其实可以不用看了，是专门针对TCGA的样本，过滤并人工分组，查看生存率的一个例子。
拓展：提取符合要求的TCGA的样本生存数据
TCGA数据库和样本命名规则可以参考：https://zhuanlan.zhihu.com/p/46520322
1.TCGA数据有可能出现一个人多次检测多个样本的情况，由于生存分析一人对应一个数据，所以一般酌情删除后面测试的数据，只保留该名患者第一次的数据。如果是用自己的数据来算，有可能出现这种情况，可以用“去重”操作
2.TCGA的样本名字中间一般用“-”间隔，但不排除自己实验得到的数据样本名用了别的间隔符，要注意修正过来避免后期没有一个样本匹配上的。
expr<-read.table("doc/JCI121476.sdt12.txt", header = TRUE, stringsAsFactors=F)
#sort to remove duplicate
expr<-expr[order(expr$Sample_ID,decreasing = T),]
expr$Sample_ID<-substr(expr$Sample_ID,1,12)
dim(expr)
expr[1:2,]
expr<-expr[!duplicated(expr$Sample_ID, fromLast=TRUE), ]
dim(expr)
expr$Sample_ID<-gsub("_", "-",expr$Sample_ID)
expr[1:2,]
3. 删除低质量样本
一般看participant来确认是哪一个人，很多数据集会只显示前四个数来区别样本。需要注意的是，最后一个数“01C”，01~09表示肿瘤，10~19表示正常对照，一般01和11最常见。还有区分转移性和复发性的数据，详见https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables/sample-type-codes
后面的vial指在一系列患者组织中的顺序，绝大多数样本该位置编码都是A; 很少数的是B，表示福尔马林固定石蜡包埋组织，一般我只保留A的数据，因为对B、C等数据没有了解，无法保证数据质量。一般B和C的不多，我都手动检查和删除了。
4.样本合并
并不是所有TCGA的样本都集合了分组数据和生存数据的，用merge来依据样本名来进行数据合并。（如图，有457个样本的该疾病的生存数据，但是只有277个样本既有herv表达量又有生存数据。）
library(RTCGA)
rownames(infoTCGA())
library(RTCGA.clinical)
library(RTCGA.mRNA)
#expr[,which(colnames(expr)=="herv_2965")]
singleherv<-expr[,c(1,which(colnames(expr)=="herv_2965"))]
#singleherv[1:4,]
clin <- survivalTCGA(COAD.clinical,extract.cols="admin.disease_code")
#Show the first few lines
dim(clin)
clin[1:4,]
colnames(clin)[2] <- "Sample_ID"
clinsample<-merge(clin,singleherv,by="Sample_ID",all=FALSE)
colnames(clinsample)[5] <- "expr"
dim(clinsample)
clinsample[1:4,]
进行生存分析
library(survival)
library(ggplot2)
library(ggpubr)
library(magrittr)
library(survminer)
#分组
library(dplyr)
a1<-as.numeric(mean(clinsample$expr))
g1<-clinsample[which(clinsample$expr<a1),]
g1$group=1
#g1[1:4,]
g3<-clinsample[which(clinsample$expr>=a1),]
g3$group=3
#g3[1:4,]
groupedclin<-as.data.frame(rbind(g1,g3))
# How many samples of each type?
table(groupedclin$group)
# Tabulate by outcome
xtabs(~group+patient.vital_status, data=groupedclin) %>% addmargins()
- 基本来说就是依照数据，再原来的临床数据后面加一个分组，分组可以用数字表示，也可以用文字表述，status 0 /1 表示生存与否。如下面：
2. 关于分组，可以采取文字（如上图），可以采取众数，可以采取平均值（如果有很多0的情况），也可以采取上/下四分位数的情况，看自身数据决定。（生存分析图只是一个辅助说明的作用，有些参数调整可能得到的图结果就不一样了，具体需要一个怎么样的图，可以再合适的范围内酌情调整优化）
#看一下数据分布
#see score distributation
hist(individual.socre$Score,breaks=100)
#取上/下四分位数或众数
quantile(individual.socre$Score)[2]
quantile(individual.socre$Score)[4]
as.numeric(strsplit(as.character(quantile(individual.socre$Score)[2]), ":")[[1]])
as.numeric(strsplit(as.character(quantile(individual.socre$Score)[4]), ":")[[1]])
median(individual.socre$Score)25%: 6479.38306575775%: 8043.83190523448
6479.383065757
8043.83190523448
7134.55590193063
画图
coxph(Surv(times, patient.vital_status)~group, data=groupedclin)
sfit <- survfit(Surv(times, patient.vital_status)~group, data=groupedclin)
summary(sfit, times=seq(0,365*5,365))
ggsurvplot(sfit, conf.int=TRUE, pval=TRUE)
5157
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
