# cox.threshold3 help get a comprimized cut-off between mutiple design object
# expr.matrix# expression object
# design.list# a list of design object.The names of them must be provided
# model # model from List.ModelFromOptimized2
# event.status#colname of status
# event.time#colname of time
# event.lower # lower limit of time
# interval.n # the cut-off selected parameter.Default is 1000
# cutoff.pval # cut-off of p value
#' @export
cox.threshold3 <- function(expr.matrix,
design.list,
model,
event.status=status,
event.time=time,
event.lower = event.lower,
interval.n = 1000,
cutoff.pval = 0.05){
## 加载必要的包
nd <- c("survival","broom","plyr")
Plus.library(nd)
## 获得risk.score矩阵和cutoff节点
print("Getting score matrixes and cutoff limits...")
L1 <- NULL
for(i in 1:length(design.list)){ #i=1
design.i <- design.list[[i]]
## 对齐
expr.matrix.i <- expr.matrix[,rownames(design.i)]
##获得患者的risk score
model.i <- model #仅1个model
genes <- as.character(model.i$ENSEMBL[is.na(model.i$SYMBOL)==F])
lg1 <- rownames(expr.matrix.i) %in% genes
expr.matrix.i1 <- subset(expr.matrix.i,subset = lg1)
if(length(genes) >= 2){expr.matrix.i1 <- expr.matrix.i1[genes,]} #基因一定要对齐!
#expr.matrix.i1 <- expr.matrix.i[genes,]
cor <- as.numeric(as.character(model.i$Cor[is.na(model.i$SYMBOL)==F]))
expr.matrix.i1 <- cbind(cor,expr.matrix.i1)
expr.matrix.i2 <- apply(expr.matrix.i1,1,function(x) x[1]*x[2:length(x)])
expr.matrix.i3 <- apply(expr.matrix.i2,1,sum)
expr.matrix.i3 <- as.data.frame(expr.matrix.i3);
colnames(expr.matrix.i3)[1] <-"risk.score"
##组合数据
prognosis <- c(event.time,event.status)
design.i1 <- design.i[,prognosis]
design.i1 <- cbind(design.i1,expr.matrix.i3)
colnames(design.i1) <- c("time","status","risk.score")
design.i1 <- design.i1[!is.na(design.i1[,1]),] #去NA值
lg1 <- design.i1[,1] > event.lower
design.i1 <- subset(design.i1,subset = lg1)
##cut-off节点
max.i <- max(design.i1$risk.score)
min.i <- min(design.i1$risk.score)
##
l <- list(max=max.i,min = min.i,design.i = design.i1)
L1 <- c(L1,list(l))
names(L1)[i] <- names(design.list)[i]
}
## 获得cutoff向量
min1 <- NULL;max1 <- NULL
for(i in 1:length(L1)){
min1 <- c(min1,L1[[i]][["min"]])
max1 <- c(max1,L1[[i]][["max"]])
}
min2 <- max(min1);max2 <- min(max1)
interval <- max2 - min2
int1 <- interval/interval.n
cutoff <- c(min2,rep(0,interval.n))#length(cutoff)
for(j in 1:interval.n){cutoff[j+1] <- cutoff[j]+int1}
cutoff <- cutoff[c(-1,-(interval.n+1))]
##计算某个cutoff值的p值
P <- NULL
print("Getting p statistics for every design object...")
for(i in 1:length(L1)){ # i=1
design.i2 <- L1[[i]][["design.i"]]
get1 <- function(cutoff.i){
rs <- ifelse(design.i2$risk.score >= cutoff.i,1,0);
design.i3 <- cbind(design.i2,rs)
sdf <- survdiff(Surv(time, status) ~ rs, data = design.i3)
p.test <- glance(sdf)
return(p.test)
}
cutoff.a <- as.matrix(cutoff)
p.test <- adply(cutoff.a,1,get1,.id = "cutoff")
p.test$cutoff <- cutoff
p.test <- subset(p.test,select = c("cutoff","statistic","p.value"))
p.test <- p.test[order(p.test$p.value,decreasing = F),]
print(paste0(names(design.list)[i],": p statistics has been gotten!"))
p.test$group <- names(design.list)[i]
P <- rbind(P,p.test)
}
## 获取系列cut-off值
get2 <- function(P1){
p.val <- P1[,"p.value"]
lg1 <- all(p.val < cutoff.pval)
d <- data.frame(logic = lg1,mean.p= mean(p.val))
return(d)
}
df1 <- ddply(P,"cutoff",.fun = get2)
## 优化cutoff值的结果整理
cutoff.ok <- df1[df1$logic == T,]
cutoff.ok <- cutoff.ok[order(cutoff.ok$mean.p,decreasing = F),]
## 输出结果
print(paste0("Get ",length(cutoff.ok$cutoff)," optimized cut-off values for design list."))
l <- list(
cutoff = cutoff.ok$cutoff,
data.cutoff = cutoff.ok,
meta.data = P
)
print("All done!")
return(l)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.