笔记内容:
- 用ggplot2绘制scatterplot with confidence ellipses
- envfit的R实现及释意
- 用ggpubr绘制scatterplot with confidence ellipses
- ggsci的调色板彩蛋!!QAQ
在微生物β-diversity分析中常用距离矩阵(unifrac)做PcoA聚类分析,以观察不同组间物种构成的差异。结合针对距离矩阵的MANOVA统计检验结果,总结出组间物种构成差异的显著性。通过将各点PcoA结果坐标绘制在二维平面图上,并加上confidence ellipses进行可视化。
##使用内置数据集iris演示
library(ggplot2)
library(ade4)
data(iris)
iris <- data.frame(iris)
iris_sub <- iris[,1:4]
group <- as.factor(iris$Species)
mds <- cmdscale(dist(iris_sub), k = 2, eig = TRUE)
mds_point <- data.frame(mds$points) # 得到各样本的坐标
colnames(mds_point) <- c('X1','X2')
eig <- mds$eig
用ggplot2绘制scatterplot with confidence ellipses
color <- c(brewer.pal(3,"Set1"))
ggplot(mds_point, aes(x = X1, y = X2, color = group)) +
geom_point(aes(color = group), size = 4, alpha = 0.6) +
stat_ellipse(aes(x = X1, y = X2, fill = group), geom = "polygon", alpha = 1/2, levles = 0.95) +
# geom用于设置填充形状,alpha设置透明度。不设置则为实心填充,遮盖椭圆中的点, levels设置confidence ellipses的置信区间, 在0-1范围内。levels越小椭圆面积越小,涵盖的点越集中。
# 不需要填充的时候去掉fill及对fill的补充参数geom,alpha等即可
scale_fill_manual(values= color) +
scale_color_manual(values = color)
# 颜色可以自己设置,或者直接用scale_color_brewer()
置信椭圆的算法复杂,背后有很多繁杂的数学原理,在这里不深究,只是了解置信椭圆在本例中用于高维数据中,像95%置信区间一样,在一组数据中随机抽取一个样本,其落在这个区间(域)内的概率为95%。在图中展现了数据集中的范围。再结合MANOVA统计检验结果,能更直观的看出有差异的组别。
envfit的R实现及释意
library(vegan)
fit <- envfit(mds, iris_sub,permutations = 999)
fit_val <- scores(fit, display = c("vectors"))
fit_val <- fit_val*vegan::ordiArrowMul(fit_val, fill = 1.5)
## fit的结果output:
## P值表示显著性,Dim1, Dim2的坐标用于接下来把向量添加到PcoA二维图中。
> fit$vectors
Dim1 Dim2 r2 Pr(>r)
Sepal.Length 0.48219 0.87607 0.9579 0.001 ***
Sepal.Width -0.11499 0.99337 0.8400 0.001 ***
Petal.Length 0.98013 -0.19836 0.9981 0.001 ***
Petal.Width 0.97852 -0.20615 0.9366 0.001 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Permutation: free
Number of permutations: 999
##ggplot添加向量
ggplot(mds_point, aes(x = X1, y = X2, color = group)) +
geom_point(size = 4, alpha = .6) +
stat_ellipse(aes(x = X1, y = X2, fill = group), geom = "polygon", alpha = 1/2) +
scale_fill_manual(values=color) +
scale_color_manual(values = color) +
geom_segment(data=data.frame(fit_val),
aes(x=0,y=0,xend=Dim1, yend=Dim2),
arrow=arrow(length=unit(0.2,"cm")),
color='black',alpha=1) +
geom_label_repel(data=data.frame(fit_val), aes(Dim1, Dim2, label=rownames(fit_val)),
color='black',alpha=1,
segment.color = 'grey35',
point.padding = unit(0.1,"lines")) +
labs(x = paste("PCoA 1 (", format(100*eig[1]/sum(eig), digits = 4), "%)",sep = ""),
y = paste("PCoA 2 (", format(100*eig[2]/sum(eig), digits = 4), "%)",sep = ""))
每个向量代表一个数据类型为连续的环境变量,在这里即4中不同的width和length。
每个向量的方向环境变量发生变化最快的方向;向量的长度表示环境变量与分组的关系。其与PcoA绘制在一起的意义为结合关键环境变量解读,通过向量的方向,长短及显著性,了解哪些环境变量的变化会趋向构成组间差异,且这样的差异有意义。起到筛选关键环境变量的作用。在下图中,可以解读为Petal.Length和Petal.Width值越大,样本越趋近versicolor和virginica这两个物种类别。
用ggpubr绘制scatterplot with confidence ellipses
library(ggpubr)
mds_point <- cbind(mds_point, iris$Species)
colnames(mds_point)[3] <- "group"
ggscatter(mds_point, x= "X1", y = "X2",
color = "group", palette = "Set1", # 任何存在的palette都可以,不仅仅是brewer.pal中的,可以调用ggsci中的各种sci-fi主题颜色包!
ellipse = TRUE, # 设置是否需要confidence ellipses
mean.point = TRUE, star.plot = TRUE, # 设置confidence ellipses中心是否与所有点连线
ellipse.level = 0.95, # 设置confidence level可以调整椭圆的大小
ggtheme = theme_minimal()) +
labs(x = paste("PCoA 1 (", format(100*eig[1]/sum(eig), digits = 4), "%)",sep = ""),
y = paste("PCoA 2 (", format(100*eig[2]/sum(eig), digits = 4), "%)",sep = "")) +
ggtitle("PcoA")
ggpubr中的ggscatter(), ggboxplot()...等等函数,在绘制较复杂设置较多的图时比较方便,参数一目了然。如果用ggplot2可能会稍繁琐一些。要根据实际情况选择使用哪种。
彩蛋!!!!!!!!!QAQ
ggsci包就是各种sci-fi主题的调色盘!!!!有人做了我一直想做的事!!!!!!
除了辛普森一家,星际迷航,也有各种大牌期刊里图表的配色,与ggplot2,ggpubr无缝兼容,直接使用scale_color_rickandmorty()
及scale_fill_rickandmorty()
参考链接:
http://ggplot2.tidyverse.org/reference/stat_ellipse.html
https://cran.r-project.org/web/packages/ggsci/vignettes/ggsci.html
vegan包中对ordination的文档:
https://cran.r-project.org/web/packages/vegan/vignettes/intro-vegan.pdf