• 主页
  • 课程

    关于课程

    • 课程归档
    • 成为一名讲师
    • 讲师信息
    同等学历教学

    同等学历教学

    免费
    阅读更多
  • 特色
    • 展示
    • 关于我们
    • 问答
  • 事件
  • 个性化
  • 博客
  • 联系
  • 站点资源
    有任何问题吗?
    (00) 123 456 789
    weinfoadmin@weinformatics.cn
    注册登录
    恒诺新知
    • 主页
    • 课程

      关于课程

      • 课程归档
      • 成为一名讲师
      • 讲师信息
      同等学历教学

      同等学历教学

      免费
      阅读更多
    • 特色
      • 展示
      • 关于我们
      • 问答
    • 事件
    • 个性化
    • 博客
    • 联系
    • 站点资源

      未分类

      • 首页
      • 博客
      • 未分类
      • (未测试)对TCGA数据库的任意癌症中任意基因做生存分析

      (未测试)对TCGA数据库的任意癌症中任意基因做生存分析

      • 发布者 weinfoauthor
      • 分类 未分类
      • 日期 2020年2月6日
      • 评论 0评论

      参考网址:https://cloud.tencent.com/developer/article/1167640

      本教程目录:

      • 首先使用cgdsr获取表达数据集临床信息
      • 临床资料解读
      • 简单的KM生存分析
      • 有分类的KM生存分析
      • 根据基因表达量对样本进行分组做生存分析
      • cox生存分析
      • 某基因突变与否也可以用来分组
      • 基因的拷贝数也可以进行分组
      • 批量进行分组并且做生存分析

      生存分析一般来说是针对RNA表达数据,可以说mRNA-seq的转录组数据,也可以说miRNA-seq数据,或者基因表达芯片的表达量值。

      生存分析,大多就是说的KM方法估计生存函数,并且画出生存曲线,然后还可以根据分组检验一下它们的生存曲线是否有显著的差异。

      在R中,有个包survival做生存分析就很方便!只需要记住和熟练使用三个函数:

      • Surv:用于创建生存数据对象
      • survfit:创建KM生存曲线或是Cox调整生存曲线
      • survdiff:用于不同组的统计检验

      首先使用cgdsr获取表达数据集临床信息

      既然是要说明如何对任意癌症的任意基因做生存分析,那么我们首先需要理解cgdsr下载TCGA任意数据的用法(见之前的教程),下面的例子是获取TCGA数据库的乳腺癌的BRCA1和BRCA2基因的表达,以及涉及到的病人的临床资料。

      rm(list = ls())
      library(cgdsr)
      library(DT) 
      mycgds <- CGDS("http://www.cbioportal.org/public-portal/")
      mycancerstudy = 'brca_tcga'  
      ## 下面的代码可以不需要运行,因为已经保存好了用来做生存分析的数据。
      ### 但是需要看懂代码,这样才能做任意癌症的任意基因的任意数据的生存分析;
      
      
      if(F){
        getCaseLists(mycgds,mycancerstudy)[,1]
        getGeneticProfiles(mycgds,mycancerstudy)[,1]
        mycaselist ='brca_tcga_rna_seq_v2_mrna'  
        mygeneticprofile = 'brca_tcga_rna_seq_v2_mrna'  
        choose_genes=c('BRCA1','BRCA2')
        ## get expression data
        expr=getProfileData(mycgds,choose_genes,
                            mygeneticprofile,mycaselist)
        ## get mutation data  
        mut_df <- getProfileData(mycgds, 
                                 caseList ="brca_tcga_sequenced", 
                                 geneticProfile = "brca_tcga_mutations",
                                 genes = choose_genes
        )
        mut_df <- apply(mut_df,2,as.factor)
        mut_df[mut_df == "NaN"] = ""
        mut_df[is.na(mut_df)] = ""
        mut_df[mut_df != ''] = "MUT" 
        ## get copy number data  
        cna <- getProfileData(mycgds, 
                              caseList ="brca_tcga_sequenced", 
                              geneticProfile = "brca_tcga_gistic", 
                              genes = choose_genes)
        rn=rownames(cna)
        cna <- apply(cna,2,function(x)
          as.character(factor(x, levels = c(-2:2), 
      
                              labels = c("HOMDEL", "HETLOSS", "DIPLOID", "GAIN", "AMP"))))
      
        cna[is.na(cna)] = ""
        cna[cna == 'DIPLOID'] = ""
        rownames(cna)=rn
        # Get clinical data for the case list
        myclinicaldata = getClinicalData(mycgds,mycaselist)
        save(expr,myclinicaldata,cna,mut_df,file='survival_input.Rdata')
      
      
      }
      
      load(file='survival_input.Rdata')
      DT::datatable(expr)

      上述代码取决于网速,我已经下载整理好了:survival_input.Rdata 数据,避免每次重复这个教程重新下载的尴尬

      DT::datatable(myclinicaldata,
                        extensions = 'FixedColumns',
                        options = list(
                          #dom = 't',
                          scrollX = TRUE,
                          fixedColumns = TRUE
                        ))## Warning in instance$preRenderHook(instance): It seems your data is too
      ## big for client-side DataTables. You may consider server-side processing:
      ## http://rstudio.github.io/DT/server.html

      可以看到所谓的表达矩阵就是每个基因在各个样本的表达量,只不过是要注意单位,可以是RPKM,TPM等。

      生存分析所需的临床信息一般是生存时间和是否死亡的状态信息。

      这两个数据,也可以自己从其它途径(比如很多人喜欢excel表格)整理,然后读入到R语言里面再来做生存分析。

      临床资料解读

      临床信息是非常复杂的,但是我们并不需要理解那么多,下面是我挑出的一些比较重要的,需要彻底理解:

      choose_columns=c('AJCC_METASTASIS_PATHOLOGIC_PM','AJCC_NODES_PATHOLOGIC_PN','AJCC_PATHOLOGIC_TUMOR_STAGE','AJCC_TUMOR_PATHOLOGIC_PT',
                       "AGE",'SEX','VITAL_STATUS','OS_STATUS','OS_MONTHS',
                       "DFS_MONTHS","DFS_STATUS")
      choose_clinicaldata=myclinicaldata[,choose_columns]
      colnames(choose_clinicaldata)[1:4]=c('M','N','STAGE','T')
      
      str(choose_clinicaldata,  no.list = T, vec.len = 2)
      ## 'data.frame':    1100 obs. of  11 variables:
      ##  $ M           : chr  "M0" "MX" ...
      ##  $ N           : chr  "N0" "N3" ...
      ##  $ STAGE       : chr  "Stage IIA" "Stage IIIC" ...
      ##  $ T           : chr  "T2" "T3" ...
      ##  $ AGE         : int  62 59 70 42 69 ...
      ##  $ SEX         : chr  "Female" "Female" ...
      ##  $ VITAL_STATUS: chr  "Alive" "Alive" ...
      ##  $ OS_STATUS   : chr  "LIVING" "LIVING" ...
      ##  $ OS_MONTHS   : num  10.3 26 ...
      ##  $ DFS_MONTHS  : num  10.3 26 ...
      ##  $ DFS_STATUS  : chr  "DiseaseFree" "DiseaseFree" ...

      虽然上面我挑出的临床信息还有很多,但是我们只需要用到OS_MONTHS和OS_STATUS就可以来估计KM生存函数,画出最简单的生存曲线!

      • 无病生存期(Disease-free survival,DFS)的定义是指从随机化开始至疾病复发或由于疾病进展导致患者死亡的时间。该指标也常作 为抗肿瘤药物III期临床试验的主要终点。某些情况下,DFS与OS相比,作为终点比较难以记录,因为它要求认真随访,及时发现疾病复发,而且肿瘤患者的 死亡原因也很难确定。肿瘤患者常有合并症(如,心血管病),这些合并症可能会干扰对DFS的判断。并且,肿瘤患者常死于医院外,不能常规进行尸检。
      • 总生存期(Overall survival,OS) 的定义是指从随机化开始至因任何原因引起死亡的时间。该指标常常被认为是肿瘤临床试验中最佳的疗效终点。如果在生存期上有小幅度的提高,可以认为是有意义的临床受益证据。作为一个终点,生存期应每天进行评价,可通过在住院就诊时,通过与患者直接接触或者通过电话与患者交谈,这些相对比较容易 记录。确认死亡的日期通常几乎没有困难,并且死亡的时间有其独立的因果关系。当记录至死亡之前的失访患者,通常截止到最后一次有记录的、与患者接触的时间。

      事件时间分布如下:死亡(DECEASED)和存活(LIVING)

      dat=choose_clinicaldata[choose_clinicaldata$OS_MONTHS > 0,]
      table(dat$OS_STATUS)
      ## 
      ## DECEASED   LIVING 
      ##      154      930
      
      attach(dat)
      
      ## The following object is masked from package:base:
      ## 
      ##     T
      
      library(ggplot2)
      ggplot(dat,       
             aes(x = OS_MONTHS, group = OS_STATUS,colour = OS_STATUS,           
                 fill = OS_STATUS
      )) + geom_density(alpha = 0.5)
      detach(dat)

      简单的KM生存分析

      library(survival)
      dat=choose_clinicaldata[choose_clinicaldata$OS_MONTHS > 0,]
      table(dat$OS_STATUS)
      
      ## 
      ## DECEASED   LIVING 
      ##      154      930
      
      attach(dat)
      
      ## The following object is masked from package:base:
      ## 
      ##     T
      
      ## 估计KM生存曲线
      my.surv <- Surv(OS_MONTHS,OS_STATUS=='DECEASED')
      ## The status indicator, normally 0=alive, 1=dead. 
      ## Other choices are TRUE/FALSE (TRUE = death) or 1/2 (2=death). 
      kmfit1 <- survfit(my.surv~1)
      # summary(kmfit1)
      plot(kmfit1)
      detach(dat)

      上面的生存分析并没有指定样本如何进行分组,是最简单版本的生存分析了。

      我们很容易看到,随着感癌症的时间延长,病人的死亡率到了一定程度就不增加了,有些病人熬过了这个癌症,或者说,到我们拿到数据为止,他们还没有死亡!

      有分类的KM生存分析

      首先根据性别分组,看看不同性别的乳腺癌患者的生存曲线

      library(survival)
      dat=choose_clinicaldata[choose_clinicaldata$OS_MONTHS > 0,]
      attach(dat)
      ## The following object is masked from package:base:
      ## 
      ##     T
      table(SEX,OS_STATUS)
      ##         OS_STATUS
      ## SEX      DECEASED LIVING
      ##   Female      153    919
      ##   Male          1     11
      my.surv <- Surv(OS_MONTHS,OS_STATUS=='DECEASED')
      kmfit2 <- survfit(my.surv~SEX)
      plot(kmfit2,col = rainbow(length(unique(SEX))))

      这个其实没啥必要,因为乳腺癌本身在女性发病率远高于男性,尤其是TCGA里面,就12个男性样本

      # 看这个因子的不同水平是否有显著差异,其中默认用是的logrank test 方法。
      survdiff(my.surv~(SEX))
      ## Call:
      ## survdiff(formula = my.surv ~ (SEX))
      ## 
      ## n=1084, 3 observations deleted due to missingness.
      ## 
      ##               N Observed Expected (O-E)^2/E (O-E)^2/V
      ## SEX=Female 1072      153    152.8  0.000262    0.0337
      ## SEX=Male     12        1      1.2  0.033305    0.0337
      ## 
      ##  Chisq= 0  on 1 degrees of freedom, p= 0.854
      # 来检测自己感兴趣的因子是否受其它因子(age,gender等等)的影响。
      # 用strata来控制协变量的影响, 比如年龄等因素
      survdiff(my.surv~SEX+strata(AGE))
      ## Call:
      ## survdiff(formula = my.surv ~ SEX + strata(AGE))
      ## 
      ## n=1084, 3 observations deleted due to missingness.
      ## 
      ##               N Observed Expected (O-E)^2/E (O-E)^2/V
      ## SEX=Female 1072      153  153.181  0.000214    0.0437
      ## SEX=Male     12        1    0.819  0.040052    0.0437
      ## 
      ##  Chisq= 0  on 1 degrees of freedom, p= 0.834
      detach(dat)

      这个分析意义真不大,因为乳腺癌女性患者数量要远超过男性患者,但是有非常多有意义的分类呀!

      然后根据癌症的stage来进行分类

      library(survival)
      dat=choose_clinicaldata[choose_clinicaldata$OS_MONTHS > 0,]
      attach(dat)
      ## The following object is masked from package:base:
      ## 
      ##     T
      table(STAGE,OS_STATUS)
      ##             OS_STATUS
      ## STAGE        DECEASED LIVING
      ##                     4      7
      ##   Stage I          13     77
      ##   Stage IA          3     82
      ##   Stage IB          0      5
      ##   Stage II          0      6
      ##   Stage IIA        34    319
      ##   Stage IIB        34    221
      ##   Stage III         2      0
      ##   Stage IIIA       25    129
      ##   Stage IIIB        8     17
      ##   Stage IIIC        9     57
      ##   Stage IV         15      5
      ##   Stage X           7      5
      my.surv <- Surv(OS_MONTHS,OS_STATUS=='DECEASED')
      kmfit2 <- survfit(my.surv~STAGE)
      plot(kmfit2,col = rainbow(length(unique(STAGE))))
      # 检验显著性
      survdiff(my.surv~(STAGE))
      ## Call:
      ## survdiff(formula = my.surv ~ (STAGE))
      ## 
      ## n=1084, 3 observations deleted due to missingness.
      ## 
      ##                    N Observed Expected (O-E)^2/E (O-E)^2/V
      ## STAGE=            11        4    2.890     0.426     0.443
      ## STAGE=Stage I     90       13   20.141     2.532     2.949
      ## STAGE=Stage IA    85        3   11.594     6.370     7.015
      ## STAGE=Stage IB     5        0    0.825     0.825     0.831
      ## STAGE=Stage II     6        0    1.246     1.246     1.266
      ## STAGE=Stage IIA  353       34   45.149     2.753     3.921
      ## STAGE=Stage IIB  255       34   37.445     0.317     0.420
      ## STAGE=Stage III    2        2    0.209    15.375    15.415
      ## STAGE=Stage IIIA 154       25   20.460     1.007     1.169
      ## STAGE=Stage IIIB  25        8    3.874     4.393     4.598
      ## STAGE=Stage IIIC  66        9    4.997     3.208     3.360
      ## STAGE=Stage IV    20       15    2.389    66.572    67.799
      ## STAGE=Stage X     12        7    2.781     6.403     6.632
      ## 
      ##  Chisq= 112  on 12 degrees of freedom, p= 0
      # 用strata来控制协变量的影响, 比如年龄等因素
      survdiff(my.surv~STAGE+strata(AGE))
      ## Call:
      ## survdiff(formula = my.surv ~ STAGE + strata(AGE))
      ## 
      ## n=1084, 3 observations deleted due to missingness.
      ## 
      ##                    N Observed Expected (O-E)^2/E (O-E)^2/V
      ## STAGE=            11        4    3.943  8.19e-04   0.00214
      ## STAGE=Stage I     90       13   21.507  3.36e+00   6.59955
      ## STAGE=Stage IA    85        3   10.870  5.70e+00   7.79916
      ## STAGE=Stage IB     5        0    0.395  3.95e-01   0.48088
      ## STAGE=Stage II     6        0    1.146  1.15e+00   1.56141
      ## STAGE=Stage IIA  353       34   42.376  1.66e+00   2.93671
      ## STAGE=Stage IIB  255       34   36.959  2.37e-01   0.50934
      ## STAGE=Stage III    2        2    0.545  3.88e+00   4.13017
      ## STAGE=Stage IIIA 154       25   20.534  9.71e-01   1.53489
      ## STAGE=Stage IIIB  25        8    2.585  1.13e+01  14.06707
      ## STAGE=Stage IIIC  66        9    5.054  3.08e+00   3.62911
      ## STAGE=Stage IV    20       15    4.698  2.26e+01  69.87385
      ## STAGE=Stage X     12        7    3.389  3.85e+00   6.94665
      ## 
      ##  Chisq= 117  on 12 degrees of freedom, p= 0
      detach(dat)

      可以看到癌症诊断里面的stage种类太多,生存分析曲线可视化效果不好。后面我们会介绍一个比较好的生存分析结果可视化包。

      根据基因表达量对样本进行分组做生存分析

      首先探索一下基因的表达量情况,这里选取的是BRCA1,BRCA2这两个基因。

      library(survival)
      dat=cbind(choose_clinicaldata[,c('OS_STATUS','OS_MONTHS')],expr[rownames(choose_clinicaldata),])
      dat=dat[dat$OS_MONTHS > 0,]
      dat=dat[!is.na(dat$OS_STATUS),]
      dat$OS_STATUS=as.character(dat$OS_STATUS)
      # ggplot(dat,aes(x=OS_STATUS,y=BRCA1))+geom_boxplot()
      library(ggpubr)
      ## Loading required package: magrittr
      p <- ggboxplot(dat, x="OS_STATUS", y="BRCA1", color = "OS_STATUS",
                     palette = "jco", add = "jitter")
      p+stat_compare_means(method = "t.test") 
      p <- ggboxplot(dat, x="OS_STATUS", y="BRCA2", color = "OS_STATUS",
                     palette = "jco", add = "jitter")
      p+stat_compare_means(method = "t.test")

      有表达差异不一定代表有生存分析的显著性,下面分别进行生存分析检验。

      dat$brca1_group=ifelse(dat$BRCA1 > median(dat$BRCA1),'high','low')
      dat$brca2_group=ifelse(dat$BRCA2 > median(dat$BRCA2),'high','low')
      attach(dat)
      table(brca1_group,brca2_group)
      ##            brca2_group
      ## brca1_group high low
      ##        high  375 167
      ##        low   167 375
      library(survival)
      my.surv <- Surv(OS_MONTHS,OS_STATUS=='DECEASED')
      kmfit1 <- survfit(my.surv~brca1_group,data = dat) 
      plot(kmfit1,col = rainbow( 2 ))
      kmfit2 <- survfit(my.surv~brca2_group,data = dat) 
      plot(kmfit2,col = rainbow( 2 ))
      library("survminer")
      ggsurvplot(kmfit1,conf.int =F, pval = T,risk.table =T, ncensor.plot = TRUE)
      ggsurvplot(kmfit2,conf.int =F, pval = T,risk.table =T, ncensor.plot = TRUE)

      可以看到这个survminer包对生存分析可视化的效果很赞,之所以可以显示P值,是因为我们的survfit函数已经做了检验,返回的kmfit这样的对象里面本身就含有非常丰富的信息,大家可以自行摸索。

      当然,即使某个基因在KM生存分析显示有显著性的差异,也还需要排除其它因素的影响,而不仅仅是该基因表达的高低决定了生存的好坏。所以还需要进行下面的 cox生存分析

      cox生存分析

      Cox比例风险回归模型(Cox’s proportional hazards regression model),简称Cox回归模型。该模型由英国统计学家D.R.Cox于1972年提出,主要用于肿瘤和其它慢性病的预后分析,也可用于队列研究的病因探索。

      参考:

      http://www.sthda.com/english/wiki/cox-proportional-hazards-model

      https://www.r-bloggers.com/cox-proportional-hazards-model/

      library(survival)
      library("survminer")
      
      dat=choose_clinicaldata[choose_clinicaldata$OS_MONTHS > 0,]
      dat=dat[!is.na(dat$OS_STATUS),]
      dat$OS_STATUS=as.character(dat$OS_STATUS)
      
      str(dat,  no.list = T, vec.len = 2)
      ## 'data.frame':    1084 obs. of  11 variables:
      ##  $ M           : chr  "M0" "MX" ...
      ##  $ N           : chr  "N0" "N3" ...
      ##  $ STAGE       : chr  "Stage IIA" "Stage IIIC" ...
      ##  $ T           : chr  "T2" "T3" ...
      ##  $ AGE         : int  62 59 70 42 69 ...
      ##  $ SEX         : chr  "Female" "Female" ...
      ##  $ VITAL_STATUS: chr  "Alive" "Alive" ...
      ##  $ OS_STATUS   : chr  "LIVING" "LIVING" ...
      ##  $ OS_MONTHS   : num  10.3 26 ...
      ##  $ DFS_MONTHS  : num  10.3 26 ...
      ##  $ DFS_STATUS  : chr  "DiseaseFree" "DiseaseFree" ...
      attach(dat)
      ## The following objects are masked from dat (pos = 4):
      ## 
      ##     OS_MONTHS, OS_STATUS
      ## The following object is masked from package:base:
      ## 
      ##     T
      my.surv <- Surv( dat$OS_MONTHS,dat$OS_STATUS=='DECEASED')
      ## 可以针对多个变量分别批量进行cox分析,或者直接把多个变量一起放在cox函数里面。
      m=coxph(my.surv ~ AGE + SEX + N + T + M + STAGE, data =  dat)
      ## Warning in fitter(X, Y, strats, offset, init, control, weights = weights, :
      ## Loglik converged before variable 5,9,18,22,25,32,35,36 ; beta may be
      ## infinite.
      ggsurvplot(survfit(m, data =  dat), color = "#2E9FDF",
                 ggtheme = theme_minimal())
      ## Warning: Now, to change color palette, use the argument palette= '#2E9FDF'
      ## instead of color = '#2E9FDF'
      beta <- coef(m)
      se <- sqrt(diag(vcov(m)))
      HR <- exp(beta)
      HRse <- HR * se
      
      detach(dat)

      某基因突变与否也可以用来分组

      library(survival)
      length(intersect(rownames(choose_clinicaldata),rownames(mut_df)))
      ## [1] 979
      dat=cbind(choose_clinicaldata[rownames(mut_df),c('OS_STATUS','OS_MONTHS')],mut_df)
      dat=dat[dat$OS_MONTHS > 0,]
      dat=dat[!is.na(dat$OS_STATUS),]
      dat$OS_STATUS=as.character(dat$OS_STATUS) 
      
      attach(dat)
      ## The following objects are masked from dat (pos = 4):
      ## 
      ##     BRCA1, BRCA2, OS_MONTHS, OS_STATUS
      table(BRCA1,BRCA2)
      ##      BRCA2
      ## BRCA1     MUT
      ##       939  14
      ##   MUT  12   0
      library(survival)
      my.surv <- Surv(OS_MONTHS,OS_STATUS=='DECEASED')
      kmfit1 <- survfit(my.surv~BRCA1,data = dat)  
      kmfit2 <- survfit(my.surv~BRCA2,data = dat)  
      library("survminer") 
      ggsurvplot(kmfit1,conf.int =F, pval = T,risk.table =T, ncensor.plot = TRUE)
      ggsurvplot(kmfit2,conf.int =F, pval = T,risk.table =T, ncensor.plot = TRUE)

      基因的拷贝数也可以进行分组

      library(survival)
      length(intersect(rownames(choose_clinicaldata),rownames(cna)))
      ## [1] 979
      dat=cbind(choose_clinicaldata[rownames(cna),c('OS_STATUS','OS_MONTHS')],cna)
      dat=dat[dat$OS_MONTHS > 0,]
      dat=dat[!is.na(dat$OS_STATUS),]
      dat$OS_STATUS=as.character(dat$OS_STATUS) 
      
      attach(dat)
      ## The following objects are masked from dat (pos = 3):
      ## 
      ##     BRCA1, BRCA2, OS_MONTHS, OS_STATUS
      ## The following objects are masked from dat (pos = 5):
      ## 
      ##     BRCA1, BRCA2, OS_MONTHS, OS_STATUS
      table(BRCA1,BRCA2)
      ##          BRCA2
      ## BRCA1         AMP GAIN HETLOSS HOMDEL
      ##           288   2   21     118      8
      ##   AMP       5   0    0       7      2
      ##   GAIN     59   4   30      84      3
      ##   HETLOSS 110   6   43     163      4
      ##   HOMDEL    2   0    1       4      1
      library(survival)
      my.surv <- Surv(OS_MONTHS,OS_STATUS=='DECEASED')
      kmfit1 <- survfit(my.surv~BRCA1,data = dat)  
      kmfit2 <- survfit(my.surv~BRCA2,data = dat)  
      library("survminer") 
      ggsurvplot(kmfit1,conf.int =F, pval = T,risk.table =T, ncensor.plot = TRUE)
      ggsurvplot(kmfit2,conf.int =F, pval = T,risk.table =T, ncensor.plot = TRUE)

      批量进行分组并且做生存分析

      这个代码非常经典,所以我可以放在了 生信技能树»生信技能树›互动作业›脚本能力实践›生信人必练的200个数据处理任务 上面,希望所有学习生物信息学的朋友都掌握:http://www.biotrainee.com/thread-929-1-1.html

      (阅读原文可以直达)

      ## 首先模拟 50个病人的 200个基因的表达数据。
      ## 50 patients and 200 genes 
      dat=matrix(rnorm(10000,mean=8,sd=4),nrow = 200)
      rownames(dat)=paste0('gene_',1:nrow(dat))
      colnames(dat)=paste0('sample_',1:ncol(dat))
      os_years=abs(floor(rnorm(ncol(dat),mean = 50,sd=20)))
      os_status=sample(rep(c('LIVING','DECEASED'),25))
      
      library(survival)
      my.surv <- Surv( os_years,os_status=='DECEASED')
      ## The status indicator, normally 0=alive, 1=dead. Other choices are TRUE/FALSE (TRUE = death) or 1/2 (2=death). 
      ## And most of the time we just care about the time od DECEASED . 
      
      fit.KM=survfit(my.surv~1)
      fit.KM
      ## Call: survfit(formula = my.surv ~ 1)
      ## 
      ##       n  events  median 0.95LCL 0.95UCL 
      ##      50      25      64      56      NA
      plot(fit.KM)
      log_rank_p <- apply(dat, 1, function(values1){
        group=ifelse(values1>median(values1),'high','low')
        kmfit2 <- survfit(my.surv~group)
        #plot(kmfit2)
        data.survdiff=survdiff(my.surv~group)
        p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
      })
      names(log_rank_p[log_rank_p<0.05])
      ## [1] "gene_63"  "gene_100" "gene_102" "gene_125" "gene_146" "gene_149"
      ## [7] "gene_150" "gene_193"
      gender= ifelse(rnorm(ncol(dat))>1,'male','female')
      age=abs(floor(rnorm(ncol(dat),mean = 50,sd=20)))
      ## gender and age are similar with group(by gene expression)
      
      cox_results <- apply(dat , 1, function(values1){
        group=ifelse(values1>median(values1),'high','low')
        survival_dat <- data.frame(group=group,gender=gender,age=age,stringsAsFactors = F)
        m=coxph(my.surv ~ age + gender + group, data =  survival_dat)
      
        beta <- coef(m)
        se <- sqrt(diag(vcov(m)))
        HR <- exp(beta)
        HRse <- HR * se
      
        #summary(m)
        tmp <- round(cbind(coef = beta, se = se, z = beta/se, p = 1 - pchisq((beta/se)^2, 1),
                           HR = HR, HRse = HRse,
                           HRz = (HR - 1) / HRse, HRp = 1 - pchisq(((HR - 1)/HRse)^2, 1),
                           HRCILL = exp(beta - qnorm(.975, 0, 1) * se),
                           HRCIUL = exp(beta + qnorm(.975, 0, 1) * se)), 3)
        return(tmp['grouplow',])
      
      })
      DT::datatable(cox_results ,
                        extensions = 'FixedColumns',
                        options = list(
                          #dom = 't',
                          scrollX = TRUE,
                          fixedColumns = TRUE
                        ))
      ## 有统计学显著性的如下:
      cox_results[,cox_results[4,]<0.05]
      ##        gene_26 gene_63 gene_94 gene_100 gene_102 gene_103 gene_146
      ## coef     0.853   1.109  -0.903    0.935   -0.995    0.890   -1.372
      ## se       0.421   0.437   0.431    0.426    0.464    0.429    0.478
      ## z        2.025   2.534  -2.096    2.196   -2.146    2.073   -2.869
      ## p        0.043   0.011   0.036    0.028    0.032    0.038    0.004
      ## HR       2.346   3.030   0.405    2.548    0.370    2.435    0.254
      ## HRse     0.988   1.326   0.175    1.086    0.171    1.046    0.121
      ## HRz      1.363   1.531  -3.406    1.426   -3.676    1.373   -6.154
      ## HRp      0.173   0.126   0.001    0.154    0.000    0.170    0.000
      ## HRCILL   1.028   1.285   0.174    1.106    0.149    1.050    0.099
      ## HRCIUL   5.354   7.142   0.943    5.874    0.918    5.649    0.647
      ##        gene_149 gene_150 gene_155 gene_193
      ## coef      1.011    1.004   -0.924   -1.072
      ## se        0.442    0.479    0.440    0.440
      ## z         2.290    2.096   -2.101   -2.437
      ## p         0.022    0.036    0.036    0.015
      ## HR        2.749    2.728    0.397    0.342
      ## HRse      1.214    1.306    0.175    0.151
      ## HRz       1.441    1.323   -3.456   -4.368
      ## HRp       0.150    0.186    0.001    0.000
      ## HRCILL    1.157    1.068    0.168    0.145
      ## HRCIUL    6.531    6.973    0.940    0.811

      请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!

      • 分享:
      weinfoauthor
      weinfoauthor

      1233

      上一篇文章

      (未测试)R语言运行相关性分析
      2020年2月6日

      下一篇文章

      环状RNA 4分文章路线及代码-GEO芯片数据
      2020年2月7日

      你可能也喜欢

      2-1675088548
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      30 1月, 2023
      9-1675131201
      如何快速批量修改 Git 提交记录中的用户信息
      26 1月, 2023
      8-1678501786
      肿瘤细胞通过改变CD8+ T细胞中的丙酮酸利用和琥珀酸信号来调控抗肿瘤免疫应答。
      7 12月, 2022

      留言 取消回复

      要发表评论,您必须先登录。

      搜索

      分类

      • R语言
      • TCGA数据挖掘
      • 单细胞RNA-seq测序
      • 在线会议直播预告与回放
      • 数据分析那些事儿分类
      • 未分类
      • 生信星球
      • 老俊俊的生信笔记

      投稿培训

      免费

      alphafold2培训

      免费

      群晖配置培训

      免费

      最新博文

      Nature | 单细胞技术揭示衰老细胞与肌肉再生
      301月2023
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      301月2023
      如何快速批量修改 Git 提交记录中的用户信息
      261月2023
      logo-eduma-the-best-lms-wordpress-theme

      (00) 123 456 789

      weinfoadmin@weinformatics.cn

      恒诺新知

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      链接

      • 课程
      • 事件
      • 展示
      • 问答

      支持

      • 文档
      • 论坛
      • 语言包
      • 发行状态

      推荐

      • iHub汉语代码托管
      • iLAB耗材管理
      • WooCommerce
      • 丁香园论坛

      weinformatics 即 恒诺新知。ICP备案号:粤ICP备19129767号

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      要成为一名讲师吗?

      加入数以千计的演讲者获得100%课时费!

      现在开始

      用你的站点账户登录

      忘记密码?

      还不是会员? 现在注册

      注册新帐户

      已经拥有注册账户? 现在登录

      close
      会员购买 你还没有登录,请先登录
      • ¥99 VIP-1个月
      • ¥199 VIP-半年
      • ¥299 VIP-1年
      在线支付 激活码

      立即支付
      支付宝
      微信支付
      请使用 支付宝 或 微信 扫码支付
      登录
      注册|忘记密码?