“码”辅导——微生物样本有其他来源的污染怎么办？Decontam开源包帮你忙！

2022-06-30 19:42:27, 星标关注上海欧易生物医学科技有限公司

编者按

想从环境中获取完全没有污染的微生物样本基本上是不可能实现的，不规范的实验操作又容易引入外源污染物，目前控制污染的策略基本是从实验角度出发，如更严格的采样方法、更完善的实验环境、实验操作等。但这些需要耗费更多的精力、更高的成本，同时只能适当地降低污染影响。那还有什么其他方案呢？

接下来轮到Decontam开源包登场（Benjamin Callahan, Nicole Davis; 2018）。

Decontam是一种算法，使用分类器根据样本的模式对污染及非污染进行鉴别。Decontam包基于R语言开发，使用简单，运行快速，污染消除效果好。目前这个包提供了两种分类模式：一是基于样本浓度；二是基于阴性对照。其中阴性对照运用较为方便，本次分享采用这种模式。

不废话，直接上硬菜！我们一起来看看怎么用吧！

备菜

环境：R（推荐安装RStudio）、phyloseq包、decontam包

输入文件：OTU丰度表、OTU分类信息表、对照样本信息表

OTUs.xls OTU丰度表，每行是OTU ID，每列是实验样本的OTU reads数

OTUs_tax.xls OTU分类信息表，丰度表后多一列各OTU对应的物种分类信息

mapping.txt 对照样本信息表，每行是Sample ID，每列是样本的对应信息

（阴性对照样本是Control）

开整

# 安装载入R包source("https://bioconductor.org/biocLite.R")biocLite("phyloseq")library(phyloseq)install.packages("BiocManager")BiocManager::install("decontam")library(decontam)

# 读取表格数据，第一行是表头，第一列是行名。丰度表和分类信息表转化成矩阵形式otu_p <- as.matrix(read.table("分析的文件路径/OTUs.xls", header = TRUE, sep = "\\t", row.names = 1))tax_p <- as.matrix(read.table("分析的文件路径/OTUs_tax.xls", header = TRUE, sep = "\\t", row.names = 1))sammapp <- read.table("分析的文件路径/mapping.txt",sep="\\t", header = TRUE, row.names = 1)

# 读入数据整合成 phyloseq 对象，也就是 去污染时需要的数据格式# 生成otu_table类型的对象，OTU 在每个样本中的丰度信息OTU <- otu_table(otu_p, taxa_are_rows = TRUE)# taxonomyTable，OTU 的物种信息，行名必须匹配otu_table OTU名TAX <- tax_table(tax_p)# sample_data，所有样本的的表型信息，行名必须匹配otu_table 样本名sammapp <- sample_data(sammapp)# 将数据组合成一个 phyloseq 对象，isContaminant的输入格式physeq_p <- phyloseq(OTU, TAX, sammapp)

# 根据对照样本信息，生成名为is.neg的列sample_data(physeq_p)$is.neg <- sample_data(physeq_p)$Group == "Control"

# 进行污染判定，阈值可通过threshold调整contamdf.prev <- isContaminant(physeq_p, method="prevalence", neg="is.neg")# contamdf.prev <- isContaminant(physeq_p, method="prevalence", neg="is.neg",threshold=0.1)

# 查看污染对应的OTU IDphyseq1.contam <- prune_taxa(contamdf.prev$contaminant, physeq_p)rownames(otu_table(physeq1.contam))# 去除污染对应的OTU IDphyseq1.noncontam <- prune_taxa(!contamdf.prev$contaminant, physeq_p)

# 准备输出结果# 提取OTU丰度表otu_final = otu_table(physeq1.noncontam)# OTU丰度表，第一列给列名 OTU_IDotu_final = as.data.frame(otu_final)otu_colname = colnames(otu_final)otu_final = cbind(rownames(otu_final),otu_final)colnames(otu_final) = c("OTU_ID", otu_colname)
# 提取OTU分类信息表tax_final = tax_table(physeq1.noncontam)# OTU分类信息表，第一列给列名 OTU_IDtax_final = as.data.frame(tax_final)tax_colname = colnames(tax_final)tax_final = cbind(rownames(tax_final),tax_final)colnames(tax_final) = c("OTU_ID", tax_colname)

# 输出去污染后的OTU结果write.table(otu_final, "输出文件路径/OTUs_out.xls", sep = "\\t", quote = F, col.names = T, row.names = F)write.table(tax_final, "输出文件路径/OTUs_tax_out.xls", sep = "\\t", quote = F, col.names = T, row.names = F)