PCoA
library(vegan)
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-7
library(ape)
#Bray-curtis距离， cailliez校正
bray_dist<-vegdist(data,method = "bray")
df.pcoa<-pcoa(bray_dist,correction = "cailliez")
#获取用于画图的数据：df.pcoa$vectors
df.pcoa$vectors
## Axis.1 Axis.2 Axis.3 Axis.4 Axis.5
## MJ1 -0.468536986 0.254232456 0.163063529 0.113124480 -0.046735729
## MJ3 -0.412837254 0.209158470 -0.080265790 -0.201727392 0.273803473
## MJ7 -0.368316566 0.227760238 0.306012188 0.160730011 -0.143716890
## MJ2 -0.278161938 -0.116786189 0.001253774 0.035095275 0.095260931
## MJ4 -0.190678292 -0.034562408 -0.438212428 0.030592969 -0.231058789
## MJ5 -0.220250335 -0.086547474 -0.353659407 -0.050563161 0.017468519
## MJ14 0.032116417 -0.474817988 0.153381377 -0.028722486 0.015929577
## MJ15 -0.005839456 -0.396258090 0.216315602 -0.089883675 0.005778204
## MJ17 0.028827218 -0.407750135 0.040404508 -0.007491264 -0.035657471
## MJ19 0.348789590 0.037480413 -0.089637446 0.166955633 0.042807000
## MJ20 0.287351145 0.017315178 0.002231703 -0.042294698 -0.007040406
## MJ23 0.313945322 0.002292314 -0.045139473 0.197103786 0.056686863
## MJ27 0.304591746 0.180479194 0.037565111 -0.205487189 -0.087929570
## MJ29 0.285166963 0.341082655 0.109250561 -0.259934653 -0.099785323
## MJ30 0.343832425 0.246921365 -0.022563810 0.182502365 0.144189612
## Axis.6 Axis.7 Axis.8 Axis.9 Axis.10
## MJ1 -0.018726542 -0.05175946 -0.047040783 0.013675167 0.0639330095
## MJ3 0.041046975 0.09586258 0.010686725 -0.022888520 -0.0458229881
## MJ7 -0.005297990 0.02584892 0.004675678 0.027026570 -0.0600571378
## MJ2 -0.136643001 -0.07132687 0.073430768 -0.111888242 0.0535648028
## MJ4 0.133872058 0.05454079 0.075320075 -0.048617361 -0.0101863347
## MJ5 -0.077660220 -0.07623058 -0.082930462 0.137693613 0.0138302050
## MJ14 0.118193161 0.02277841 0.030553498 0.021202632 0.0813456294
## MJ15 0.064772629 0.04062246 0.040521479 0.070543958 -0.0271594807
## MJ17 -0.007593444 -0.03042825 -0.132549851 -0.089523551 -0.0796778939
## MJ19 -0.086362013 0.05395017 -0.060833019 -0.040843897 0.0096435916
## MJ20 -0.135517435 -0.18516598 0.110141597 0.026558811 -0.0585831617
## MJ23 -0.136080491 0.18720834 0.036786829 0.045538245 0.0021855925
## MJ27 -0.032437087 -0.01031281 -0.054484054 -0.029935802 0.0597315780
## MJ29 0.010821419 0.05337726 0.003781183 -0.001573728 0.0003093007
## MJ30 0.267611980 -0.10896495 -0.008059663 0.003032107 -0.0030567128
## Axis.11 Axis.12
## MJ1 -0.01576180 -0.0116042069
## MJ3 -0.04139736 -0.0019106363
## MJ7 -0.01048523 0.0080247336
## MJ2 0.06671996 0.0084256838
## MJ4 -0.00897824 0.0013671508
## MJ5 0.02500805 0.0020758220
## MJ14 -0.05863058 -0.0165575255
## MJ15 0.03669569 0.0277199230
## MJ17 0.01767772 -0.0163933390
## MJ19 -0.06585079 0.0249489036
## MJ20 -0.05041900 -0.0089516087
## MJ23 0.03626110 -0.0184358494
## MJ27 -0.01094794 0.0124359239
## MJ29 0.04315511 -0.0112497725
## MJ30 0.03695329 0.0001047978
df.plot<-data.frame(df.pcoa$vectors)
#获得相关信息（如坐标轴上的百分比为Rel_corr_eig）：df.pcoa$values
df.pcoa$values
## Eigenvalues Corr_eig Rel_corr_eig Broken_stick Cum_corr_eig Cum_br_stick
## 1 1.284861244 1.59914688 0.309550391 0.24462567 0.3095504 0.2446257
## 2 0.942903070 1.20127513 0.232533480 0.16770260 0.5420839 0.4123283
## 3 0.539666635 0.70983133 0.137403618 0.12924106 0.6794875 0.5415693
## 4 0.304620691 0.43430324 0.084069036 0.10360003 0.7635565 0.6451694
## 5 0.205708438 0.30630121 0.059291401 0.08436926 0.8228479 0.7295386
## 6 0.180037243 0.27614844 0.053454663 0.06898465 0.8763026 0.7985233
## 7 0.116577951 0.19279854 0.037320438 0.05616414 0.9136230 0.8546874
## 8 0.060669392 0.11893557 0.023022621 0.04517512 0.9366456 0.8998625
## 9 0.054079104 0.11205136 0.021690030 0.03555974 0.9583357 0.9354223
## 10 0.033768638 0.08732714 0.016904108 0.02701273 0.9752398 0.9624350
## 11 0.023867338 0.06989057 0.013528871 0.01932042 0.9887687 0.9817554
## 12 0.002914842 0.03801483 0.007358615 0.01232742 0.9961273 0.9940828
## 13 0.000000000 0.02000663 0.003872728 0.00591716 1.0000000 1.0000000
## 14 -0.009020685 0.00000000 0.000000000 0.00000000 1.0000000 1.0000000
## 15 -0.029475407 0.00000000 0.000000000 0.00000000 1.0000000 1.0000000
x_label<-round(df.pcoa$values$Rel_corr_eig[1]*100,2)
y_label<-round(df.pcoa$values$Rel_corr_eig[2]*100,2)
x_label
## [1] 30.96
y_label
## [1] 23.25
#添加分组信息
class = data.frame(class = c(rep('Seawater', 3),rep('7 days', 3),rep('14 days', 3),rep('21 days', 3),rep('1 month', 3) ))
data = cbind(df.plot, class)
Plot!
library(ggplot2)
#先调整一下图例顺序
data$class = factor(data$class, levels = c('Seawater', '7 days', '14 days', '21 days', '1 month'))
ggplot(data,aes(x=Axis.1,y=Axis.2, color = class, shape = class))+ geom_point(alpha=1, size=7) + theme_bw() + theme(panel.grid = element_blank()) + labs(x=paste0("PCoA1 (",x_label,"%)"), y=paste0("PCoA2 (",y_label,"%)")) + geom_vline(xintercept = 0, color = 'gray', size = 0.4) + geom_hline(yintercept = 0, color = 'gray', size = 0.4) + theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent', size = 1.2),legend.title=element_blank()) + theme(axis.text = element_text(colour = 'black', size = 20), axis.ticks = element_line(colour = 'black', size = 0.6), axis.title = element_text(size = 22), legend.key.size = unit(0.9, 'cm'), legend.text = element_text(size=18), axis.ticks.x = element_line(size=0.8), axis.ticks.y = element_line(size=0.8), axis.ticks.length = unit(0.18, 'cm') ) + scale_shape_manual(values = c(7,8,9,13,14)) + scale_colour_manual(values=c("coral4", "darkblue", "firebrick1", "darkorchid1","green4"))
Change the color and shape!
library(ggplot2)
#注意这里是fill而不是color
ggplot(data,aes(x=Axis.1,y=Axis.2, fill = class, shape = class))+ geom_point(alpha=0.8, size=7) + theme_bw() + theme(panel.grid = element_blank()) + labs(x=paste0("PCoA1 (",x_label,"%)"), y=paste0("PCoA2 (",y_label,"%)")) + geom_vline(xintercept = 0, color = 'gray', size = 0.4) + geom_hline(yintercept = 0, color = 'gray', size = 0.4) + theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent', size = 1.2),legend.title=element_blank()) + theme(axis.text = element_text(colour = 'black', size = 20), axis.ticks = element_line(colour = 'black', size = 0.6), axis.title = element_text(size = 22), legend.key.size = unit(0.9, 'cm'), legend.text = element_text(size=18), axis.ticks.x = element_line(size=0.8), axis.ticks.y = element_line(size=0.8), axis.ticks.length = unit(0.18, 'cm') ) + scale_shape_manual(values = c(21,22,23,24,25)) + scale_fill_manual(values=c("lightgoldenrod1", "cornflowerblue", "indianred1", "mediumturquoise","darkolivegreen3"))
Add labels!
data_labels = data.frame(rownames(data))
colnames(data_labels) = 'datalabels'
data = cbind(data, data_labels)
ggplot(data,aes(x=Axis.1,y=Axis.2, fill = class, shape = class))+ geom_point(alpha=0.8, size=7) + theme_bw() + theme(panel.grid = element_blank()) + labs(x=paste0("PCoA1 (",x_label,"%)"), y=paste0("PCoA2 (",y_label,"%)")) + geom_vline(xintercept = 0, color = 'gray', size = 0.4) + geom_hline(yintercept = 0, color = 'gray', size = 0.4) + theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent', size = 1.2),legend.title=element_blank()) + theme(axis.text = element_text(colour = 'black', size = 20), axis.ticks = element_line(colour = 'black', size = 0.6), axis.title = element_text(size = 22), legend.key.size = unit(0.9, 'cm'), legend.text = element_text(size=18), axis.ticks.x = element_line(size=0.8), axis.ticks.y = element_line(size=0.8), axis.ticks.length = unit(0.18, 'cm') ) + scale_shape_manual(values = c(21,22,23,24,25)) + scale_fill_manual(values=c("lightgoldenrod1", "cornflowerblue", "indianred1", "mediumturquoise","darkolivegreen3")) + geom_text(aes(label=datalabels))
#让标签不重叠！
library(ggrepel)
ggplot(data,aes(x=Axis.1,y=Axis.2, fill = class, shape = class))+ geom_point(alpha=0.8, size=7) + theme_bw() + theme(panel.grid = element_blank()) + labs(x=paste0("PCoA1 (",x_label,"%)"), y=paste0("PCoA2 (",y_label,"%)")) + geom_vline(xintercept = 0, color = 'gray', size = 0.4) + geom_hline(yintercept = 0, color = 'gray', size = 0.4) + theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent', size = 1.2),legend.title=element_blank()) + theme(axis.text = element_text(colour = 'black', size = 20), axis.ticks = element_line(colour = 'black', size = 0.6), axis.title = element_text(size = 22), legend.key.size = unit(0.9, 'cm'), legend.text = element_text(size=18), axis.ticks.x = element_line(size=0.8), axis.ticks.y = element_line(size=0.8), axis.ticks.length = unit(0.18, 'cm') ) + scale_shape_manual(values = c(21,22,23,24,25)) + scale_fill_manual(values=c("lightgoldenrod1", "cornflowerblue", "indianred1", "mediumturquoise","darkolivegreen3")) + geom_text_repel(aes(label=datalabels))
