生存分析的基本流程

发布于 2022-07-03  125 次阅读


准备

来源

https://mp.weixin.qq.com/s/ykQm3OygcfPgC8P0q4cwZA

https://blog.csdn.net/qq_19600291/article/details/120465166

https://cloud.tencent.com/developer/article/1708598

整合数据

clinical <- read.csv('../tcga/PRAD/prad_survival.csv')
group <- read.csv('Est.prop.weighted.csv')
mergeID <- intersect(clinical$bcr_patient_barcode, group$X)
clinical <- clinical[clinical$bcr_patient_barcode %in% mergeID,]
group <- group[group$X %in% mergeID,]
rownames(clinical) <- clinical$bcr_patient_barcode
rownames(group) <- group$X
df <- clinical[mergeID, c('os_time', 'os_status')]
df <- cbind(df, group[mergeID, 'iCAF'])
colnames(df) <- c('os_time','os_status','group_o')
df[['group']] <- ifelse(df$group_o < median(df$group_o), 0, 1)
df

生存分析

library(survival)
library(survminer)
library(ggpubr)
require(survRM2)
f_surv <- function(df, isplot=T, timeN='os_time', statusN='os_status', groupN='group', tau=5){
    res <- list()
    my.surv <<- Surv(df[[timeN]]/365, df[[statusN]]==0)
    my.surv.f <<- paste0("my.surv~", groupN)
    my.surv.df <<- df
    kmfit <- surv_fit(formula(my.surv.f), data = my.surv.df)
    sdiff <- survdiff(formula(my.surv.f), data = my.surv.df)
    res$pv <- 1 - pchisq(sdiff$chisq, length(sdiff$n) - 1)
    if(isplot){
        res$plot <- ggsurvplot(kmfit, conf.int =F, pval = T, risk.table =T, ncensor.plot = TRUE)
    }
    res$RMS <- rmst2(df[[timeN]]/365, df[[statusN]]==0, df[[groupN]], tau=tau)
    res
}
df <- readRDS('test.rds')
test <- f_surv(df)
test$RMS

一枚爱好探索的医学生