准备
https://occdn.limour.top/1895.html
https://occdn.limour.top/1945.html
https://occdn.limour.top/1820.html
来源
https://mp.weixin.qq.com/s/ykQm3OygcfPgC8P0q4cwZA
https://blog.csdn.net/qq_19600291/article/details/120465166
https://cloud.tencent.com/developer/article/1708598
整合数据
1
2
3
4
5
6
7
8
9
10
11
12
| 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
|
生存分析
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
| 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
|