1.导入数据

library("MicrobiotaProcess")
library("phyloseq")
otu <- "./table-dada2_pe.qza"
rep <- "rep-seqs-dada2_pe.qza"
tree <- "rooted-tree_pe.qza"
tax <- "taxonomy_pe.qza"
sample <- "metadata_pe.txt"
qiimedata <- import_qiime2(otuqza=otu, taxaqza=tax,refseqqza=rep,
                          mapfilename=sample,treeqza=tree)

2.抽平

set.seed(100)#设置一个随机种子便于重复
qiimedata <-  rarefy_even_depth(qiimedata, rngseed=1, 
                                sample.size=0.95*min(sample_sums(qiimedata)), 
                                replace=F)
qiimedata@tax_table %>% data.frame() %>% View()

3. 提取特定属Genus的tax

3.1 乳酸菌

Lactobacillus = subset_taxa(rare.data, Genus == "Lactobacillus")
sample_sums(Lactobacillus)
#SRR10097305 SRR10097306 SRR10097307 SRR10097308 SRR10097309 
#          0           0           0           0           0 
#SRR10097310 SRR10097311 SRR10097312 SRR10097313 SRR10097314 
#          6          65         194           0          44 
#SRR10097315 SRR10097316 SRR10097317 SRR10097318 SRR10097319 
#        142           9          41           0           0 
#SRR10097320 SRR10097321 SRR10097323 SRR10097324 SRR10097325 
#         13          98          46           0           0 
#SRR10097326 SRR10097327 SRR10097328 SRR10097329 SRR10097330 
#          0           0           0           0           0

可视化

p <-  plot_bar(Lactobacillus, fill="Species")+
  facet_wrap(~group,scale="free_x")+
  scale_fill_aaas()+xlab(NULL)
p

plot_bar.png

plot_bar这个函数虽然可以绘制每个样品中物种丰度，但是不能进行组间显著性检验。

方法二：合并相同的Genus

psdata = tax_glom(qiimedata,taxrank = "Genus")
p = plot_bar(psdata, fill="Genus")
p$data
p$data$group <- factor(p$data$group,levels = c("Pre","Post"))

提取特定的Genus

data <- subset(p$data,Genus %in%c("Lactobacillus"))
p<-ggplot(data,aes(x=group,y = Abundance , fill = group)) +
  geom_boxplot(alpha = 0.9)+
  xlab(NULL)+
  ylab("Abundance")+theme_bw()+
  facet_wrap(~Genus,scales= "free",nrow = 1)+#free每个画框自由缩放
  scale_fill_aaas()
p
mycompare=list(c("Pre","Post"))
p<-p + stat_compare_means(comparisons=mycompare,
                          label = "p.signif",
                          method = 'wilcox')
p

Lactobacillus.png

提取多个Genus

data <- subset(p$data,Genus %in%c("Lactobacillus","Bacillus",
                                  "Streptococcus", "Lactococcus",
                                  "Enterococcus","Clostridium",
                                  "Bifidobacterium","Bacteroides"))

多个Genus.png