第二步：整洁数据做映射操作，确定x,y,color,size,shape,alpha等
第三步：选择合适的几何对象（根据画图的目的、变量的类型和个数）
第四步：坐标系和刻度配置
第五步：标签信息和图例信息
第六步：选择合适的主题
ggplot2的语法包括10个部件。
数据（data）
映射（mapping）
几何对象（geom）
标度（scale）
统计变换（stats）
坐标系（coord）
位置调整（Position adjustments）
分面（facet）
主题（theme）
输出（output）
前3个是必须的，其它部件ggplot2会自动配置，也可以手动配置
ggplot2基本绘图模板：
注意：
1）添加图层的加号（+）只能放在行末尾
2）红色方框里面mapping是全局域，绿色方框里面mapping是局部域，执行先后顺序，先局部域，后全局域
ggplot2画图必要部件-数据，映射和几何对象
2.1 数据
数据（Data）用于画图的整洁数据
library(tidyverse
ggplot()先只提供数据，创建一个空图形。
# ggplot()先提供整洁数据，生成一个空图形
2映射
映射，把数据变量集与图形属性库建立关联。
最常用的映射有：
x：x轴
y：y轴
color：颜色
size：大小
shape：形状
fill：填充
alpha：透明度
以mpg数据集为例，把变量displ和hwy分别映射到x和y，变量drv映射到color，此时图形就有了坐标轴和网格线，color需要在有了几何对象后才能体现出来。
# 映射操作
ggplot(data = mpg, mapping = aes(x = displ,
y = hwy, color = drv))
2.3 几何对象
几何对象是表达数据的视觉对象
不同类型的几何对象是从不同的角度表达数据。
pgglot2提供了50多种“几何对象”，均以geom_xxxx()的方式命名，常用的有：
几何对象很简单，只需要添加图层即可。
例如，以mpg数据集为例，画散点图。
ggplot(data = mpg, mapping = aes(x = displ,
y = hwy,
color = drv)) +
geom_point()层依次叠加，在上图的基础上，再添加一个几何对象：光滑曲线。
#继续增加一个几何对象：光滑曲线
# 写法1
ggplot(data = mpg, mapping = aes(x = displ,
y = hwy,
color = drv)) +
geom_point() +
geom_smooth(se=FALSE)
# 写法2
ggplot(data = mpg, mapping = aes(x = displ, y = hwy)) +
geom_point(aes(color = drv)) +
geom_smooth(se=FALSE)
思考题：
1）写法1和写法2的差异？（全局域和局部域的使用差异）
2）写法2若是要实现写法1的功能，怎么编写代码？
03
标度
ggplot2会自动根据输入变量选择最优的坐标刻度方法，若要手动设置或调整，就需要使用标度函数。
标度函数用来控制几何对象中的标度映射（x轴，y轴或者由color，fill，shape，size产生的图例）。
ggplot2提供丰富的标度函数，常用的有：
拓展功能：scales包提供很多设置刻度标签风格的函数，比如百分数、科学计数法法、美元格式等。
3.1 修改坐标轴刻度及标签
连续变量使用scale_*_continuous()函数，参数breaks设置各个刻度的位置，参数labels设置各个刻度对应的标签。
离散变量使用scale_*_discrete()函数，修改离散变量坐标轴的标签。
时间变量使用scale_x_date()函数设置日期刻度，参数date_breaks设置刻度间隔，date_labels设置标签的日期格式
以mpg数据集为例，修改连续变量坐标轴刻度及标签。
# scale_y_continuous函数
# 对比分析和观察
# 图1
ggplot(mpg, aes(displ, hwy)) +
geom_point()
# 图2
ggplot(mpg, aes(displ, hwy)) +
geom_point() +
scale_y_continuous(breaks = seq(15, 40, by = 10))
# 图3
ggplot(mpg, aes(displ, hwy)) +
geom_point() +
scale_y_continuous(breaks = seq(15, 40, by = 10),
labels = c(" 一五 "," 二五 "," 三五 "))
以mpg数据集为例，修改离散变量的标签
# scale_x_discrete函数
# 对比分析和观察
# 图1
ggplot(mpg, aes(x = drv)) +
geom_bar()
# 图2
ggplot(mpg, aes(x = drv)) +
geom_bar() +
scale_x_discrete(labels = c("4" = " 四驱 ", "f" = " 前驱 ",
"r" = " 后驱 "))
以ggplot2自带的economics数据集为例，修改日期变量。
# scale_x_date函数
# 以ggplot2自带的economics为例
economics %>% glimpse()
# 图1
ggplot(tail(economics, 45), aes(date, uempmed / 100)) +
geom_line()
# 图2
ggplot(tail(economics, 45), aes(date, uempmed / 100)) +
geom_line() +
scale_x_date(date_breaks = "6 months", date_labels = "%Y-%b") +
scale_y_continuous(labels = scales::percent)
3.2 修改坐标轴标签、图例名及图例位置
用labs()函数参数x，y或者xlab()，ylab()，设置x轴，y轴标签。
若用参数color生成了图例，可以在labs()函数用参数color修改图例名。
用theme图层的参数legend.position设置图例的位置。
以mpg数据为例。
# 修改坐标轴标签，图例名和图例位置
mpg
# 图1
ggplot(mpg, aes(displ, hwy)) +
geom_point(aes(color = drv)) +
labs(x = " 引擎大小 (L)", y = " 高速燃油率 (mpg)",
color = " 驱动类型 ") +
theme(legend.position = "top")
# 图2
ggplot(mpg, aes(displ, hwy)) +
geom_point(aes(color = drv)) +
xlab(" 引擎大小 (L)") +
ylab(" 高速燃油率 (mpg)") +
labs(color = " 驱动类型 ") +
theme(legend.position = "top")
# 图3 不需要图例
ggplot(mpg, aes(displ, hwy)) +
geom_point(aes(color = drv)) +
xlab(" 引擎大小 (L)") +
ylab(" 高速燃油率 (mpg)") +
theme(legend.position = "none")
3.3 设置坐标轴的范围
用coord_cartesian()函数参数xlim和ylim，或者用xlim()，ylim()设置x轴和y轴的范围。
以mpg数据集为例。
# 修改坐标轴的范围
# 图1 coord_cartesian()的参数xlim和ylim
ggplot(mpg, aes(displ, hwy)) +
geom_point(aes(color = drv)) +
coord_cartesian(xlim = c(5, 7), ylim = c(10, 30))
# 图2 xlim()和ylim()函数
ggplot(mpg, aes(displ, hwy)) +
geom_point(aes(color = drv)) +
xlim(5, 7) +
ylim(10, 30)
3.4 变换坐标轴
用scale_x_log10()函数变换坐标系，可以保持原始数据的坐标刻度。
# 修改坐标轴的范围
# 图1 coord_cartesian()的参数xlim和ylim
ggplot(mpg, aes(displ, hwy)) +
geom_point(aes(color = drv)) +
coord_cartesian(xlim = c(5, 7), ylim = c(10, 30))
# 图2 xlim()和ylim()函数
ggplot(mpg, aes(displ, hwy)) +
geom_point(aes(color = drv)) +
xlim(5, 7) +
ylim(10, 30)
3.5 设置图形标题
用labs()函数设置图形标题。
参数title 设置正标题
参数subtitle 设置副标题
参数caption 设置脚注标题（默认右下角）
# 设置标题
# mpg数据集为例
p <- ggplot(mpg, aes(displ, hwy)) +
geom_point(aes(color = drv)) +
geom_smooth(se = FALSE) +
labs(title = " 燃油效率与引擎大小的关系图 ",
subtitle = " 两座车 ( 跑车 ) 因重量小而符合预期 ",
caption = " 数据来自 fueleconomy.gov")
p
标题若要居中，采用theme图层设置。
p + theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
3.6 设置color、fill颜色
数据的某个维度信息可以通过颜色来表示。
可以直接使用颜色值，建议使用RColorBrewer（调色板）或者colorspace包。
1）连续变量
- 用scale_color_gradient()设置二色渐变色。
# 连续变量
# 图1 scale_color_gradient()函数
ggplot(mpg, aes(displ, hwy, color = hwy)) +
geom_point() +
scale_color_gradient(low = "green", high = "red")
- 用scale_color_distiller()设置调色板中的颜色
# 图2 scale_color_distiller()函数
ggplot(mpg, aes(displ, hwy, color = hwy)) +
geom_point() +
scale_color_distiller(palette = "Set1")
2）离散变量
- 用scale_color_manual()手动设置颜色，还可以修改图例及其标签信息
# 离散变量
# 图1 scale_color_manual()函数
ggplot(mpg, aes(displ, hwy, color = drv)) +
geom_point() +
scale_color_manual(" 驱动方式 ",
values = c("red", "blue", "green"),
breaks = c("4", "f", "r"))
ggplot(mpg, aes(displ, hwy, color = drv)) +
geom_point() +
scale_color_manual(" 驱动方式 ",
values = c("red", "blue", "green"),
labels = c(" 四驱 ", " 前驱 ", " 后驱 "))
-用scale_fill_brewer()调用调色板中的颜色
# 图2 scale_fill_brewer()函数
ggplot(mpg, aes(x = class, fill = class)) +
geom_bar() +
scale_fill_brewer(palette = "Dark2")
.7 添加文字标注
ggrepel包提供了geom_label_repel()函数或者geom_text_repel()函数，为图形添加文字标注。
操作步骤：
第一步：先准备好标记点的数据
第二步：增加文字标注图层，包括标记点的数据和标注的文字给label参数
# 设置文字标注信息
library(ggrepel)
# 选取每种车型 hwy 值最大的样本
best_in_class <- mpg %>%
group_by(class) %>%
slice_max(hwy, n = 1)
best_in_class %>% select(class, model, hwy)
ggplot(mpg, aes(displ, hwy)) +
geom_point(aes(color = class)) +
geom_label_repel(data = best_in_class,
aes(label = model))
04
计变换、坐标系和位置调整
.1 统计变换
统计变换是构建新的统计量而画图。
例如，条形图或直方图，是对数据分组的频数做画图；平滑曲线是对数据拟合模型的预测值画图。
gplot2可以把统计变换直接融入画图中，不必先在对数据做统计变换后再画图。
gplot2提供30多种统计，均以stats_xxx()的方式命名。
1）可在几何对象中直接使用的统计变换，直接使用几何对象就可以了。
能在几何对象创建的，而需要单独使用。
mpg数据集为例。
stat_summary()做统计绘图并汇总。
# 图1 stat_summary()做统计绘图并汇总
p <- ggplot(mpg, aes(x = class, y = hwy)) +
geom_violin(trim = FALSE, alpha = 0.5, color = "green")
p
p + stat_summary(fun = mean,
fun.min = function (x) {mean(x) - sd(x)},
fun.max = function (x) {mean(x) + sd(x)},
geom = "pointrange",
color = "red")
tat_smooth()添加光滑曲线，与geom_smooth()相同。
参数method设置平滑曲线的拟合方法，如lm线性回归、glm广义线性回归、loess多项式回归、gam广义加法模型(mgcv包)、rlm稳健回归(MASS包)等。
参数formula指定平滑曲线方程，如y ~ x, y ~ poly(x, 2)， y ~ log(x)等。
参数se设置是否绘制置信区间。
# 图2 stat_smooth()添加平滑曲线
ggplot(mpg, aes(displ, hwy)) +
geom_point() +
stat_smooth(method = "lm",
formula = y ~ splines::bs(x, 3),
se = FALSE)
ggplot(mpg, aes(displ, hwy)) +
geom_point() +
geom_smooth(method = "lm",
formula = y ~ splines::bs(x, 3),
se = FALSE)
4.2 坐标系
ggplot2默认是直角坐标系。
- coord_cartesian()
常用的其它坐标系：
以mpg数据集为例，坐标轴翻转。
# 图1 坐标轴翻转coord_flip()
p <- ggplot(mpg, aes(class, hwy)) +
geom_boxplot()
p
p + coord_flip()
直角坐标下条形图转换为极坐标下玫瑰图。
# 图2 直角坐标条形图-->极坐标玫瑰图
p <- ggplot(mpg, aes(class, fill = drv)) +
geom_bar()
p
p + coord_polar()
4.3 位置调整
条形图的位置调整
# 图1：条形图条形位置调整
ggplot(mpg, aes(class, fill = drv)) +
geom_bar()
ggplot(mpg, aes(class, fill = drv)) +
geom_bar(position = "dodge")
ggplot(mpg, aes(class, fill = drv)) +
geom_bar(position = position_dodge(preserve = "single"))
散点图的散点位置调整
# 图1：散点图的散点位置调整
ggplot(mpg, aes(displ, hwy)) +
geom_point()
ggplot(mpg, aes(displ, hwy)) +
geom_point(position = "jitter")
用patchwork包排布多个图形
library(patchwork)
p1 <- ggplot(mpg, aes(displ, hwy)) +
geom_point()
p2 <- ggplot(mpg, aes(drv, displ)) +
geom_boxplot()
p3 <- ggplot(mpg, aes(drv)) +
geom_bar()
p1 | (p2 / p3)
p1 | p2 | p3
p1 / p2 / p3
p1 / (p2 | p3)
05
分面
利用分类变量把图形分成若干“子图”（面），实际上就是对数据分组后再画图，属于数据分析里面细分和下钻的思想。
5.1 用facet_wrap()函数
封装分面，先生成一维的面板系列，再封装到二维中。
语法形式：~ 分类变量 或者 ~ 分类变量1 + 分类变量2
参数scales设置是否共用坐标刻度，fixed 默认 共用， free 不共用，还可以额通过free_x，free_y单独设置。
cor()函数可以提供双变量之间的相关系数，还可以用scatterplotMatrix()函数生成散点图矩阵 不过R语言没有直接给出偏相关的函数； 我们要是做的话，要先调用cor.test()对变量进行Pearson相关性分析， 得到简单相关系数，然后做t检验，判断显著性。孟德尔随机化（Mendelian randomization,MR）是以孟德尔独立分配定律为基础进行流行病学研究设计和数据分析，论证病因假说的一种方法。由基因型决定中间表型（暴露）的差异， 因果方向明确。通过引入一个称之为工具变量的中间变量，来分析暴露因素和结局之间的因果关系
2.孟德尔随机化 vs RCT
孟德尔随机化的目的不是估计遗传效应的大小，而是估计暴露对结果的因果效应，所以与遗传变异相关的结局的平均变化幅度可能与干预措施导致的变化幅度不同
即使遗传变异与结果之间的关联程度很小，暴露的人群归因风险也不一定很低，因为暴露可能会比遗传变异解释更大的变化程度（例如，他汀类药物对低密度脂蛋白胆固醇水平的影响比低密度脂蛋白胆固醇水平与HMGCR基因变异的关联要大几倍，因此对后续结果的影响更大。）
孟德尔随机化要求大样本研究，变异发生率不能太小（最小等位基因频率MAF>5%)
3.工具变量
工具变量本身是一个计量经济学的概念，在孟德尔随机中，遗传变异被用作工具变量评估暴露对结局的因果效应，遗传变异满足工具变量的基本条件总结为(孟德尔随机化核心假设)：
关联性假设——遗传变异与暴露有关
独立性假设——该遗传变异与暴露-结果关联的任何混杂因素均不相关
排他性假设——该遗传变异不会影响结果，除非可能通过与暴露的关联来实现
某研究组想了解非洲村落里的儿童补充维生素A和其死亡情况的关联，如果仅仅利用维生素A的服用情况和死亡情况去判断两者的关联，那极有可能会产生很大的偏倚，这是因为维生素A的服用情况和很多潜在因素相关，比如家庭的经济困难程度、家庭成员以及实验儿童的依从性，而这些潜在的因素也可能对儿童的身体健康有很大的影响。因此，在研究起始设计中，研究者便利用工具变量来解决这个问题。
在这里，工具变量Z是指服用维生素A这个任务，类似于随机抽签。这样的话工具变量Z便只和X服用维生素A这个行为相关，与除X以外的混杂因素不相关。
4.应用范围
行为因素与健康：基因变异引起各个倾向某行为，决定暴露状态。如ALDH2变异引起乙醛代谢障碍，改变饮酒行为，不同ALDH基因型代表饮酒量多少；
机体代谢产物与疾病关系，估计长期效应。代谢产物是基因表达的中间表型，酶的底物或者体外难测量的代谢指标：如LDL受体基因变异引起家族高胆固醇血症，比较不同基因型之间CHD发病情况的差异，可模拟血胆固醇水平和CHD发病关系；
子宫内环境暴露于子代健康关系。
5.发文分析
孟德尔随机化研究均发表在影响因子5分以上的期刊中
6.基础分析流程——TwoSampleMR
找工具变量，我们要的是基因作为工具变量，这些基因都是从别人的研究中挑出来的，所有的基因研究有个专门的库叫做genome wide association studies (GWAS)。我们需要做的就是从这个库中挑出来我们自己需要的和我们暴露相关的基因变量SNPs。
估计工具变量对结局的作用，工具变量对结局的作用也是从所有的研究中估计出来的整体效应，这样可以拒绝单个研究的偏倚。
合并多个SNP的效应量，这个效应量是我们得到暴露和结局因果效应的前提。
处理数据，用合并后的数据进行孟德尔随机化分析和相应的敏感性分析。
7.TwoSampleMR代码实现
安装相关R包
install.packages('devtools')
library('devtools')
install_github("MRCIEU/TwoSampleMR") #安装TwoSampleMR包
library('TwoSampleMR')
devtools::install_github("mrcieu/ieugwasr",force = TRUE)
获取MR base的表型ID，将结果保存为pheno_info.csv这个文件
ao <-available_outcomes(access_token=NULL) #获取GWAS数据，但近期Google限制，容易被墙
write.csv(ao,'pheno_info.csv',row.names=F）#将数据写入本地存储
查看pheno_info.csv文件，获取与暴露相关的工具变量的信息以及结局信息。这里选择暴露为obesity class 2 （ID = 91）， 结局为 type 2 diabetes (ID = 1090)
exp_dat <- extract_instruments(outcomes=91,access_token=NULL)
obesity_exp_dat <- clump_data(exp_dat)
t2d_out_dat <- extract_outcome_data(snps=obesity_exp_dat$SNP, outcomes=1090, access_token=NULL)#提取结果信息
dat <- harmonise_data(exposure_dat =obesity_exp_dat, outcome_dat= t2d_out_dat)#数据合并，计算基因对结局的合并效应量
孟德尔随机化
results <- mr(dat)
OR值
OR <- generate_odds_ratios(results)
异质性检验
heterogeneity<- mr_heterogeneity(dat)
多效性检验
pleiotropy<- mr_pleiotropy_test(dat)
逐个剔除检验
leaveoneout<- mr_leaveoneout(dat)
散点图
mr_scatter_plot(results,dat)
森林图
results_single<- mr_singlesnp(dat)
mr_forest_plot(results_single)
漏斗图
mr_funnel_plot(results_single)
实例解析
2022年10月10日
西安交通大学生物医学信息与基因组学中心杨铁林教授团队在Nature Neuroscience （IF=28.771）期刊发表了题为：Mendelian randomization analyses support causal relationships between brain imaging-derived phenotypes and risk of psychiatric disorders 的文章。
研究背景
精神类疾病是一组脑功能紊乱的复杂疾病，会导致情感、认知和行为受到干扰和破坏。全球约有数亿人患有不同的精神障碍，被列为严重的公共卫生问题。近年来，脑影像学数据在脑疾病和功能的研究中受到广泛关注。以核磁共振成像为代表的脑影像技术，可用于活体无创定量评估人脑结构、连接和功能的特性。
虽然已有大量的观察性研究证据表明，精神疾病患者与健康正常人的脑影像表型存在显著差异，但脑影像学数据与精神障碍发病机制的因果关系尚不明确，探讨脑影像表型对精神疾病的因果作用具有重要的生物学和临床研究意义。
研究方法和结果
该研究基于大规模基因组数据，对常见的10种精神类疾病（包括注意力缺陷多动症、神经性厌食症、焦虑症、孤独症、双相情感障碍、抑郁症、强迫症、创伤后应激障碍、精神分裂症、抽动症）和587个关键的脑磁共振成像（MRI）结构表型进行了因果关系评估。
正向孟德尔随机化结果发现，脑白质纤维束的上额枕束的FA值和上放射冠的ICVF值、胼胝体内矢状层的MD值、第三脑室的体积等9个脑影像表型是精神分裂症、神经性厌食症和双相情感障碍的风险因素。进一步通过反向孟德尔随机化分析显示，发现精神分裂症的发生会导致额下回眶部的表面积和体积的增加。
该研究将基因组信息作为纽带，使脑影像表型和精神疾病联系起来，避免了观察性研究中由于药物或环境、生活方式等改变引起的样本检测数据偏差的缺点，确保了研究结果的稳健性。
![[object HTMLTableSectionElement]是什么](/aiimages/%5Bobject+HTMLTableSectionElement%5D%E6%98%AF%E4%BB%80%E4%B9%88.png)
