你好,

根据最近发表的文章(https://genomebiology.biomedcentral.com/articles/10.1186/s13059-020-02104-1)我一直在尝试将TMM规范化(使用edger::calcNormFactors)应用到我使用过的许多不同的微生物组数据集。在使用它们时，我注意到，当我打乱原始读取计数矩阵中的样本(列)的顺序时，一些数据集的规范化因子无法再现。我用9个不同的数据集测试了这种行为，发现当矩阵> 0.79稀疏度时通常会出现这种情况。我想弄清楚为什么会这样以及为什么这会影响标准化因子的重现性。我已经包括一个dropbox的链接，包含感兴趣的数据集以及rmarkdown文件与所有的代码。非常感谢您的指导。

https://www.dropbox.com/sh/6fo3fipyc0p36f3/AACnNhDmOk5U285g_u9g0fTba?dl=0

这里有几行代码可以重现这个问题。

下载TSV文件https://www.dropbox.com/s/w6l11rfyh8z19wl/BISCUIT_ASVs_table.tsv?dl=1count_matrix.tsv <- as.matrix(read.table("count_matrix. tsv")Tsv "， sep="\t"， row.names = 1, comment。char =""， skip=1, header=T, check.names =F, quote="")) count_norm1 <- calcNormFactors(饼干计数，方法="TMM") #shuffle列order shuffled_df <- shuffled_count [，sample(c(1:38)， 38, replace=F)] shuffled_norm <- calcNormFactors(shuffled_df，sorted_norm1 <- count_norm1[sort(names(count_norm1))] sorted_shuffle_norm <- shuffled_norm[sort(names(shuffled_norm))] identical(sorted_norm1, sorted_shuffle_norm) corr .test(sorted_norm1, sorted_shuffle_norm))

谢谢,雅各接近

calcNormFactors应该返回相同的归一化因子，不管样本的顺序。

我很欣赏您提供的代码和数据，这很好，但我没有时间来详细分析您提供的复杂分析会话。您应该能够用一个计数矩阵和几行代码来演示这个问题。你能做到吗?代码行可以复制到你的问题中——我们不需要html或标记。

20分钟后

好的，我看到问题了:

> x < - read.delim(“ArcticFireSoils_ASVs_table.tsv”,跳过= 2,头= FALSE, row.names = 1) >暗(x) [1] 37244 148 > set.seed j(20200912) > < -样本(148)> f1 < - calcNormFactors (x) > f2 < - calcNormFactors (x [j]) >总结(f1 [j] - f2)分钟。1日瞿瞿。中值的意思是3日,马克斯。-0.55507 -0.13682 -0.01441 0.01381 0.11001 1.83245

我来看看。

又过了20分钟

这个问题是由上四分位归一化方法引起的。TMM方法要求选择其中一个样本作为参考样本，与其他所有样本进行比较。默认情况下，使用上四分位数方法来选择参考样本。的calcNormFactors帮助页面说:

如果refColumn未指定，则使用其上四分位数最接近平均上四分位数的库。

但是对于您的数据，超过75%的数据值为零，因此上四分位数为零，因此上四分位数的归一化因子是未定义的:

> summary(calcNormFactors(x,method="upperquartile")) Max. Max. Max. Max. Max. Max. Max. Max. Max. Max. Max. Max。NA的NA NA NA NA NA NA NA NA NA NA 148

因此，TMM只好选择第一列作为参考样本。如果您打乱列，那么第一列将是不同的，从而导致不同的参考样本和不同的归一化因素。

您可以通过以一致和健壮的方式选择参考样本来解决这个问题。我还建议改用TMMwsp代替普通的TMM，因为它对稀疏数据更可靠:

> Ref1 < - which.max (colSums (sqrt (x))) > f1 < - calcNormFactors (x,方法=“TMMwsp refCol = Ref1) > Ref2 < - which.max (colSums (sqrt (x [j]))) > f2 < - calcNormFactors (x [j],方法=“TMMwsp refCol = Ref2) >总结(f1 [j] - f2)分钟。1日瞿瞿。中值的意思是3日,马克斯。0 0 0 0 0 0