{"id":291,"date":"2022-01-24T11:35:08","date_gmt":"2022-01-24T03:35:08","guid":{"rendered":"https:\/\/linguopeng.top\/?p=291"},"modified":"2022-01-24T11:39:10","modified_gmt":"2022-01-24T03:39:10","slug":"%e5%be%ae%e7%94%9f%e7%89%a9%e8%8f%8c%e7%be%a4cda%e4%b8%8erda%e5%88%86%e6%9e%90","status":"publish","type":"post","link":"https:\/\/linguopeng.top\/?p=291","title":{"rendered":"\u5fae\u751f\u7269\u83cc\u7fa4CDA\u4e0eRDA\u5206\u6790"},"content":{"rendered":"\n
# \u9996\u5148\u8981\u5b89\u88c5devtools\u5305\uff0c\u4ec5\u9700\u5b89\u88c5\u4e00\u6b21\ninstall.packages(\"devtools\")\n# \u52a0\u8f7ddevtools\u5305\nlibrary(devtools)\n# \u4e0b\u8f7dggvegan\u5305\ndevtools::install_github(\"gavinsimpson\/ggvegan\")\nlibrary(ggvegan)\notu.tab <- read.csv(\"otutab.txt\", row.names = 1, header=T, sep=\"\\t\")\nenv.data <- read.csv(\"new_meta.txt\", row.names = 1, fill = T, header=T, sep=\"\\t\")\n#transform data\notu <- t(otu.tab)\n#data normolization (Legendre and Gallagher,2001)\n##by log \u73af\u5883\u56e0\u5b50\u8fdb\u884clog\u8f6c\u5316\uff0c\u4ee5\u51cf\u5c11\u540c\u4e00\u79cd\u73af\u5883\u56e0\u5b50\u4e4b\u95f4\u672c\u8eab\u6570\u503c\u5927\u5c0f\u9020\u6210\u7684\u5f71\u54cd\u3002\nenv.data.log <- log1p(env.data)##\n##delete NA \u5220\u6389\u7a7a\u503c\u884c\nenv <- na.omit(env.data.log)\n###hellinger transform OTU\u8fdb\u884chellinger\u8f6c\u5316\uff0c\u4f7f\u6570\u636e\u5747\u4e00\u6027\u66f4\u597d\u3002\notu.hell <- decostand(otu, \"hellinger\")\n#DCA analysis   DCA\u5206\u6790\u770b\u5206\u6790\u7ed3\u679c\u4e2dLengths of gradient\u7684\u7b2c\u4e00\u8f74\u4e2d\u7684Axis lengths\u5927\u5c0f\uff0c\n#\u5982\u679c\u5927\u4e8e4.0\uff0c\u5c31\u5e94\u8be5\u9009CCA\uff0c\u5982\u679c\u57283.0-4.0\u4e4b\u95f4\uff0c\u9009RDA\u548cCCA\u5747\u53ef\uff0c\u5982\u679c\u5c0f\u4e8e3.0\uff0cRDA\u7684\u7ed3\u679c\u8981\u597d\u4e8eCCA\u3002\n#\u672c\u4f8b\u4e2d\uff0c\u6211\u4eec\u6570\u636e\u7684Axis lengths \u5927\u5c0f\u4e3a0.63859\uff0c\u6240\u6709\u6211\u4eec\u5e94\u8be5\u9009\u62e9RDA\u8fdb\u884c\u5206\u6790\nsel <- decorana(otu.hell)\nsel\notu.tab.0 <- rda(otu.hell ~ 1, env) #no variables\n#Axis \u7b2c\u4e00\u9879\u5927\u4e8e\u56db\u5e94\u8be5\u7528CCA\u5206\u6790\notu.tab.1<- rda(otu.hell ~ ., env)\n#\u6211\u4eec\u5728\u7b5b\u9009\u5b8cRDA\u548cCCA\u5206\u6790\u540e\uff0c\u6211\u4eec\u9700\u8981\u5bf9\u6240\u6709\u73af\u5883\u56e0\u5b50\u8fdb\u884c\u5171\u7ebf\u6027\u5206\u6790\uff0c\u5229\u7528\u65b9\u5dee\u81a8\u80c0\u56e0\u5b50\u5206\u6790\nvif.cca(otu.tab.1)\n#\u5220\u9664\u6389\u5171\u7ebf\u6027\u7684\u73af\u5883\u56e0\u5b50\uff0c\u5220\u6389\u6700\u5927\u7684\u53d8\u91cf\uff0c\u76f4\u5230\u6240\u6709\u7684\u53d8\u91cf\u90fd\u5c0f\u4e8e10\notu.tab.1 <- rda(otu.hell ~ N+P+K+Ca+Mg+pH+Al+Fe+Mn+Zn+Mo, env.data.log)\notu.tab.1\nvif.cca(otu.tab.1)\n#\u8fdb\u4e00\u6b65\u7b5b\u9009\notu.tab.1 <- rda(otu.hell ~ N+P+K+Mg+pH+Al+Fe+Mn+Zn+Mo, env.data.log)\nvif.cca(otu.tab.1)\n#test again\notu.tab.1 <- rda(otu.hell ~ N+P+K+Mg+pH+Fe+Mn+Zn+Mo, env.data.log)\n#\u65b9\u5dee\u81a8\u80c0\u56e0\u5b50\u5206\u6790,\u76ee\u524d\u6240\u6709\u53d8\u91cf\u90fd\u5df2\u7ecf\u5c0f\u4e8e10\nvif.cca(otu.tab.1)\n##\u7528step\u6a21\u578b\u68c0\u6d4b\u6700\u4f4eAIC\u503c\n#mod.u <- step(otu.tab.0, scope = formula(otu.tab.1), test = \"perm\")\n# \"perm\"\u589e\u52a0P\u503c\u7b49\u53c2\u6570\nmod.d <- step(otu.tab.0, scope = (list(lower = formula(otu.tab.0), \n                                       upper = formula(otu.tab.1))))\nmod.d\n##\u672c\u5904\u7b5b\u9009\u7684\u7ed3\u679c\uff0c\u627e\u5230\u4e00\u4e2aMg\u73af\u5883\u56e0\u5b50\u9002\u5408\u6a21\u578b\u6784\u5efa\uff0c\u4e3a\u4e86\u4e0b\u4e00\u6b65\u753b\u56fe\uff0c\u6211\u4eec\n#\u4fdd\u7559\u4e0a\u8ff0\u7b5b\u9009\u7684\u7ed3\u679c\u6240\u6709\u975e\u5171\u7ebf\u6027\u7684\u73af\u5883\u56e0\u5b50\n#choose variables for best model and rda analysis again#\n(otu.rda.f <- rda(otu.hell ~ N+P+K+Mg+pH+Fe+Mn+Zn+Mo, env))\nanova(otu.rda.f)\nanova(otu.rda.f, by = \"term\")\nanova(otu.rda.f, by = \"axis\")\n#\u8ba1\u7b97db-rda\notu.tab.bray <- vegdist(otu.hell, \"bray\")\n#dbrda() for dbRDA # cca() for CCA\notu.tab.b<- capscale(otu.tab.bray ~ ., env)\n##\u7ed8\u5236db-RDA\u56fe\nplot(otu.tab.b ,display=c(\"si\",\"bp\",\"sp\"))\n\n## \u7528ggvegan\u7ed8\u5236RDA\u56fe\np1<- autoplot(otu.rda.f, arrows = TRUE,axes = c(1, 2), geom = \"text\", \n              layers = c(\"sites\"), \n              legend.position = \"right\", title = \"db-RDA\")\np1\np2<- autoplot(otu.rda.f, arrows = TRUE,axes = c(1, 2), geom = \"text\", \n              layers = c(\"biplot\"), \n              legend.position = \"right\", title = \"db-RDA\")\np2\np3<- autoplot(otu.rda.f, arrows = TRUE,axes = c(1, 2), geom = \"text\", \n              layers = c(\"sites\",\"biplot\"), \n              legend.position = \"right\", title = \"db-RDA\")\np3\np4<- autoplot(otu.rda.f, arrows = TRUE,axes = c(1, 2), \n              geom = \"text\", layers = c( \"species\",\"sites\", \n                                         \"biplot\", \"centroids\"), legend.position = \"right\", \n              title = \"db-RDA\")\np4\n## \u6dfb\u52a0\u56fe\u5c42\np5<-p4 + theme_bw()+theme(panel.grid.major =element_blank(), \n                          panel.grid.minor = element_blank())\np5\n###################<\/code><\/pre>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
###################\n\nlibrary(vegan)\nlibrary(tidyverse)\nlibrary(ggrepel)\n\ngenus <- read.delim(\"cpeBA.xls\",check.names=F,row.names = 1) %>%\n  t() %>%\n  decostand(method = \"hellinger\") %>% \n  as.data.frame() %>% \n  rownames_to_column(\"sample\") %>%\n  arrange(sample) %>% column_to_rownames(\"sample\")\n\nenvi <- read.delim(\"envir.xls\",check.names=F) %>%\n  arrange(sample) %>% \n  column_to_rownames(\"sample\") %>% \n  mutate(Depth=as.numeric(Depth)) %>% select(-Chla)\n\ngrp <- read.delim(\"group.txt\",check.names = F) \n\n\np <- rda(genus~.,data=envi) %>% summary()\np\nsp <- as.data.frame(p$species[,1:2])*2\nst <- as.data.frame(p$sites[,1:2])\nyz <- as.data.frame(p$biplot[,1:2]) \n\ncolors <- c(\"#FED439FF\",\"#709AE1FF\",\"#8A9197FF\",\"#D2AF81FF\",\n            \"#D5E4A2FF\",\"#197EC0FF\",\"#F05C3BFF\",\"#46732EFF\",\n            \"#71D0F5FF\",\"#370335FF\",\"#075149FF\",\"#C80813FF\",\n            \"#91331FFF\",\"#1A9993FF\",\"#FD8CC1FF\")\n\nggplot()+\n  geom_point(data= sp %>% rownames_to_column(\"genus\") %>% as.data.frame(),\n             aes(RDA1,RDA2,fill=genus,color=genus),\n             size=4,pch=21,color=\"white\")+\n  scale_fill_manual(values = colors)+\n  geom_text_repel(data = sp,aes(RDA1,RDA2),label=row.names(sp))+\n  geom_segment(data=yz,aes(x=0,y =0,xend=RDA1,yend=RDA2), \n               arrow=arrow(angle=30,length=unit(0.25,\"cm\"),\n                           type=\"closed\"),linetype=1,\n               size=0.8,colour = \"#00A08A\")+\n  geom_text_repel(data = yz,aes(RDA1,RDA2,label=row.names(yz)))+\n  labs(x=paste(\"RDA1(\",format(100*p$cont[[1]][2,1],\n                              digits=4),\"%)\",sep=\"\"),\n       y=paste(\"RDA2(\",format(100*p$cont[[1]][2,2], \n                              digits=4),\"%)\",sep=\"\"))+\n  geom_hline(yintercept=0,linetype=3,size=0.8)+ \n  geom_vline(xintercept=0,linetype=3,size=0.8)+\n  theme_bw()+\n  theme(panel.grid=element_blank(),\n        legend.title =element_blank(),\n        plot.margin=unit(c(0.3,0.3,0.3,0.3),units=,\"cm\"),\n        panel.grid.major=element_blank(), \n        panel.grid.minor=element_blank(),\n        panel.background=element_rect(colour=\"black\",\n                                      size=1,fill=\"white\"),\n        axis.title.y=element_text(size=11,color=\"black\",\n                                  margin=margin(r=8),face=\"plain\"),\n        axis.title.x=element_text(size=11,color=\"black\",\n                                  margin=margin(t=8),face=\"plain\"),\n        axis.text=element_text(size=11,color=\"black\",face=\"plain\"))+\n  annotate(\"text\",x=0.2,y=0.6,label=expression(NO[3]^-\"\"),\n           parse = TRUE,size=4,colour=\"black\")+\n  annotate(\"text\",x=0.2,y=0.4,label=expression(NH[4]^+\"\"),\n           parse = TRUE,size=4,colour=\"black\")+\n  annotate(\"text\",x=0.6,y=0.4,label=expression(SiO[4]^4^- \"\"),\n           parse = TRUE,size=4,colour=\"black\")+\n  annotate(\"text\",x=0.7,y=0.5,label=expression(PO[4]^3^- \"\"),\n           parse = TRUE,size=4,colour=\"black\")+\n  annotate(\"text\",x=-0.6,y=0.5,label=expression(NO[2]^- \"\"),\n           parse = TRUE,size=4,colour=\"black\")+\n  theme(legend.position=c(0,1),legend.justification = c(0,1),\n        legend.background=element_blank(),\n        legend.box.background=element_rect(colour = \"black\"),\n        legend.title = element_blank(), # \u56fe\u4f8b\u6807\u9898\u4e3a\u7a7a\n        legend.key=element_blank(),   # \u56fe\u4f8b\u952e\u4e3a\u7a7a\n        legend.text = element_text(color=\"black\",size=10), # \u5b9a\u4e49\u56fe\u4f8b\u6587\u672c\n        legend.spacing.x=unit(0,'cm'), # \u5b9a\u4e49\u6587\u672c\u4e66\u5e73\u8ddd\u79bb\n        legend.key.width=unit(0.5,'cm'), # \u5b9a\u4e49\u56fe\u4f8b\u6c34\u5e73\u5927\u5c0f\n        legend.key.height=unit(0.5,'cm'))<\/code><\/pre>\n\n\n\n
\"\"<\/div><\/figure>\n","protected":false},"excerpt":{"rendered":"","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1],"tags":[3],"class_list":["post-291","post","type-post","status-publish","format-standard","hentry","category-uncategorized","tag-r"],"_links":{"self":[{"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/posts\/291","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/linguopeng.top\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=291"}],"version-history":[{"count":2,"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/posts\/291\/revisions"}],"predecessor-version":[{"id":302,"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/posts\/291\/revisions\/302"}],"wp:attachment":[{"href":"https:\/\/linguopeng.top\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=291"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/linguopeng.top\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=291"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/linguopeng.top\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=291"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}