搜索
查看: 10758|回复: 9

对有临床信息的表达矩阵批量做生存分析

[复制链接]

633

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2017-2-22 10:18:32 | 显示全部楼层 |阅读模式
R里面实现生存分析非常简单!

用my.surv <- surv(OS_MONTHS,OS_STATUS=='DECEASED')构建生存曲线。
用kmfit2 <- survfit(my.surv~TUMOR_STAGE_2009)来做某一个因子的KM生存曲线。
用 survdiff(my.surv~type, data=dat)来看看这个因子的不同水平是否有显著差异,其中默认用是的logrank test 方法。
用coxph(Surv(time, status) ~ ph.ecog + tt(age), data=lung) 来检测自己感兴趣的因子是否受其它因子(age,gender等等)的影响。


包括KM和cox的,我就不解释原理了,直接上代码:
[AppleScript] 纯文本查看 复制代码
rm(list=ls())
## 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
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])

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',])
  
})
cox_results[,cox_results[4,]<0.05]


这个R代码是可以直接运行的,里面是我模拟的测试数据,需要一步步运行,才能更好的理解!

PS: 里面的os_years应该修改为os_month,因为总的生存期不应该长达几十年,而且因为ag和os_years都是随机生成的,可能会出现很不符合自然科学的现象。但是这个是模拟数据,请大家不用较真。

如果想知道原理,请在本论坛搜索生存分析即可!





上一篇:建议:请放弃GTF文件格式,改用GenePred格式处理转录组注释!
下一篇:hapmapsnp6就存储着一堆hapmap计划包括的snp6芯片数据而已
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

6

主题

36

帖子

465

积分

中级会员

Rank: 3Rank: 3

积分
465
发表于 2017-2-23 15:34:22 | 显示全部楼层
谢楼主,最近一直卡在这里
回复 支持 反对

使用道具 举报

0

主题

2

帖子

83

积分

注册会员

Rank: 2

积分
83
发表于 2017-2-25 02:09:24 | 显示全部楼层
COX这个模型还是加上coxzph检验,筛除不符合等比例检验的或者多家几个因素,做个stepAIC筛选下变量
回复 支持 反对

使用道具 举报

4

主题

24

帖子

184

积分

注册会员

Rank: 2

积分
184
发表于 2017-3-4 14:13:51 | 显示全部楼层
这些东西我曾经看懂过,每次再看就又不懂了,得花好多时间才能再次看懂。被逼无奈,用rpy2写成函数放我的Python package里了。

输入:clin表格,每行是一个病人。第一列"_E"是event status,第二列"_T"是survival time。如果后面还有列的话,会全部加进去做多元CoxPH.          data表格:每行是一个病人。每列是一个基因的数值,比如表达量。
输出:rst表格:每行是一个基因。第一列是”coef”,第二列是“pvalue”。

[Python] 纯文本查看 复制代码
import pandasimport rpy2.robjects as robjects
R = robjects.r
robjects.numpy2ri.activate()
pandas2ri.activate()
def CoxPH(clin,data):
    '''
    Cox (Proportional Hazards) Regression analysis.
    Parameters:
        clin: pandas.DataFrame
            Row name: patient ID.
            Required columns: "_E": event, "_T": time
            Optional columns: eg. age, gender et. al. for multivariate CoxPH regression.
        data: pandas.DataFrame
            Expression values. Row: patient ID, Col: gene names.
    Returns:
        rst: pandas.DataFrame
            CoxPH reuslt. coef, pvalue
    '''
    R.assign('clin',clin)
    strata = '+' + '+'.join(clin.columns[2:]) if clin.shape[1]>2 else "" # Multi-variate CoxPH if more clinical features
    R("library('survival')")
    rst = []
    for name in data.columns:
        try:
            cdata = clin.merge(data.ix[:,[name]],left_index=True,right_index=True)
            R.assign('data',cdata)
            R('coef<- coxph(Surv(X_T, X_E) ~ {0}{1}, data=cdata)'.format(name,strata))
            R('res <-summary(coef)$coefficients[\'{0}\',c(1,5)]'.format(name)) # coef & pvalue
            rst.append(list(R['res']))
        except:
            rst.append((numpy.nan,numpy.nan))
    coef, pval = zip(*rst)
    return pandas.DataFrame({'coef':coef,'pvalue':pval},index=data.columns)




Survdiff是同理,就不列出来了。需要的话就扔进自己的工具包吧。
回复 支持 反对

使用道具 举报

2

主题

17

帖子

145

积分

注册会员

Rank: 2

积分
145
发表于 2017-5-11 11:24:53 | 显示全部楼层
tsznxx 发表于 2017-3-4 14:13
这些东西我曾经看懂过,每次再看就又不懂了,得花好多时间才能再次看懂。被逼无奈,用rpy2写成函数放我的Py ...

请问您的代码可以用于R语言吗?
回复 支持 反对

使用道具 举报

0

主题

1

帖子

31

积分

新手上路

Rank: 1

积分
31
发表于 2017-6-14 20:32:47 | 显示全部楼层
Jimmy老师,您好!我想做某一个基因根据其表达量的median分组,在同一个图上做两条生存曲线,并把log_rank P-value, HR, 95%CI放到注解里面,我不知道怎么写代码,请Jimmy老师不吝赐教~
回复 支持 反对

使用道具 举报

4

主题

20

帖子

395

积分

中级会员

Rank: 3Rank: 3

积分
395
发表于 2017-6-14 22:07:19 | 显示全部楼层
看看这个代码,感觉自己又进步了一些,虽然我离老司机们还是那么遥远
回复 支持 反对

使用道具 举报

2

主题

51

帖子

476

积分

中级会员

Rank: 3Rank: 3

积分
476
发表于 2017-6-29 21:05:54 | 显示全部楼层
之前把生存数据和变量放在一个 table 了,想了半天 survfit 的公式要怎么写,看到 Jimmy 的代码豁然开朗:直接把数据拆开,先生成 surv 对象啊
回复 支持 反对

使用道具 举报

1

主题

5

帖子

54

积分

注册会员

Rank: 2

积分
54
发表于 2019-9-25 20:44:56 | 显示全部楼层
这个做完批量生存分析后只能得到p值,怎么用一些自动的方法筛选出哪些基因是高表达对生存不利的呢
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-11-12 10:12 , Processed in 0.038234 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.