引言
????在對(duì)樣本的微生物群落組成多樣性分析過程中，除了Alpha多樣性與Beta多樣性的分析餐抢，最重要的一個(gè)部分就是對(duì)樣本中微生物具體的群落結(jié)構(gòu)組成情況進(jìn)行分析皇筛。一般我們從門水平及屬水平上計(jì)算優(yōu)勢(shì)物種的相對(duì)豐度并通過表格或者柱狀堆積圖進(jìn)行展示属拾。
正文—代碼
設(shè)置工作目錄并加載包
rm(list=ls())#clear Global Environment
setwd('D:\\桌面\\物種豐度計(jì)算及可視化')#設(shè)置工作路徑
#安裝包
install.packages("reshape2")
install.packages("ggplot2")
install.packages("ggprism")
install.packages("plyr")
#加載R包
library (reshape2)
library(ggplot2)
library(ggprism)
library(plyr)
門水平的物種豐度計(jì)算及可視化
1桃移、讀取數(shù)據(jù)
df1 <- read.table(file="Phylum.txt",sep="\t",header=T,check.names=FALSE)
2、對(duì)數(shù)據(jù)進(jìn)行處理
##利用循環(huán)將其中的重復(fù)門數(shù)據(jù)進(jìn)行求和
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蹂析、計(jì)算物種總豐度并進(jìn)行降序排列
df3$rowsum <- apply(df3,1,sum)
df4 <- df3[order (df3$rowsum,decreasing=TRUE),]
df5 = df4[,-6]
#求物種相對(duì)豐度
df6 <- apply(df5,2,function(x) x/sum(x))
#導(dǎo)出數(shù)據(jù)
write.table (df9, file ="phylum.csv",sep =",", quote =FALSE)
#變量格式轉(zhuǎn)換,寬數(shù)據(jù)轉(zhuǎn)化為長(zhǎng)數(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
屬水平的物種豐度計(jì)算及可視化
1、加載數(shù)據(jù)
df1 <- read.table(file="Genus.txt",sep="\t",header=T,check.names=FALSE)
2电抚、對(duì)數(shù)據(jù)進(jìn)行處理
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惕稻、計(jì)算物種豐度并降序排列
df3$rowsum <- apply(df3,1,sum)
df4 <- df3[order (df3$rowsum,decreasing=TRUE),]
df5 = df4[,-6]
#求物種相對(duì)豐度
df6 <- apply(df5,2,function(x) x/sum(x))
#由于之間已經(jīng)按照每行的和進(jìn)行過升序排列，所以可以直接取前10行
df7 <- df6[1:10,]
df8 <- 1-apply(df7, 2, sum) #計(jì)算剩下物種的總豐度
#合并數(shù)據(jù)
df9 <- rbind(df7,df8)
row.names(df9)[11]="Others"
#導(dǎo)出數(shù)據(jù)
write.table (df9, file ="genus.csv",sep =",", quote =FALSE)
#變量格式轉(zhuǎn)換,寬數(shù)據(jù)轉(zhuǎn)化為長(zhǎng)數(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ù)及源碼請(qǐng)?jiān)谖⑿殴娞?hào)后臺(tái)回復(fù)Re_abundance獲人跻恕！
