# cox生存分析结果也可以火山图可视化

### 首先明白火山图的本质

• X轴必须是0附近的值，比如上下调基因的变化倍数，可以log转换，这样看起来对称
• Y轴必须是统计学相关P值，进行 -log10的转换，这样与你的X轴数值离散开来就有火山喷发的感觉

### 然后认识cox结果

``````# 这里是核心代码
# 会根据每个基因的表达量中位值进行分组，然后对高低表达某基因的两组病人进行cox分析
cox_results <-apply(exprSet , 1 , function(gene){
# gene= as.numeric(exprSet[1,])
group=ifelse(gene>median(gene),'high','low')
survival_dat <- data.frame(group=group,stage=phe\$stage,
stringsAsFactors = F)
m=coxph(mySurv ~ stage+ 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=t(cox_results)
# 把结果保存一下，后续绘制火山图
save(cox_results,file = 'cox_results.Rdata')
``````

``````step00-install-packages.R
step01-getData-from-GDC.R
step01-getData-from-RTCGA.R
step01-getData-from-Xena.R
step01-getData-from-firehose.R
step02-DEG-3-packages.R
step03-batch-logRank.R
step04-batch-coxp.R
step05-lasso.R
step06-coxph-forest.R
step07-risk-score-distribution.R
step08-Random-foreast.R
step09-miRNA-downstream.R
step10-maftools.R
step11-boxplot.R
step12-correlation.R
step13-split-cohort.R
step14-timeROC.R
step15-choose_lncRNA.R
step16-clinical-tables.R
step17-mutation-signatures.R
step17-others.R
step18-SVM.R
``````

### 最后对那两列数据画散点图

``````load(file = 'cox_results.Rdata')
if(F){
need_DEG=as.data.frame(cox_results[,c(4,5)])
colnames(need_DEG)=c('pvalue','log2FoldChange')
need_DEG\$log2FoldChange=log( need_DEG\$log2FoldChange)
boxplot(need_DEG\$log2FoldChange)
fivenum(need_DEG\$log2FoldChange)
need_DEG\$log2FoldChange[abs(need_DEG\$log2FoldChange)>3]=3
logFC_cutoff <- with(need_DEG,mean(abs( log2FoldChange)) + 2*sd(abs( log2FoldChange)) )
# logFC_cutoff=1
logFC_cutoff
need_DEG\$change = as.factor(ifelse(need_DEG\$pvalue < 0.05 & abs(need_DEG\$log2FoldChange) > logFC_cutoff,
ifelse(need_DEG\$log2FoldChange > logFC_cutoff ,'UP','DOWN'),'NOT')
)
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
'\nThe number of up gene is ',nrow(need_DEG[need_DEG\$change =='UP',]) ,
'\nThe number of down gene is ',nrow(need_DEG[need_DEG\$change =='DOWN',])
)
library(ggplot2)
g = ggplot(data=need_DEG,
aes(x=log2FoldChange, y=-log10(pvalue),
color=change)) +
geom_point(alpha=0.4, size=1.75) +
theme_set(theme_set(theme_bw(base_size=20)))+
xlab("log2 fold change") + ylab("-log10 p-value") +
ggtitle( this_tile ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
scale_colour_manual(values = c('blue','black','red')) ## corresponding to the levels(res\$change)
print(g)
}
``````