扩增子数据怎么去除样本污染？decontam和microDecon两款工具任您挑选！

创新多组学技术服务

2023.6.13

作者欧易生物

TA的动态

利用扩增子调查微生物群落，是微生物组研究中最基本的方法，而高通量测序技术不仅增强了检测物种的能力，同时也带来污染序列的风险。污染序列对低微生物量的样本产生不成比例的影响，并可能导致关于血液和身体组织等低微生物量环境中存在细菌的争议性说法[1]。

有研究表明，污染是无处不在的。即使在不应该含DNA/RNA的地方，污染也会发生[2]。污染一方面可以通过实验技术去除，另一方面也可以在生物信息分析中检测与去除。decontam[1]与microDecon[2]两款R语言软件包均可实现污染序列的检测，不同的是:decontam会直接去除被认为是污染的ASV/OTU（以下统称ASV）；microDecon只会从每个ASV中减去污染序列数。

decontam使用统计学方法检测污染序列并去除，因此使用decontam包进行去污染操作至少需要2个负对照样本。decontam包有两种运行模式，以下仅介绍其prevalence方法。它假设样本总DNA(T)是污染DNA(C)与样本真实DNA(S)的混合，即T = C + S，因此，在负对照样本中S~0。对于每个ASV，构造2×2列联表并行卡方检验，对于小样本则进行Fisher精确检验。计算得到的P值与用户自定义的P*作比较，如P < P* 则认为该ASV是污染。

使用decontam需要准备：（1）样本-分组文件；（2）ASV丰度表；（3）ASV分类信息表（如有）。

# 加载必要包library(phyloseq)library(decontam)
# 读取数据文件mapp <- read.table("样本-分组文件", sep="\\t", row.names=1, header=T, comment.char="")head(mapp)##              Group## Ind1   True Sample## Ind2   True Sample## Ind3   True Sample## Blank1    Negative## Blank2    Negative## Blank3    Negativeotu <- as.matrix(read.table("ASV丰度表", sep="\\t", row.name=1, header=T, check.names=F))head(otu)##   Ind1 Ind2 Ind3 Blank1 Blank2 Blank3## 1  240  194  186    182    206    203## 2  218  192  241    224    187    205## 3  170  213  203    210    208    206## 4  221  200  190    192    205    183## 5  214  186  186    204    188    195## 6  193  195  221    193    220    213## 如没有ASV分类信息表则不执行下一行代码tax <- as.matrix(read.table("ASV分类信息表", sep="\\t", row.names=1, header=T, check.names=F))head(tax)##   Ind1  Ind2   Ind3   Blank1 Blank2 Blank3 taxonomy                       ## 1  240   194    186      182    206    203   k__Bacteria;p__Proteobacteria## 2  218   192    241      224    187    205   k__Bacteria;p__Proteobacteria## 3  170   213    203      210    208    206   k__Bacteria;p__Proteobacteria## 4  221   200    190      192    205    183   k__Bacteria                  ## 5  214   186    186      204    188    195   k__Bacteria;p__Firmicutes    ## 6  193   195    221      193    220    213   k__Bacteria# 转为phyloseq对象## ASV丰度表中行为ASV，列为样本otu <- otu_table(otu, taxa_are_rows=T)## 如没有ASV分类信息表则不执行下一行代码tax <- tax_table(tax)mapp <- sample_data(mapp)ps <- phyloseq(otu, mapp, tax)## 如没有ASV分类信息表则改为执行如下代码ps <- phyloseq(otu, mapp)

以下使用decontam官方数据进行演示：

ps <- readRDS(system.file("extdata", "MUClite.rds", package="decontam"))# 有569个样本、1951个ASVps## phyloseq-class experiment-level object## otu_table()   OTU Table:         [ 1951 taxa and 569 samples ]## sample_data() Sample Data:       [ 569 samples by 6 sample variables ]## tax_table()   Taxonomy Table:    [ 1951 taxa by 6 taxonomic ranks ]# 标记负对照样本sample_data(ps)$is.neg <- sample_data(ps)$Sample_or_Control == "Control Sample"head(sample_data(ps))##                    X.SampleID PlateNumber Subject     Habitat quant_reading## P1101C01701R00 P1101C01701R00           3    1101      Tongue          2829## P1101C01702R00 P1101C01702R00           3    1101 Cheek_Right          5808## P1101C01703R00 P1101C01703R00           3    1101  Cheek_Left          5052## P1101C08701R00 P1101C08701R00           3    1101      Tongue          1227## P1101C08702R00 P1101C08702R00           3    1101 Cheek_Right          3872## P1101C08703R00 P1101C08703R00           3    1101  Cheek_Left          5872##                Sample_or_Control is.neg## P1101C01701R00       True Sample  FALSE## P1101C01702R00       True Sample  FALSE## P1101C01703R00       True Sample  FALSE## P1101C08701R00       True Sample  FALSE## P1101C08702R00       True Sample  FALSE## P1101C08703R00       True Sample  FALSE# 使用prevalence方法去污染，使用默认P*值=0.1contamdf.prev <- isContaminant(ps, method="prevalence", neg="is.neg", threshold=0.1)
# 有124个ASV被鉴定为污染table(contamdf.prev$contaminant)## ## FALSE  TRUE ##  1827   124# 剔除污染ASV，以便后续输出keep <- prune_taxa(!contamdf.prev$contaminant, ps)

microDecon作者McKnight认为，decontam存在删除研究系统中真实存在的低丰度ASV的风险。因此，microDecon是去除ASV中的污染序列数，它的假设条件是所有样本都会受到来自同一来源的相同比例的污染。microDecon去除污染步骤可简单归纳为4步：

（1）将数据子集化为仅包括污染ASV（在负对照中出现的ASV）；

（2）估计重叠ASV的数量（在负对照和样本中均出现的ASV），并使用该估计来确定要用作constant的最佳ASV（使用的算法基于通过多次模拟开发的回归方程）；
（3）将负对照中每个ASV的read数除以负对照中constant的reads数；
（4）将这些值乘以样本中constant的reads数。

microDecon文章中去污染算法示意图：

使用microDecon需要准备：ASV丰度表。以下使用microDecon官方数据进行演示：

# 加载microDecon包library(microDecon)
# ASV丰度表格式要求## 第一列为ASV ID## 第二列为负对照## 负对照后一列为样本## 如有分类信息则最后一列必须为分类信息Example_1##   OTU_ID Blank1 Blank2 Blank3 Ind1 Ind2 Ind3                         taxa## 1   OTU1      0      0      0   60   64   40 K_Bacteria; P_Actinobacteria## 2   OTU2    200    220    180  660  520  480 K_Bacteria; P_Proteobacteria## 3   OTU3   1000    800   1300 1440 1000  700 K_Bacteria; P_Proteobacteria## 4   OTU4     50     30     70   70   48   35  K_Bacteria; P_Bacteroidetes## 5   OTU5      0      0      0 2400 1900 2100                   K_Bacteria## 6   OTU6     25     10     30   30   20   15                   K_Bacteria# 去污染主函数## 需要显式指定的参数## numb.blanks —— 负对照的数量## numb.ind —— 向量，指定每个自定义组内的样本数，Example_1中样本只有1组3个样本## numb.tax —— 布尔值，数据中是否含有分类信息result <- decon(data=Example_1, numb.blanks=3, numb.ind=3, taxa=T)
# 结果为1个列表，包含5个数据框## 去除污染reads数后的数据为result$decon.table##   OTU_ID Mean.blank Ind1 Ind2 Ind3                         taxa## 1   OTU1     0.0000   60   64   40 K_Bacteria; P_Actinobacteria## 2   OTU2   208.0753  359  318  329 K_Bacteria; P_Proteobacteria## 5   OTU5     0.0000 2400 1900 2100                   K_Bacteria

参考文献

[1] DAVIS N M, PROCTOR D, HOLMES S P, et al. Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data[J]. Microbiome, 2017.

[2] MCKNIGHT D T, HUERLIMANN R, BOWER D S, et al. microDecon: A highly accurate read‐subtraction tool for the post‐sequencing removal of contamination in metabarcoding studies[J]. Environmental DNA, 2019, 1(1).