R packages：phyloseq提取某一ASV或OTU

导入数据

library("ggpubr")#用于事后检验标记
library("MicrobiotaProcess")
library("ggsci")
otu <- "1_table-dada2_pe.qza"
rep <- "2_rep-seqs-dada2_pe.qza"
tree <- "3_rooted-tree_pe.qza"
tax <- "4_taxonomy_pe.qza"
sample <- "5_metadata.txt"
qiimedata <- import_qiime2(otuqza=otu,taxaqza=tax,refseqqza=rep,
                           mapfilename=sample,treeqza=tree)
library(phyloseq)
set.seed(100)#设置一个随机种子便于重复
qiimedata <-  rarefy_even_depth(qiimedata, rngseed=1, 
                                sample.size=0.95*min(sample_sums(qiimedata)), 
                                replace=F)

统计所有ASV的序列长度

stat <-data.frame(qiimedata@refseq)[,1]%>%nchar()
plot(table(stat),xlab = "Length/bp",
     ylab = "Count")

Rplot.png

ASV长度和在样本中总Count的关系

taxa_sums<-taxa_sums(qiimedata)%>%data.frame()
plot(stat,taxa_sums[,1],xlab = "Length/bp",ylab = "Count")

Rplot01.png

提取名字前缀为f69045904的序列的完整名字

otu <- qiimedata@refseq %>% data.frame() %>% row.names()
name <- otu[grep("^f69045904",otu)]
name

#[1] "f69045904699eefc4f5f8b0c3f809d7e"

提取这条序列

seq <- data.frame(qiimedata@refseq)[grep("^ef57f22f3235b2585",otu),]
seq

#[1] "TGGGGGATATTGCACAATGGGGGGAACCCTGATGCAGCGACGCCGCGTGGGTGAAGAAGCGCCTCGGCGCGTAAAGCCCTGTCAGCAGGGAAGAAAATGACGGTACCTGAAGAAGAAGCCCCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGGGCAAGCGTTATCCGGATTTACTGGGTGTAAAGGGGGCGCAGACGGCGATGCAAGCCAGGAGTGAAAGCCCGGGGCCCAACCCCGGGACTGCTCTTGGAACTGCGTGGCTGGAGTGCAGGAGGGGCAGGCGGAATTCCTGGTGTAGCGGTGAAATGCGTAGATATCAGGAGGAACACCGGTGGCGAAGGCGGCCTGCTGGACTGCAACTGACGTTGAGGCCCGAAAGCGTGGGGAGCAAA"

计算这条序列的长度

len <- data.frame(qiimedata@refseq)["f69045904699eefc4f5f8b0c3f809d7e",]%>%nchar()
len <-seq %>%nchar()
len

#[1] 400

提取这条序列的物种分类学信息

tax <- data.frame(qiimedata@tax_table)[grep("^740e5",otu),]
tax

提取这条序列在所有样本中的丰度信息

testdata = subset_taxa(qiimedata,qiimedata@refseq %>% 
                         data.frame() %>% 
                         row.names()%in%c("ef57f22f3235b25850f7e57adf243a6e",
                                          "b258b594ebabb590b38f9f57ca23cfce",
                                          "29795bb279254860b9b8484fa2a7be60"
                                          ))
p <-  plot_bar(testdata,fill = "Kingdom")+
  facet_wrap(~group,scale="free_x")+
  scale_fill_aaas()+xlab(NULL)
p

plot_bar.png

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

p

ggplot2.png