R使用笔记: scatterplot with confidence ellipses;envfit的实现及释意

笔记内容:

  • 用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统计检验结果,能更直观的看出有差异的组别。

ggplot2: stat_ellipse()
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这两个物种类别。

ggplot2 + envfit
用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可能会稍繁琐一些。要根据实际情况选择使用哪种。

ggpubr: ggscatter()

彩蛋!!!!!!!!!QAQ

ggsci包就是各种sci-fi主题的调色盘!!!!有人做了我一直想做的事!!!!!!
除了辛普森一家,星际迷航,也有各种大牌期刊里图表的配色,与ggplot2,ggpubr无缝兼容,直接使用scale_color_rickandmorty()scale_fill_rickandmorty()

"rickandmorty" palette

参考链接:
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

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 205,033评论 6 478
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 87,725评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 151,473评论 0 338
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,846评论 1 277
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,848评论 5 368
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,691评论 1 282
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,053评论 3 399
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,700评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 42,856评论 1 300
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,676评论 2 323
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,787评论 1 333
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,430评论 4 321
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,034评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,990评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,218评论 1 260
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,174评论 2 352
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,526评论 2 343

推荐阅读更多精彩内容