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
网友评论