「微生物组」基于LorMe包进行微生物组数据分析—微生物群落结构组成分析！！！

###安装LorMe包install.packages("LorMe")##加载MicrobiomeStat包library(LorMe)# Lightening One-Code Resolving Microbial Ecology Solutionlibrary(ggplot2)# Create Elegant Data Visualisations Using the Grammar of Graphicslibrary(dplyr)# A Grammar of Data Manipulation

##加载示例数据集data(testotu)# Load the standard Qiime output feature tablehead(testotu)# First column -ID, last column -taxonomic annotation, others - the featuretable.feature_table<- testotu[, -c(1,22)]tax_anno<- testotu[, c(1,22)]
#创建分组数据groupinformation <- data.frame(group = c(rep("a", 10), rep("b", 10)),factor1 = rnorm(10),factor2 = rnorm(mean = 100, 10),  subject = factor(c(1:10, 1:10)), group2 = c(rep("e", 5), rep("f", 5), rep("e", 5), rep("f", 5)))
###对数据进行封装test_object <- tax_summary(groupfile = groupinformation,inputtable = feature_table,reads = TRUE,taxonomytable = tax_anno)
###对一些默认设置进行配置my_col=color_scheme("Plan1")my_order=c("b","a")my_facet_order=c("e","f")test_object_plan2 <- object_config(taxobj = test_object,treat_location = 1,rep_location = 4,facet_location = 5,  subject_location = NULL, treat_col = my_col,treat_order = my_order,  facet_order = my_facet_order)

#Taxonomy tidying...#Genarating tables...#Done##Contained elements:#[1]"Groupfile""Base""Base_percent"#[4]"Base_taxonomy""Phylum""Phylum_percent"#[7]"Phylum_taxonomy""Genus""Genus_percent"#[10]"Genus_taxonomy""parameters"#Warning message:#Expected 7 pieces. Missing pieces filled with `NA`in261 rows [4,#6, 9, 12, 18, 20, 24, 32, 37, 43, 44, 45, 46, 58, 59, 68, 70, 71,#72, 76, ...].

##绘制门水平排名前10的物种require(magrittr)phylum10 <- community_plot(taxobj = Two_group,taxlevel ="Phylum",#分析水平设置n =10,#保留排名前10的物种rmprefix ="p__",#去除一些物种名中附带的前缀palette ="Paired")

##柱状堆积图绘制phylum10$barplotphylum10$mean_barplot##面积图phylum10$areaplot##冲积图phylum10$alluvialplot

genus10 <- community_plot(  taxobj = Two_group,  taxlevel = "Genus",  n = 10,  rmprefix = "g__",  palette = "Paired")##柱状堆积图绘制genus10$barplotgenus10$mean_barplot##面积图genus10$areaplot##冲积图genus10$alluvialplot

##提取数据中属水平物种排名前10的相对丰度数据data <- genus10$Grouped_Top10Genusdata$SampleID <- factor(data$SampleID,levels = data$SampleID[1:16])##颜色col <- genus10$filled_color#绘图ggplot(data, aes( x = SampleID, y=100 * rel_abundance,fill = tax))+  geom_bar(stat="summary",fun=mean,position = 'stack')+  scale_y_continuous(expand = c(0,0))+# 调整y轴属性，使柱子与X轴坐标接触  scale_x_discrete(expand = c(0,0))+# 调整y轴属性，使柱子与X轴坐标接触  labs(x=NULL,y="Relative Abundance(%)",#设置X轴和Y轴的名称以及添加标题       fill="Genus")+  guides(fill=guide_legend(keywidth = 1, keyheight = 1)) +#修改图例的框大小  theme_bw()+  theme(panel.grid = element_blank(),        axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),        axis.title.y = element_text(size = 12),        axis.text = element_text(size = 9))+  scale_fill_manual(values = col)+  facet_grid(~Group, scales = 'free',space = "free")

##合并ggplot(data, aes( x = Group, y=100 * rel_abundance,fill = tax))+  geom_bar(stat="summary",fun=mean,position = 'stack')+  scale_y_continuous(expand = c(0,0))+# 调整y轴属性，使柱子与X轴坐标接触  scale_x_discrete(expand = c(0,0))+# 调整y轴属性，使柱子与X轴坐标接触  labs(x=NULL,y="Relative Abundance(%)",#设置X轴和Y轴的名称以及添加标题       fill="Genus")+  guides(fill=guide_legend(keywidth = 1, keyheight = 1)) +#修改图例的框大小  theme_bw()+  theme(panel.grid = element_blank(),        axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),        axis.title.y = element_text(size = 12),        axis.text = element_text(size = 9))+  scale_fill_manual(values = col)