让我对phyloseq欲罢不能的10个理由
只有一个原始数据:文末获得:
理由1:--一条函数载入微生物组分析的全部文件
如何提取phyloseq中的信息
理由2:一条代码过滤掉序列数少的样本
理由3:一条代码:根据分组文件的分组或者ID列,或者任何一列过滤样本。
理由4:一条代码提取序列数大于某个阈值的otu
理由5:一条代码根据最小序列数抽平
理由6 一条代码标准化
理由7:一条代码提取目标OTU
理由8:一条代码计算样本序列数量
理由9:一条代码计算超过20中距离
理由10:一条代码计算七:种排序
欢迎加入微生信生物
这就是你认识的微生信生物
理由1:—一条函数载入微生物组分析的全部文件
将otu表格,注释文件,分组文件,进化树,代表序列文件都可以保存在一个rds文件中,并且作为一个整体来活动。
给大家举个例子:
让我们过滤到了某个投,则其相应的注释信息和进化树中的位置都会消失,代表序列也会消失,这是连带处理,非常方便。
###导入原始ps文件
ps = readRDS("./ps_liu.rds")
如何提取phyloseq中的信息
otu = as.data.frame(otu_table(ps))
tax = as.data.frame(tax_table(ps))
map = as.data.frame(sample_data(ps))
write.table("ID\t", file="otu.txt",append = FALSE, quote = FALSE, sep="\t",eol = "", na = "NA", dec = ".", row.names = F,col.names = F)
write.table(otu, file="otu.txt",append = T, quote = FALSE, sep="\t",eol = "\n", na = "NA", dec = ".", row.names = TRUE,col.names = TRUE)
write.table("ID\t", file="tax.txt",append = FALSE, quote = FALSE, sep="\t",eol = "", na = "NA", dec = ".", row.names = F,col.names = F)
write.table(tax, file="tax.txt",append = T, quote = FALSE, sep="\t",eol = "\n", na = "NA", dec = ".", row.names = TRUE,col.names = TRUE)
write.table("ID\t", file="map.txt",append = FALSE, quote = FALSE, sep="\t",eol = "", na = "NA", dec = ".", row.names = F,col.names = F)
write.table(map, file="map.txt",append = T, quote = FALSE, sep="\t",eol = "\n", na = "NA", dec = ".", row.names = TRUE,col.names = TRUE)
理由2:一条代码过滤掉序列数少的样本
##选择count数量大于2000的样品,也就是说将序列数量过少的样品去除,一般不操作,我们肯定弄好才分析。# ps1 <- prune_samples(sample_sums(ps) >=2000,ps1);ps0
理由3:一条代码:根据分组文件的分组或者ID列,或者任何一列过滤样本。
ps_sub <- subset_samples(ps,SampleType %in% c("OF"));ps_sub
ps_sub <- subset_samples(ps,ID %in% c("SRR8247296"));ps_sub
理由4:一条代码提取序列数大于某个阈值的otu
ps_sub = filter_taxa(ps, function(x) sum(x ) > 1 , TRUE);ps_sub #筛选序列数量大于1的
#是否需要保存
# saveRDS(ps_sub,"./ps_sub.rds")
理由5:一条代码根据最小序列数抽平
当然可以根据自己需求来指定抽平数量
#
set.seed(72) # setting seed for reproducibility
bulk.physeq.rare = rarefy_even_depth(ps)
理由6 一条代码标准化
# 相对丰度标准化编号:ps1_rela
ps1_rela = transform_sample_counts(ps, function(x) x / sum(x) );ps1_rela
saveRDS(ps1_rela, "./ps_rela.rds")
#对数标准化编号:ps1_log(添加1是为了方式0的影响)
ps1_log = transform_sample_counts(ps, function(x) log(1 + x) );ps1_log
saveRDS(ps1_log, "./ps_log.rds")
理由7:一条代码提取目标OTU
ps2 <- ps %>%
subset_taxa(
Kingdom == "Bacteria" #&
# Genus == "Fusarium"
# Species %in%c("Fusarium_oxysporum","Fusarium_keratoplasticum")
#row.names(tax_table(ps1_rela ))%in%c("SH010924.07FU_KF986690_reps_singleton","SH020983.07FU_JN235282_refs")
)
ps2
理由8:一条代码计算样本序列数量
sample_sum(ps)
理由9:一条代码计算超过20中距离
distance :指定需要的距离
ordinate(ps, distance=dist)
理由10:一条代码计算七种排序
ordinate(ps, method="PCoA", distance=dist)