對(duì)有臨床信息的表達(dá)矩陣批量做生存分析

ypgao 2017-07-30

展開全文

R里面實(shí)現(xiàn)生存分析非常簡(jiǎn)單！

用my.surv <- surv(OS_MONTHS,OS_STATUS=='DECEASED')構(gòu)建生存曲線。
用kmfit2 <- survfit(my.surv~TUMOR_STAGE_2009)來(lái)做某一個(gè)因子的KM生存曲線,。
用 survdiff(my.surv~type, data=dat)來(lái)看看這個(gè)因子的不同水平是否有顯著差異，其中默認(rèn)用是的logrank test 方法,。
用coxph(Surv(time, status) ~ ph.ecog + tt(age), data=lung) 來(lái)檢測(cè)自己感興趣的因子是否受其它因子(age,gender等等)的影響。

包括KM和cox的，我就不解釋原理了,，直接上代碼：

[AppleScript] 純文本查看 復(fù)制代碼

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]

這個(gè)R代碼是可以直接運(yùn)行的，里面是我模擬的測(cè)試數(shù)據(jù),，需要一步步運(yùn)行,，才能更好的理解！

PS: 里面的os_years應(yīng)該修改為os_month,，因?yàn)榭偟纳嫫诓粦?yīng)該長(zhǎng)達(dá)幾十年,，而且因?yàn)閍g和os_years都是隨機(jī)生成的，可能會(huì)出現(xiàn)很不符合自然科學(xué)的現(xiàn)象,。但是這個(gè)是模擬數(shù)據(jù),，請(qǐng)大家不用較真。

如果想知道原理,，請(qǐng)?jiān)诒菊搲阉魃娣治黾纯桑?br>

本站是提供個(gè)人知識(shí)管理的網(wǎng)絡(luò)存儲(chǔ)空間,，所有內(nèi)容均由用戶發(fā)布，不代表本站觀點(diǎn),。請(qǐng)注意甄別內(nèi)容中的聯(lián)系方式,、誘導(dǎo)購(gòu)買等信息，謹(jǐn)防詐騙,。如發(fā)現(xiàn)有害或侵權(quán)內(nèi)容,，請(qǐng)點(diǎn)擊一鍵舉報(bào)。

轉(zhuǎn)藏 分享

QQ空間 QQ好友新浪微博微信

獻(xiàn)花（0） +1

來(lái)自： ypgao > 《數(shù)據(jù)庫(kù)》

舉報(bào)/認(rèn)領(lǐng)