当前位置:网站首页>Microbial ecological sequencing analysis -- CCA analysis

Microbial ecological sequencing analysis -- CCA analysis

2022-06-21 06:02:00 Silent voice

Microbial ecological sequencing analysis ——CCA analysis

library(vegan)
library(ggplot2)
library(permute)
library(lattice)
sa4=read.table("spes.csv",header=T,row.names=1,sep=",")
S8.p=read.table("S8.p.csv",header=T,sep=",",row.names=1)# Select physical and chemical factors with significant correlation 
group <- read.table("name.csv", header=F,sep=",",colClasses=c("character","character"))# Used to generate legend , Classify the samples 
cca<-cca(sa4,S8.p, scale=T)
sam <- data.frame(cca$CCA$u[,c(1,2)], group$V2)  # Extract sample scores 
colnames(sam) <- c("CCA1","CCA2","group")
spec <- cca$CCA$v[,c(1,2)]  # Species score 
spec <- as.data.frame(spec)
# It can be used to display species data in the figure 
#spec<-data.frame(spece=row.names(cca$CCA$v),RDA1=cca$CCA$v[,1],RDA2=cca$CCA$v[,2])  
env <- cca$CCA$biplot[,c(1,2)] # Environmental factor score 
env <- as.data.frame(env)
cca1 =round(cca$CCA$eig[1]/sum(cca$CCA$eig)*100,2) # First axis label 
cca2 =round(cca$CCA$eig[2]/sum(cca$CCA$eig)*100,2) # Second axis label 
# Creation of drawing objects 
p <- ggplot(data=sam,aes(CCA1,CCA2))
# Geometry object 
p <- p + geom_point(aes(colour=group,shape=group),size=5) + 
# Show the species in the figure 
#geom_point(data=spec,aes(shape=spece),size=2) + scale_shape_manual(values=seq(0,20))+
  # Add sample label   
	#geom_text(aes(label=rownames(sam)),
     #         size=4,hjust=0.5,vjust=-0.7,position = "jitter") + 
    scale_shape_manual(values=c(19,19,19,19,19)) + 
    labs(title="CCA Plot",x=paste("CCA1 ",cca1," %"),y=paste("CCA2 ",cca2," %")) +
    theme(text=element_text(family="serif"))
# Remove background and gridlines 
p=p+theme_bw() + theme(panel.grid=element_blank()) 
p=p + geom_segment(data = env,aes(x=0,y=0,xend = env[,1], yend = env[,2]), colour="purple", size=1.5, arrow=arrow(angle=35, length=unit(0.3, "cm"))) +
    geom_text(data=env, aes(x=env[,1], y=env[,2], label=rownames(env)),
              size=5, colour="purple",hjust = (1 - 2 * sign(env[ ,1])) / 3,
              angle = (180/pi) * atan(env[ ,2]/env[ ,1]))

I see that many students ask questions about the imported data format , I wrote a detailed ranking analysis article again , It was posted on the official account “ When scientific research is in progress ” On , Article title :R mapping -RDA Sorting analysis , There are screenshots of data format , You can scan the QR code and pay attention to view .
 Insert picture description here

原网站

版权声明
本文为[Silent voice]所创,转载请带上原文链接,感谢
https://yzsam.com/2022/172/202206210550435216.html