1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
| library("survival")
library("survminer")
f_DEG_coxph_oneGene <- function(trainSet, geneN){
if(!(geneN %in% colnames(trainSet$data))){
return(data.frame(row.names = geneN, mean=NA, lower=NA, upper=NA, Pvalue=1, VarName=geneN))
}
df <- cbind(trainSet$data[geneN],trainSet$meta[c('pfs_status', 'pfs_time')])
coxmf <- paste0("Surv(pfs_time, pfs_status==1)~", '`',geneN,'`')
res.cox <- coxph(formula(coxmf), data = df)
tmp_res <- summary(res.cox)
res <- as.data.frame(tmp_res$conf.int)
res <- res[,c(1,3,4)]
colnames(res) <- c('mean', 'lower', 'upper')
res[['Pvalue']] <- tmp_res$coefficients[,'Pr(>z)']
res[['VarName']] <- rownames(res)
res
}
f_DEG_coxph_GeneList <- function(dat, geneNs){
if(length(geneNs) == 1){
return(f_DEG_coxph_oneGene(dat,geneNs))
}
y <- Surv(round(dat$meta$pfs_time,2), dat$meta$pfs_status == 1)
y <- data.matrix(y)
x <- dat$data[,geneNs]
x <- as.matrix(x)
df <- cbind(x,y)
df <- as.data.frame(df)
coxmf <- paste0("Surv(time, status==1)~", paste(paste('`',geneNs,'`',sep = ''), collapse = '+'))
res.cox <- coxph(formula(coxmf), data = df)
tmp_res <- summary(res.cox)
res <- as.data.frame(tmp_res$conf.int)
res <- res[,c(1,3,4)]
colnames(res) <- c('mean', 'lower', 'upper')
res[['Pvalue']] <- tmp_res$coefficients[,'Pr(>z)']
res[['VarName']] <- rownames(res)
res
}
f_DEG_coxph <- function(trainSet, geneList, FPv=0.05, strict=T){
res <- NULL
for(gene in geneList){
res <- rbind(res,f_DEG_coxph_oneGene(trainSet, gene))
}
if(!strict){
return(res)
}
res <- subset(res, Pvalue<FPv)$VarName
if(length(res) == 0){
return(data.frame(mean=NA, lower=NA, upper=NA, Pvalue=1, VarName=NA))
}
res <- f_DEG_coxph_GeneList(trainSet, res)
res[order(res$Pvalue),]
}
f_DEG_coxph_Sets <- function(trainSets, geneList, FPv=0.05, strict=T){
res <- list()
for(Name in names(trainSets)){
res[[Name]] <- f_DEG_coxph(trainSets[[Name]], geneList, FPv = FPv, strict = strict)
}
res
}
x <- readRDS('../A_ref_A_fiig.2_A/hub.rds')
train_sets <- readRDS('../../../COXPH/train_sets.rds')
tmp <- f_DEG_coxph_Sets(train_sets, x)
|