引言

????在對樣本的微生物群落組成多樣性分析過程中，除了Alpha多樣性與Beta多樣性的分析，最重要的一個部分就是對樣本中微生物具體的群落結(jié)構(gòu)組成情況進行分析。一般我們從門水平及屬水平上計算優(yōu)勢物種的相對豐度并通過表格或者柱狀堆積圖進行展示。

正文—代碼

設置工作目錄并加載包

rm(list=ls())#clear Global Environment
setwd('D:\\桌面\\物種豐度計算及可視化')#設置工作路徑
#安裝包
install.packages("reshape2")
install.packages("ggplot2")
install.packages("ggprism")
install.packages("plyr")
#加載R包
library (reshape2)
library(ggplot2)
library(ggprism)
library(plyr)

門水平的物種豐度計算及可視化

1、讀取數(shù)據(jù)

df1 <- read.table(file="Phylum.txt",sep="\t",header=T,check.names=FALSE)

2、對數(shù)據(jù)進行處理

##利用循環(huán)將其中的重復門數(shù)據(jù)進行求和
data<-aggregate(E ~ Tax,data=df1,sum)
colnames(data)[2]<-"example"
for (i in colnames(df1)[2:length(colnames(df1))]){
  data1<-aggregate(df1[,i]~Tax,data=df1,sum)
  colnames(data1)[2]<-i  
  data<-merge(data,data1,by="Tax")
}
df2<-data[,-2]
rownames(df2)=df2$Tax
df3=df2[,-1]

3、計算物種總豐度并進行降序排列

df3$rowsum <- apply(df3,1,sum)
df4 <- df3[order (df3$rowsum,decreasing=TRUE),]
df5 = df4[,-6]
#求物種相對豐度
df6 <- apply(df5,2,function(x) x/sum(x))
#導出數(shù)據(jù)
write.table (df9, file ="phylum.csv",sep =",", quote =FALSE)
#變量格式轉(zhuǎn)換,寬數(shù)據(jù)轉(zhuǎn)化為長數(shù)據(jù),方便后續(xù)作圖
df_phylum <- melt(df6)
names(df_phylum)[1:2] <- c("Taxonomy","sample")

4、繪圖

p1 <- ggplot(df_phylum, aes( x = sample,y=100 * value,fill = Taxonomy))+
  geom_col(position = 'stack', width = 0.6)+#geom_bar(position = "stack", stat = "identity", width = 0.6) 
  scale_y_continuous(expand = c(0,0))+# 
  labs(x="Samples",y="Relative Abundance(%)",
       fill="Taxonomy")+
  guides(fill=guide_legend(keywidth = 1, keyheight = 1)) +
  theme_prism(palette = "candy_bright",
              base_fontface = "plain", 
              base_family = "serif",
              base_size = 16,
              base_line_size = 0.8, 
              axis_text_angle = 45)+ 
  scale_fill_prism(palette = "candy_bright")
p1

image.png

屬水平的物種豐度計算及可視化

1、加載數(shù)據(jù)

df1 <- read.table(file="Genus.txt",sep="\t",header=T,check.names=FALSE)

2、對數(shù)據(jù)進行處理

data<-aggregate(E ~ Tax,data=df1,sum)
colnames(data)[2]<-"example"
for (i in colnames(df1)[2:length(colnames(df1))]){
  data1<-aggregate(df1[,i]~Tax,data=df1,sum)
  colnames(data1)[2]<-i  
  data<-merge(data,data1,by="Tax")
  }
df2<-data[,-2]
rownames(df2)=df2$Tax
df3=df2[,-1]

3、計算物種豐度并降序排列

df3$rowsum <- apply(df3,1,sum)
df4 <- df3[order (df3$rowsum,decreasing=TRUE),]
df5 = df4[,-6]
#求物種相對豐度
df6 <- apply(df5,2,function(x) x/sum(x))
#由于之間已經(jīng)按照每行的和進行過升序排列，所以可以直接取前10行
df7 <-  df6[1:10,]
df8 <- 1-apply(df7, 2, sum) #計算剩下物種的總豐度
#合并數(shù)據(jù)
df9 <- rbind(df7,df8)
row.names(df9)[11]="Others"
#導出數(shù)據(jù)
write.table (df9, file ="genus.csv",sep =",", quote =FALSE)
#變量格式轉(zhuǎn)換,寬數(shù)據(jù)轉(zhuǎn)化為長數(shù)據(jù),方便后續(xù)作圖
df_genus <- melt(df9)
names(df_genus)[1:2] <- c("Taxonomy","sample")

4、繪圖

p2 <- ggplot(df_genus, aes( x = sample,y=100 * value,fill = Taxonomy))+
  geom_col(position = 'stack', width = 0.6)+#geom_bar(position = "stack", stat = "identity", width = 0.6) 
  scale_y_continuous(expand = c(0,0))+
  labs(x="Samples",y="Relative Abundance(%)",
       fill="Taxonomy")+
  guides(fill=guide_legend(keywidth = 1, keyheight = 1)) +
  theme_prism(palette = "candy_bright",
              base_fontface = "plain", 
              base_family = "serif", 
              base_size = 16, 
              base_line_size = 0.8,
              axis_text_angle = 45)+ 
  scale_fill_prism(palette = "colors")
p2

image.png

拼圖

library("cowplot")
plot_grid(p1,p2, labels=c('A','B'), ncol=1, nrow=2)

image.png

作圖數(shù)據(jù)及源碼請在微信公眾號后臺回復Re_abundance獲取！