让我对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)

欢迎加入微生信生物

这就是你认识的微生信生物

微生信生物

赞 (0)

相关推荐