R语言与医学统计图形-【33】生存曲线、森林图、曼哈顿图

1.生存曲线

基础包survival+扩展包survminer。

survival包内置肺癌数据集lung。

library(survival)
library(survminer)
str(lung)

#拟合模型
fit <- survfit(Surv(time,status)~sex,data=lung)
#绘制生存曲线
ggsurvplot(fit,
           pval = TRUE, #添加log rank检验的p值
           conf.int = TRUE, #添加置信区间
           risk.table = TRUE, #添加风险表
           risk.table.col='strata', #根据数据分组为风险表添加颜色
           linetype = 'strata', #不同组别生存曲线线型
           surv.median.line = 'hv', #标注中位生存时间
           ggtheme = theme_bw(), #主题风格
           palette = c("#E7B800","#2E9FDF")) #颜色风格

image.png

生存曲线进一步修饰。

ggsurvplot(fit,
           pval = TRUE,
           conf.int = TRUE,
           conf.int.style='ribbon', #置信区间风格
           xlab='Time in days', 
           break.time.by=200, #将x轴按200为间隔切分
           ggtheme = theme_light(),
           risk.table = 'abs_pct', #风险表中添加绝对数和相对数
           risk.table.y.text.col=TRUE, #风险表文字颜色
           risk.table.y.text=FALSE, #以条图展示风险表标签(非文字)
           ncensor.plot=TRUE, #展示随访中不同时间点死亡和删失的情况
           surv.median.line = 'hv',
           legend.labs=c('Male','Female'), #改变图例标签
           palette = c("#E7B800","#2E9FDF"))

image.png
其他类型生存曲线绘制。
fun参数默认pct (survival probability in percentage),即累计生存率,此外fun还可设为cumhaz(cumulative hazard)累计风险,或event(cumulative events)累计的结局事件发生数,即累计死亡人数。

#累计风险
ggsurvplot(fit,
           conf.int = TRUE,
           risk.table.col="strata",
           ggtheme = theme_bw(),
           palette = c("#E7B800","#2E9FDF"),
           fun = 'cumhaz') 

#累计死亡人数
ggsurvplot(fit,
           conf.int = TRUE,
           risk.table.col="strata",
           ggtheme = theme_bw(),
           palette = c("#E7B800","#2E9FDF"),
           fun = 'event') 
#注意y轴是累计死亡人数占总人数比例

image.png

2. meta分析森林图

常见meta分析(重分析)图形,如森林图、漏斗图。

#install.packages('meta')
library(meta)
data("Olkin95") #内置

#metabin进行二分类数据的meta分析
metal <- metabin(event.e, #实验组事件发生数
                 n.e, #实验组样本量
                 event.c, #对照组事件发生数
                 n.c, #对照组样本量
                 data = Olkin95,subset = c(41,47,51,59), #选取4篇文章meta分析
                 sm="RR", #汇总指标
                 method = "I", #结果汇总方法,I表倒方差法
                 studlab = paste(author,year))#森林图上文献标识

forest(metal)

image.png

认识一下普通meta分析森林图的几大基本元素。图形的最左边,即第一部分,是纳入文献的标识,此处使用的是第一作者姓名加上发表日期,这个标识是由metabin()函数中的studlab参数来决定的。依次往右,第二部分是纳入文献中的原始数据,由于此处使用的是二分类数据,所以有实验组和对照组之分,在每个组中,又列出了事件发生数和总样本数,即events和total。图形第三部分是森林图的图形主干,展示的是每篇纳入文献的点估计和区间估计,大家注意,在图中,纳入文献的点估计是用方块表示,而在所有方块下面,就是汇总的结果,用菱形表示,垂直于菱形中央的虚线表示的是汇总结果的点估计值,菱形的延伸范围表示汇总结果的区间估计值。继续往右,第四部分展示了每篇纳入文献的RR值的点估计和区间估计值。图形第五部分展示的是不同模型下(固定效应模型和随机效应模型)的每篇纳入文献所占权重。第六部分位于图形的下方,可以看到图中展示了固定效应模型和随机效应模型下的汇总结果,并且展示了文献之间的异质性和异质性检验结果。

forest函数的参数非常多,下面对其中主要参数进行设置。

修饰的森林图。

forest(metal,
       studlab = TRUE, #添加文献标识
       layout = 'RevMan5', #布局模式,JAMA
       comb.fixed = TRUE, #展示固定效应模型结果
       comb.random = FALSE, #展示随机效应模型结果
       overall = TRUE, #展示汇总值
       lty.fixed = 2,#固定效应模型点估计值的线型,1-6
       prediction = TRUE,#展示汇总值预测区间
       xlab = 'RR estimates',
       xlab.pos = 1, #x轴标签中点的位置
       xlim = c(0.1,5), 
       ref = 1, #设置参考值
       lab.e = '实验组',
       lab.c='对照组',
       col.square='forestgreen', #每篇文献对应方形颜色
       col.square.lines='blue',#每篇文献对应区间线条颜色
       col.diamond.fixed='deeppink1',#固定效应模型汇总值的菱形颜色
       test.overall=TRUE, #汇总值统计检验
       fs.study=10, #每篇文献结果字体大小
       ff.study='italic',
       fs.study.label=11, #文献标识字体大小
       ff.study.label='bold',
       fs.axis=5,
       ff.axis='italic',
       ff.smlab='bold.italic',
       ff.fixed='plain',
       ff.hetstat='plain',
       digits=3
       )

image.png

统计结果汇总森林图

分组汇总结果。
使用forestplot包。

3.曼哈顿图

展示全基因组关联分析(GWAS)结果。

#install.packages('qqman')
library(qqman)
gwasdata <- data.frame(SNP=1:50000,
                       CHR=c(paste0('Chr',rep(1:10,each=50000))),
                       BP=rep(1:50000,10), #染色体中SNP位点的位置
                       P=abs(rnorm(50000*10,0.2,0.3))) #SNP位点的统计检验P值
library(RColorBrewer)

manhattan(gwasdata,
          col=brewer.pal(10,'Spectral'), #定义x轴染色体颜色
          cex=0.6,
          chrlabs = rep(paste0('Chr',1:10)), #自定义染色体标号
          suggestiveline = -log10(1e-5), #添加指示线,默认-log10(1e-5)
          genomewideline = -log10(1e-7))#添加P值显著线,默认-log10(5e-10)

#内置数据
str(gwasResults)
manhattan(gwasResults,
          highlight = snpsOfInterest, #需要标亮的SNP位点
          col=c("#A4D3EE","#DDA0DD"))

#单条染色体
manhattan(subset(gwasResults,CHR==2),col="#DDA0DD")

image.png

原文地址:https://www.cnblogs.com/jessepeng/p/12354311.html