排序出图三个重复带来的一系列问题的解决方案
wentao
2021/8/8
写在前面
目前patchwork包拼图很不错，我尝试将群落两两比较结果和排序图放到一起看看；
三个样本的排序图置信椭圆绘制这个问题，我来做两个方案；
三个重复的置换检验对于校正p值不要要求太高，低于0.2就可以了。
这是目前我唯一不使用phloseq封装的排序教程；
本文会交给你，如何在图层覆盖的情况下显示覆盖的内容；
library(vegan)
library(tidyverse)
pairwise.adonis1 <-function(x,factors,p.adjust.m)
{
x = x %>% as.matrix()
co = as.matrix(combn(unique(factors),2))
pairs = c()
F.Model =c()
R2 = c()
p.value = c()
for(elem in 1:ncol(co)){
ad = adonis(x[factors %in%c(as.character(co[1,elem]),as.character(co[2,elem])),factors %in%c(as.character(co[1,elem]),as.character(co[2,elem]))] ~
factors[factors %in%c(as.character(co[1,elem]),as.character(co[2,elem]))] ,permutations = 999);
pairs =c(pairs,paste(co[1,elem],'vs',co[2,elem]));
F.Model =c(F.Model,ad$aov.tab[1,4]);
R2 = c(R2,ad$aov.tab[1,5]);
p.value = c(p.value,ad$aov.tab[1,6])
}
p.adjusted =p.adjust(p.value,method=p.adjust.m)
pairw.res = data.frame(pairs,F.Model,R2,p.value,p.adjusted)
return(pairw.res)
}
#---导入数据
otu =read.csv("./otu.csv",row.names = 1) %>% t() %>% as.data.frame()
head(otu)## sample1 sample2 sample3 sample4 sample5 sample6
## OTU1 0.45724487 0.47512130 0.48076031 0.42379195 0.44256752 0.39355241
## OTU2 0.20552670 0.24318511 0.18764187 0.18471523 0.19670737 0.18976133
## OTU3 0.22758158 0.10338618 0.19905067 0.02797361 0.07637440 0.09065376
## OTU4 0.04945469 0.33300000 0.04415425 0.01939449 0.04472809 0.04478728
## OTU5 0.06160738 0.05617967 0.06110291 0.07147906 0.05521273 0.04514701
## OTU6 0.01783519 0.03173110 0.01958146 0.04212304 0.03914910 0.05081287
## sample7 sample8 sample9 sample10 sample11 sample12
## OTU1 0.38670554 0.35410840 0.41203157 0.29126527 0.33256955 0.21680160
## OTU2 0.23526525 0.20533941 0.20383444 0.20279485 0.23199267 0.25137003
## OTU3 0.14723778 0.20033434 0.05784600 0.22615791 0.07583567 0.13916428
## OTU4 0.05568320 0.07731510 0.03570940 0.08116948 0.05522967 0.05496949
## OTU5 0.02754142 0.03477204 0.04467715 0.04105715 0.03546474 0.03578279
## OTU6 0.04480170 0.03549645 0.04322292 0.05857946 0.05109446 0.07431810map = read.csv(
"./map.csv"
)
head(map)## 锘縎ample.ID Group
## 1 sample1 A
## 2 sample2 A
## 3 sample3 A
## 4 sample4 B
## 5 sample5 B
## 6 sample6 Bcolnames(map)[1] = "ID"
row.names(map) = map$ID
idx = rownames(map) %in% colnames(otu)
map1 = map[idx,]
map1## ID Group
## sample1 sample1 A
## sample2 sample2 A
## sample3 sample3 A
## sample4 sample4 B
## sample5 sample5 B
## sample6 sample6 B
## sample7 sample7 C
## sample8 sample8 C
## sample9 sample9 C
## sample10 sample10 D
## sample11 sample11 D
## sample12 sample12 Dotu = otu[,rownames(map1)]
head(otu)## sample1 sample2 sample3 sample4 sample5 sample6
## OTU1 0.45724487 0.47512130 0.48076031 0.42379195 0.44256752 0.39355241
## OTU2 0.20552670 0.24318511 0.18764187 0.18471523 0.19670737 0.18976133
## OTU3 0.22758158 0.10338618 0.19905067 0.02797361 0.07637440 0.09065376
## OTU4 0.04945469 0.33300000 0.04415425 0.01939449 0.04472809 0.04478728
## OTU5 0.06160738 0.05617967 0.06110291 0.07147906 0.05521273 0.04514701
## OTU6 0.01783519 0.03173110 0.01958146 0.04212304 0.03914910 0.05081287
## sample7 sample8 sample9 sample10 sample11 sample12
## OTU1 0.38670554 0.35410840 0.41203157 0.29126527 0.33256955 0.21680160
## OTU2 0.23526525 0.20533941 0.20383444 0.20279485 0.23199267 0.25137003
## OTU3 0.14723778 0.20033434 0.05784600 0.22615791 0.07583567 0.13916428
## OTU4 0.05568320 0.07731510 0.03570940 0.08116948 0.05522967 0.05496949
## OTU5 0.02754142 0.03477204 0.04467715 0.04105715 0.03546474 0.03578279
## OTU6 0.04480170 0.03549645 0.04322292 0.05857946 0.05109446 0.07431810bray_curtis =vegan:: vegdist(t(otu),method = "bray", na.rm=TRUE)
pcoa <- otu %>%
t() %>%
vegan:: vegdist(method = "bray", na.rm=TRUE) %>%
cmdscale(k=2, eig=T)
points <- pcoa$points %>% as.data.frame() %>%
# as.tibble() %>%
dplyr::rename(.,x = "V1",y = "V2" )
head(points)## x y
## sample1 -0.11731023 -0.04846822
## sample2 -0.08290919 -0.06069551
## sample3 -0.10191448 -0.06509878
## sample4 0.06356767 -0.05047065
## sample5 0.01970246 -0.04969312
## sample6 0.19534848 -0.07646325eig = pcoa$eig
eig## [1] 8.972991e-02 6.050659e-02 3.402134e-02 2.660735e-02 9.079289e-03
## [6] 2.808910e-03 1.395223e-03 1.505786e-04 1.387779e-17 -2.373096e-04
## [11] -2.402774e-03 -4.738457e-03points = cbind(points, map1[match(rownames(points), map1$ID), ])
ado = adonis(bray_curtis~ map1$Group,permutations = 999,method="bray")
a = round(as.data.frame(ado$aov.tab[5])[1,1],3)
R2 <- paste("adonis:R ",a, sep = "")
b = as.data.frame(ado$aov.tab[6])[1,1]
p_v = paste("p: ",b, sep = "")
title = paste(R2," ",p_v, sep = "")
title## [1] "adonis:R 0.525 p: 0.001"# anosim.result<-anosim(bray_curtis, map1$Group,permutations = 999)
# summary(anosim.result)
pair_bray_adonis = pairwise.adonis1(bray_curtis,map1$Group, p.adjust.m= "bonferroni")
head(pair_bray_adonis)## pairs F.Model R2 p.value p.adjusted
## 1 A vs B 3.481549 0.4653514 0.1 0.6
## 2 A vs C 2.499122 0.3845322 0.1 0.6
## 3 A vs D 4.810213 0.5459815 0.1 0.6
## 4 B vs C 1.792541 0.3094567 0.2 1.0
## 5 B vs D 3.058955 0.4333439 0.1 0.6
## 6 C vs D 1.479269 0.2699756 0.3 1.0library("ggalt")
p1 <- ggplot(points, aes(x=x, y=y, fill = Group)) +
geom_point(alpha=.7, size=5, pch = 21) +
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=""),
title=title) +
geom_encircle(aes(x=x, y=y,group = Group,color = Group), linetype = 2,alpha = 0.1) +
theme_classic()
p1
p0 <- ggplot(points, aes(x=x, y=y, fill = Group)) +
geom_point(alpha=.7, size=5, pch = 21) +
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=""),
title=title) +
geom_encircle(aes(x=x, y=y,group = Group,color = Group), linetype = 2,alpha = 0.1) +
theme_classic()
p0
library(car)
library(ggforce)
n = 0.8
p1 <- ggplot(points, aes(x=x, y=y, fill = Group)) +
geom_point(alpha=.7, size=5, pch = 21) +
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=""),
title=title) +
geom_mark_ellipse(aes(x=x * n, y=y * n,fill=Group,label=Group),alpha=0.1) +
theme_classic()
p1
# ggsave("./cs.pdf",p2,width = 12,height =6)
# ggsave("./cs1.pdf",p3,width = 12,height =6)library(ggpubr)
library(patchwork)
tab2 <- ggtexttable(pair_bray_adonis, rows = NULL)
p2 = p1 | tab2
p2
p3 = p0 | tab2
p3
library(concaveman)
# install.packages("concaveman")
n = 0.8
p1 <- ggplot(points, aes(x=x, y=y, fill = Group)) +
geom_point(alpha=.7, size=5, pch = 21) +
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=""),
title=title) +
geom_mark_ellipse(aes(x=x * n, y=y * n,fill=Group,label=Group),alpha=0.1) +
# coord_cartesian(clip = "off") +
# stat_ellipse(level=0.95,show.legend=F) +
# geom_encircle(aes(x=x, y=y,group = Group,color = Group), linetype = 2,alpha = 0.1) +
theme_classic()
p1
library(concaveman)
# install.packages("concaveman")
n = 1
p1 <- ggplot(points, aes(x=x, y=y, fill = Group)) +
geom_point(alpha=.7, size=5, pch = 21) +
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=""),
title=title) +
geom_mark_ellipse(aes(x=x * n, y=y * n,fill=Group,label=Group),alpha=0.1) +
# coord_cartesian(clip = "off") +
# stat_ellipse(level=0.95,show.legend=F) +
# geom_encircle(aes(x=x, y=y,group = Group,color = Group), linetype = 2,alpha = 0.1) +
theme_classic()
p1
根际互作生物学研究室 简介
根际互作生物学研究室是沈其荣教授土壤微生物与有机肥团队下的一个关注于根际互作的研究小组。本小组由袁军副教授带领，主要关注：1.植物和微生物互作在抗病过程中的作用；2 环境微生物大数据整合研究；3 环境代谢组及其与微生物过程研究体系开发和应用。团队在过去三年中在 isme J， Microbiome， PCE，SBB，Horticulture Research等期刊上发表了多篇文章。欢迎关注微生信生物公众号对本研究小组进行了解。