前言
实验过程中，背景细菌的存在可能导致交叉污染，引起假阳性的发生，从而给mNGS结果解读带来困扰。mNGS专家共识中指出，实验室应该建立背景细菌库并定时更新，用于减少假阳性。
decontam通过浓度梯度（frequency）或阴性对照的方式（Prevalence）鉴定实验过程中的背景菌。因此本文旨在介绍decontam内部细节。
正文开始前，有如下注意事项：
本文主要集中于decontam原理解读，因此使用方法和结果解读不过多留恋，读者请参考官方教程。
下文只涉及官网教程中提到的两个算法（frequency和prevalence）的核心部分，细枝末节和其他算法略过不提。
案例数据和代码均摘自作者源码，同时结合自身理解。本人非数学科班出身，若理解有误，欢迎大家批评指正。
基本使用
suppressPackageStartupMessages({
library(rio)
library(ggplot2)
library(dplyr)
library(decontam)
})
使用官方数据：
# load data
ps <- readRDS(system.file("extdata", "MUClite.rds", package="decontam"))
head(data.frame(ps@otu_table)[1:5, 1:5])
Seq1 Seq2 Seq3 Seq4 Seq5
P1101C01701R00 3502 8391 0 193 2838
P1101C01702R00 12040 152 0 7924 0
P1101C01703R00 9877 2401 0 4333 0
P1101C08701R00 4035 5706 0 257 3161
P1101C08702R00 12491 1444 0 4384 294
行为样本，列为不同的变量或病原。
# calculate contaminant
contamdf.freq <- isContaminant(ps, method="frequency", conc="quant_reading")
head(contamdf.freq)
freq prev p.freq p.prev p contaminant
Seq1 0.323002694 549 1.000000e+00 NA 1.000000e+00 FALSE
Seq2 0.098667396 538 1.000000e+00 NA 1.000000e+00 FALSE
Seq3 0.003551358 160 1.135975e-18 NA 1.135975e-18 TRUE
Seq4 0.067588419 519 9.999998e-01 NA 9.999998e-01 FALSE
Seq5 0.045174743 354 1.000000e+00 NA 1.000000e+00 FALSE
Seq6 0.040417101 538 1.000000e+00 NA 1.000000e+00 FALSE
结果最后一列即判定该变量或病原是否为背景菌。
详情见官方教程.
两个算法
frequency
思路如下：
随着样本浓度的增加，文库中的背景菌测序时被"抽中"的概率下降，因此定量结果中背景菌占比会下降。因此，将某变量的丰度结果（y）和样本浓度（x）构建的线性模型。如果该模型较零模型（
）具有显著差异（F检验，阈值默认0.1），则判定该变量为背景菌。
如何统计检验
算法基于F分布。F分布定义如下：
,则
假设有一组数据如下：
load('decontam_data.rdata')
head(freq_data)
logc logf
1 7.947679 -1.9106629
2 8.666992 -0.6715281
3 8.527539 -0.8704869
4 7.112327 -1.6684781
5 8.261526 -0.6038552
6 8.677951 -0.9100405
logc和logf分别为样本的荧光定量结果对数化和定量丰度结果对数化。
frequency构建的两个个模型如下：
# model 1
lm1 <- lm(logf~offset(-1*logc), data=freq_data) # y=-x+b
lm1
Call:
lm(formula = logf ~ offset(-1 * logc), data = freq_data)
Coefficients:
(Intercept)
6.679
模型1为。
# model 2
lm0 <- lm(logf~1, data=freq_data) # y=mean x
lm0
Call:
lm(formula = logf ~ 1, data = freq_data)
Coefficients:
(Intercept)
-1.315
零模型为。
数据和模型绘制散点图如下：
freq_data %>%
ggplot(aes(logc, logf))+
geom_point()+
geom_smooth(formula = 'y~offset(-1*x)', method = 'lm', color='red',linetype=1,se =F)+
geom_smooth(formula = 'y~1', method = 'lm', color='black',linetype=2, se=F)
红色线条为拟定的拟合模型，黑色虚线为零模型，即作者教程中展示的图片：
分别表示拟合模型和零模型的残差平方和，服从
分布，进行
检验（左侧检验，
）。
计算统计量：
dof <- nrow(freq_data)-1 # 548
SS1 <- sum(lm1$residuals^2)
SS0 <- sum(lm0$residuals^2)
F <- SS1/SS0
F
[1] 3.570293
计算上分位数：
alpha <- 0.1
q <- qf(alpha, dof, dof) # 左侧检验
q
[1] 0.8962215
自由度的
分布图如下：
x <- seq(0, 4,length = 1000)
y <- df(x, dof, dof)
plot(x, y, type="l", main=sprintf('F(%s,%s)分布',dof,dof),,xlab='', ylab='density')
abline(h=0, v=c(F,q), col="gray",lty=3:2)
x_sub <- x[x < q]
y_sub <- y[x < q]
polygon(rbind(c(q,0),cbind(x_sub,y_sub)),border=NA,col="gray")
arrows(c(0.5,F+0.2), c(1,1), c(0.85,F), c(0.7,0.7))
text(c(0.5,F+0.2), c(1.1,1.2), labels=c('alpha','F value'))
左侧检验的拒绝域位于左侧小尾巴。统计量
F大于上计算p值：
p <- pf(F, dof, dof)
p
[1] 1
因为，无法拒绝
，因此判定为病原菌。
Prevalence
frequency是根据浓度梯度构建丰度与荧光强度之间的线性关系，进而判定线性模型是否具有显著性。
更常见的做法是通过阴性对照样本判定该变量是否为背景菌。主要思路如作者所说，统计某个病原在测试组和阴性对照中中出现的比率，因此问题转换成两个样本比例的统计检验或fisher检验。
如何统计检验
定量数据转换成丰度，同时包含组别信息：
head(preva_data)
Seq1 neg
1 0.1479823 FALSE
2 0.5109272 FALSE
3 0.4187476 FALSE
4 0.1885338 FALSE
5 0.5466999 FALSE
6 0.4025079 FALSE
Seq1为待检验的病原，neg指定组别，False表示为测试样本，True表示对照样本。
统计指标为该病原是否在样本中检出，因此构建2x2列联表如下：
freq=factor(preva_data$Seq1>0, levels=c(TRUE, FALSE))
neg=factor(preva_data$neg, levels=c(TRUE, FALSE))
tab <- table(freq, neg)
tab
neg
freq TRUE FALSE
TRUE 12 537
FALSE 18 2
上述结果中，行表示是否在样本中检出，列指定是否为对照样本。阴性对照组中共30个样本，12个样本检出有该病原，因此比率为12/30。同理，测试组中99%(537/539)的样本检出该病原。
作者进行右侧检验，两个假设为（列联表中，行success表示病原，列success表示对照组，p表示病原菌的概率）：
p=prop.test(tab, alternative="greater")$p.value
p
[1] 1
右侧检验拒绝域位于右侧小尾巴。位于接受域，因此接受
，即病原为病原菌。
小结
frequency算法统计连续性数值变量，prevalence统计分类变量。基于F分布进行统计检验。
补充
测试数据
- 本文测试数据地址。
部分源码
Frequency:
lm1 <- lm(logf~offset(-1*logc), data=df)
SS1 <- sum(lm1$residuals^2)
lm0 <- lm(logf~1, data=df) # fit an intersept only
SS0 <- sum(lm0$residuals^2)
dof <- sum(freq>0)-1
pval <- pf(SS1/SS0,dof,dof)
Prevalence:
fisher.pval <- function(tab, alternative) {
excess <- fisher.test(tab, alternative="greater")$p.value +
fisher.test(tab, alternative="less")$p.value - 1
pval <- fisher.test(tab, alternative=alternative)$p.value
pval <- pval - excess/2
pval
}
pval <- tryCatch(
prop.test(tab, alternative="greater")$p.value,
warning=function(w) fisher.pval(tab, alternative="greater")
)
pval
