|
Freescience由浙江大学医学院几个硕博士发起创建，旨在最广泛分享有价值的科研技能和知识；FreeScience的宗旨：“科学自由分享、人人平等，共求真理”。 进入正题之前，咱们先来复习一下R语言（已经学过的孩子跳过直接往下翻~） 生存分析是医学中常见的分析方法，一般指对研究群体生存时间的分布情况进行描述、刻画。一般对于医学数据的生存时间进行分析和推断，研究生存时间和相关影响因素间关系及其程度大小，并对其生存时间进行预测。 生存分析，大多就是说的KM方法估计生存函数，并且画出生存曲线，然后还可以根据分组检验一下它们的生存曲线是否有显著的差异。在R中，有个包survival做生存分析就很方便！只需要记住三个函数： Surv：用于创建生存数据对象 survfit：创建KM生存曲线或是Cox调整生存曲线 survdiff：用于不同组的统计检验 我们先来举个栗子。 对于上面的数据，我们用下面的代码做生存分析！ library(survival) my.surv <- surv(os_months,os_status=""> ##这个生存对象是看看病人的总生存期与死亡状态的关系 ##这个Surv函数第一个参数必须是数值型的时间，第二个参数是逻辑向量，1,0表示死亡与否 kmfit1 <-> ## 直接对生存对象拟合生存函数 summary(kmfit1) plot(kmfit1) ##画出生存曲线 survdiff(my.surv~type, data=dat) ### 根据生存对象再加上一个分组因子来拟合生存函数，并且比较不同因子分组的生存效果。 TCGA数据里面的生存分析例子 然而，我们已经知道了生存分析，是随着时间的流逝，死亡率是如何增加的。那么如何根据某些因子把样本分组，看到他们死亡率的变化趋势显著的不同，从而说明这个因子是非常有效的分类方式，这个因子可以是一个biomarker呢？ 我们可以用cox模型来分析这个因子是如何影响生存函数的，我们在举一个TCGA数据库的栗子进行说明，这个栗子的是TCGA里面卵巢癌的数据，根据甲基化数据分成了4个组，我们就下载四个组样本的临床数据。 数据是用cgdsr下载的 （cgdsr说明：http://www./?p=1257） library(cgdsr) mycgds <->http://www./public-portal/‘) test(mycgds) all_TCGA_studies <-> all_tables <- getcaselists(mycgds,=""> all_dataset<- getgeneticprofiles(mycgds,=""> BRCA1 <- getprofiledata(mycgds,="" my_gene,="" my_dataset,=""> ov_tcga_pub_meth1<- getclinicaldata(mycgds,=""> ov_tcga_pub_meth2<- getclinicaldata(mycgds,=""> ov_tcga_pub_meth3<- getclinicaldata(mycgds,=""> ov_tcga_pub_meth4<- getclinicaldata(mycgds,=""> 根据甲基化数据，把癌症病人分成了4组，我们的临床数据记录了13项，但是我们只需要用到OS_MONTHS和OS_STATUS就可以来估计KM生存函数，画出生存曲线啦！ 无病生存期（Disease-free survival,DFS）的定义是指从随机化开始至疾病复发或由于疾病进展导致患者死亡的时间。该指标也常作 为抗肿瘤药物III期临床试验的主要终点。某些情况下，DFS与OS相比，作为终点比较难以记录，因为它要求认真随访，及时发现疾病复发，而且肿瘤患者的 死亡原因也很难确定（16）。肿瘤患者常有合并症（如，心血管病），这些合并症可能会干扰对DFS的判断。并且，肿瘤患者常死于医院外，不能常规进行尸检。 总生存期（Overall survival,OS）的定义是指从随机化开始至因任何原因引起死亡的时间。该指标常常被认为是肿瘤临床试验中最佳的疗效终点。如果在生存期上有小幅度的提高，可以认为是有意义的临床受益证据。作为一个终点，生存期应每天进行评价，可通过在住院就诊时，通过与患者直接接触或者通过电话与患者交谈，这些相对比较容易 记录。确认死亡的日期通常几乎没有困难，并且死亡的时间有其独立的因果关系。当记录至死亡之前的失访患者，通常截止到最后一次有记录的、与患者接触的时间。 library(survival) attach(ov_tcga_pub_meth1) ## 估计KM生存曲线 my.surv <- surv(os_months,os_status=""> kmfit1 <-> summary(kmfit1) plot(kmfit1) ##我们很容易看到，随着感染癌症的时间延长，病人的死亡率到了一定程度就不增加了，有些病人熬过了这个癌症，或者说，到我们拿到数据为止，他们还没有死亡。 ## 随便取一根因子来分组TUMOR_STAGE_2009估计KM生存曲线 kmfit2 <-> summary(kmfit2) plot(kmfit2,col = rainbow(length(unique(ov_tcga_pub_meth1[,13])))) detach(ov_tcga_pub_meth1) ##可以看到，我们根据病人的TUMOR_STAGE_2009把他们分成了这些组，不同组的生存曲线不一样，但是我们肉眼无法看出它们这些组直接的生存率是否有显著差异！我们需要做统计检验！ 我们就不对这个进行检验了，我们还是用下面的代码来对甲基化数据的分组来做检验，看看我们的分组是否有效！ ov_tcga_pub_meth1$sample<-> ov_tcga_pub_meth2$sample<-> ov_tcga_pub_meth3$sample<-> ov_tcga_pub_meth4$sample<-> ov_tcga_pub_meth1$type<-> ov_tcga_pub_meth2$type<-> ov_tcga_pub_meth3$type<-> ov_tcga_pub_meth4$type<-> dat=rbind(ov_tcga_pub_meth1,ov_tcga_pub_meth2,ov_tcga_pub_meth3,ov_tcga_pub_meth4) attach(dat) ## 根据分组type估计KM生存曲线 my.surv <- surv(os_months,os_status=""> kmfit3 <-> summary(kmfit3) plot(kmfit3,col = rainbow(4)) ##由图中可以看到甲基化数据分成的4个组，生存率差异还是蛮大的 ##用图形方法检验PH假设 plot(kmfit3,fun=’cloglog’,col = rainbow(4)) # 检验显著性 survdiff(my.surv~type, data=dat) # 用strata来控制协变量的影响 survdiff(my.surv~type+strata(TUMOR_STAGE_2009),data=dat) ##然后我们用这个代码来检验，是否显著，结果发现，p值并没有小于0.05，只能说是可以这样分组，但是分组的效果没有我们想象中的那么好！ names(dat) detach(dat) jimmy实名入驻研论问答平台，如何在研论上邀约jimmy回答你的问题呢？ 第一步：本公众号科研实战菜单栏进入研论问答平台
第三步：直接点击头像可以看到邀约信息和被隐藏的聊天方式哦 第四步：邀约成功后，你将看到对方的深入沟通的联系方式，同时，对方默认回答你最近的至少一个提问（由公众号短信形式提示对方） 生信菜鸟专栏 生信菜鸟专栏是生信技能树论坛的版主团队的专栏，团队成员生信技能背景丰富，文件格式，数据资源，软件使用，脚本技巧，统计绘图，组学实战均有对应人才。而本专栏将从基础到深入，为零基础的各位剖析生信技能。 |
|
|
