-
Notifications
You must be signed in to change notification settings - Fork 10
/
Fig.S4.EPICON.PCA.compositionaldata.2019.03.19.Rmd
79 lines (58 loc) · 2.19 KB
/
Fig.S4.EPICON.PCA.compositionaldata.2019.03.19.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
---
title: "Untitled"
author: "Cheng Gao"
date: "3/19/2019"
output: html_document
---
```{r, message=FALSE, warning=FALSE}
###PCA of CLR##
wd<- "/Users/chengg/Google Drive/EPICON/Mycobiome/Fungal ITS/statistic/Total.fungi"
setwd(wd)
rm(list = ls())
load("EPICON.data.preparation.RC.bNTI.ted.2019.04.19.Rdata")
library(compositions)
library(zCompositions)
source("codaSeq_functions.R")
library(vegan)
library(ggplot2)
library(colorRamps)
library(plyr)
env<-env0
gr<-fung0$Fungi=="Fungi"
otu_table0<-data.frame(fung0[gr ,c(1:1251) ]) ### raw data of all samples including pseudo samples #
ID<-data.frame(fung0[gr ,c(-1:-1251) ]) ### raw data of all samples including pseudo samples #
row.names(otu_table0)<-row.names(ID)
otu_table1<-data.frame(t(otu_table0[ID$Kingdom=="Fungi",]))
env0$Treatment1<-factor(env0$Treatment1,levels=c("Control", "Pre_flowering_drought", "Post_flowering_drought"))
subs<-env$Habitat!="Prop_root"
otu<-otu_table1[subs, ]
env1<-env0[subs, ]
d.n0 <- cmultRepl(otu, label=0, method="CZM", output="counts")
d.n0.clr <- codaSeq.clr(d.n0, samples.by.row=TRUE)
pcx <- prcomp(d.n0.clr)
mvar.clr <- mvar(d.n0.clr)
xxxx<-summary(pcx)
yy<-xxxx$importance
#write.csv(yy,"yy.csv")
sum(xxxx$importance)
pc1=round(sum(pcx$sdev[1]^2)/mvar.clr, 3)*100
pc2=round(sum(pcx$sdev[2]^2)/mvar.clr, 3)*100
envpc<-data.frame(pcx$x, env1)
```
```{r, fig.height = 7, fig.width = 10}
df <- envpc[,c("PC1", "PC2", "Habitat")]
find_hull <- function(envpc) envpc[chull(envpc$PC1, envpc$PC2), ]
hulls <- ddply(df, "Habitat", find_hull)
ggplot(data=envpc, aes(x= PC1 * -1 , y= PC2 * -1,group=Habitat))+
geom_point(aes(shape=Habitat, colour=factor(TP)),size=2) +
geom_polygon(data=hulls, alpha=.1)+
labs(x = sprintf("PC1 (%.1f%%)", pc1), y = sprintf("PC2 (%.1f%%)", pc2))+
scale_shape_manual(values=c(0, 6, 3, 1))+
scale_colour_manual(values=blue2red(18))+theme_bw()+
guides(color=guide_legend(title= "Week"),
shape=guide_legend(title= "Compartment"))+
theme(legend.title = element_text(colour="black", size=15, face="bold"),
legend.text = element_text(colour="black", size=15, face="bold"),
axis.text=element_text(size=10,face="bold"),
axis.title=element_text(size=15,face="bold"))
```