TCGA的28篇教程- 对TCGA数据库的任意癌症中任意基因做生存分析

长期更新列表: 使用R语言的cgdsr包获取TCGA数据(cBioPortal)TCGA的28篇教程- 使用R语言的RTCGA包获取TCGA数据 (离线打包版本)TCGA的28篇教程- 使用R语言的RTCGAToolbox包获取TCGA数据 (FireBrowse portal)TCGA的28篇教程- 批量下载TCGA所有数据 ( UCSC的 XENA)TCGA的28篇教程- 数据下载就到此为止吧 TCGA的28篇教程- 指定癌症查看感兴趣基因的表达量

本教程目录:

  • 首先使用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)

img

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)

img

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个男性样本!img

# 看这个因子的不同水平是否有显著差异,其中默认用是的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))))

img

# 检验显著性
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") 

img

p <- ggboxplot(dat, x="OS_STATUS", y="BRCA2", color = "OS_STATUS",
               palette = "jco", add = "jitter")
p+stat_compare_means(method = "t.test") 

img

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

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 ))

img

kmfit2 <- survfit(my.surv~brca2_group,data = dat) 
plot(kmfit2,col = rainbow( 2 ))

img

library("survminer")
ggsurvplot(kmfit1,conf.int =F, pval = T,risk.table =T, ncensor.plot = TRUE)

img

ggsurvplot(kmfit2,conf.int =F, pval = T,risk.table =T, ncensor.plot = TRUE)

img

可以看到这个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'

img

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)

img

ggsurvplot(kmfit2,conf.int =F, pval = T,risk.table =T, ncensor.plot = TRUE)

img

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

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)

img

ggsurvplot(kmfit2,conf.int =F, pval = T,risk.table =T, ncensor.plot = TRUE)

img

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

这个代码非常经典,所以我可以放在了 生信技能树»生信技能树›互动作业›脚本能力实践›生信人必练的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)

img

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

本文分享自微信公众号 - 生信技能树(biotrainee)

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2018-06-30

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏灯塔大数据

探秘 | Python 求职 Top10 城市,来看看是否有你所在的城市

导读:从智联招聘爬取相关信息后,我们关心的是如何对内容进行分析,获取用用的信息。本次以上篇文章“5分钟掌握智联招聘网站爬取并保存到MongoDB数据库”中爬取的...

35030
来自专栏*坤的Blog

hdu1011

18560
来自专栏阿凯的Excel

统计函数与通配符相爱,不是意外!

首插入音乐,功放党请慎点。 海鸟跟鱼相爱,只是一场意外! 但是统计函数和通配符相爱,却是一种必然! 统计函数何许人也:Sumif、Countif、Average...

32160
来自专栏Albert陈凯

StuQ 大数据工程师技能图谱

https://github.com/TeamStuQ/skill-map StuQ 程序员技能图谱 官网 Web 页面地址:http://skill-map...

427100
来自专栏ATYUN订阅号

赫尔辛基大学AI基础教程:搜索和解决问题(2.1节)

想象一下你在一个外国的城市,在某个地方(比如一家酒店),想用公共交通工具去另一个地方(比如一家不错的餐馆)。你是做什么?如果你会像许多人一样,掏出智能手机,输入...

15660
来自专栏算法+

MTCNN人脸检测 附完整C++代码

Joint Face Detection and Alignment using Multi-task Cascaded Convolutional Neura...

2.7K50

利用基因突变和K均值预测地区种群

这是一篇关于西北基因组中心的Deborah Siegel和华盛顿大学联合Databricks的Denny Lee,就ADAM和Spark基因组变异分析方面的合作...

300100
来自专栏about云

使用Spark MLlib给豆瓣用户推荐电影

问题导读: 1.常用的推荐算法有哪些? 2.推荐系统是什么样的流程? 3.从这个推荐系统我们能学到什么? 推荐算法就是利用用户的一些行为,通过一些数学算法,推测...

96270
来自专栏数据结构与算法

CDQZ 0003:jubeeeeeat

总时间限制: 1000ms 内存限制: 256000kB描述 众所周知,LZF很喜欢打一个叫Jubeat的游戏。这是个音乐游戏,游戏界面是4×4的方阵,会根...

33560
来自专栏崔庆才的专栏

中文分词原理及常用Python中文分词库介绍

99260

扫码关注云+社区

领取腾讯云代金券