前言
實(shí)驗(yàn)過程中，背景細(xì)菌的存在可能導(dǎo)致交叉污染，引起假陽性的發(fā)生，從而給mNGS結(jié)果解讀帶來困擾。mNGS專家共識(shí)中指出，實(shí)驗(yàn)室應(yīng)該建立背景細(xì)菌庫并定時(shí)更新，用于減少假陽性。
decontam通過濃度梯度（frequency）或陰性對(duì)照的方式（Prevalence）鑒定實(shí)驗(yàn)過程中的背景菌。因此本文旨在介紹decontam內(nèi)部細(xì)節(jié)。
正文開始前，有如下注意事項(xiàng)：
本文主要集中于decontam原理解讀，因此使用方法和結(jié)果解讀不過多留戀，讀者請(qǐng)參考官方教程。
下文只涉及官網(wǎng)教程中提到的兩個(gè)算法（frequency和prevalence）的核心部分，細(xì)枝末節(jié)和其他算法略過不提。
案例數(shù)據(jù)和代碼均摘自作者源碼，同時(shí)結(jié)合自身理解。本人非數(shù)學(xué)科班出身，若理解有誤，歡迎大家批評(píng)指正。
基本使用
suppressPackageStartupMessages({
library(rio)
library(ggplot2)
library(dplyr)
library(decontam)
})
使用官方數(shù)據(jù)：
# 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
結(jié)果最后一列即判定該變量或病原是否為背景菌。
詳情見官方教程.
兩個(gè)算法
frequency
思路如下：
隨著樣本濃度的增加，文庫中的背景菌測(cè)序時(shí)被"抽中"的概率下降，因此定量結(jié)果中背景菌占比會(huì)下降。因此，將某變量的豐度結(jié)果（y）和樣本濃度（x）構(gòu)建的線性模型。如果該模型較零模型（
）具有顯著差異（F檢驗(yàn)，閾值默認(rèn)0.1），則判定該變量為背景菌。
如何統(tǒng)計(jì)檢驗(yàn)
算法基于F分布。F分布定義如下：
,則
假設(shè)有一組數(shù)據(jù)如下：
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分別為樣本的熒光定量結(jié)果對(duì)數(shù)化和定量豐度結(jié)果對(duì)數(shù)化。
frequency構(gòu)建的兩個(gè)個(gè)模型如下：
# 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
零模型為。
數(shù)據(jù)和模型繪制散點(diǎn)圖如下：
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)
紅色線條為擬定的擬合模型，黑色虛線為零模型，即作者教程中展示的圖片：
分別表示擬合模型和零模型的殘差平方和，服從
分布，進(jìn)行
檢驗(yàn)（左側(cè)檢驗(yàn)，
）。
計(jì)算統(tǒng)計(jì)量：
dof <- nrow(freq_data)-1 # 548
SS1 <- sum(lm1$residuals^2)
SS0 <- sum(lm0$residuals^2)
F <- SS1/SS0
F
[1] 3.570293
計(jì)算上分位數(shù)：
alpha <- 0.1
q <- qf(alpha, dof, dof) # 左側(cè)檢驗(yàn)
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'))
左側(cè)檢驗(yàn)的拒絕域位于左側(cè)小尾巴。統(tǒng)計(jì)量
F大于上計(jì)算p值：
p <- pf(F, dof, dof)
p
[1] 1
因?yàn)?img class="math-inline" src="https://math.jianshu.com/math?formula=p%3E0.1" alt="p>0.1" mathimg="1">，無法拒絕，因此判定為病原菌。
Prevalence
frequency是根據(jù)濃度梯度構(gòu)建豐度與熒光強(qiáng)度之間的線性關(guān)系，進(jìn)而判定線性模型是否具有顯著性。
更常見的做法是通過陰性對(duì)照樣本判定該變量是否為背景菌。主要思路如作者所說，統(tǒng)計(jì)某個(gè)病原在測(cè)試組和陰性對(duì)照中中出現(xiàn)的比率，因此問題轉(zhuǎn)換成兩個(gè)樣本比例的統(tǒng)計(jì)檢驗(yàn)或fisher檢驗(yàn)。
如何統(tǒng)計(jì)檢驗(yàn)
定量數(shù)據(jù)轉(zhuǎn)換成豐度，同時(shí)包含組別信息：
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為待檢驗(yàn)的病原，neg指定組別，F(xiàn)alse表示為測(cè)試樣本，True表示對(duì)照樣本。
統(tǒng)計(jì)指標(biāo)為該病原是否在樣本中檢出，因此構(gòu)建2x2列聯(lián)表如下：
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
上述結(jié)果中，行表示是否在樣本中檢出，列指定是否為對(duì)照樣本。陰性對(duì)照組中共30個(gè)樣本，12個(gè)樣本檢出有該病原，因此比率為12/30。同理，測(cè)試組中99%(537/539)的樣本檢出該病原。
作者進(jìn)行右側(cè)檢驗(yàn)，兩個(gè)假設(shè)為（列聯(lián)表中，行success表示病原，列success表示對(duì)照組，p表示病原菌的概率）：
p=prop.test(tab, alternative="greater")$p.value
p
[1] 1
右側(cè)檢驗(yàn)拒絕域位于右側(cè)小尾巴。位于接受域，因此接受
，即病原為病原菌。
小結(jié)
frequency算法統(tǒng)計(jì)連續(xù)性數(shù)值變量，prevalence統(tǒng)計(jì)分類變量?；贔分布進(jìn)行統(tǒng)計(jì)檢驗(yàn)。
補(bǔ)充
測(cè)試數(shù)據(jù)
部分源碼
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
