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")) #顏色風格
生存曲線進一步修飾。
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"))
其他類型生存曲線繪制。
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軸是累計死亡人數占總人數比例
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)
認識一下普通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
)
統計結果匯總森林圖
分組匯總結果。
使用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")
