用DESeq2对表达矩阵进行归一化的详细步骤

“清晰而又吸引人——这无疑是学习R的有趣方式！”
                                                               ——Amos A. Folarin，伦敦大学学院

1.对readscount进行ln()转换

表达矩阵的行代表feature（如基因、外显子等），列代表sample；对于原始的每一个基因每一个样本对应的reads数进行以e为底的转换：ln(counts)ln()函数在R语言中是log()

log(counts)

2.对每一行的值进行均值计算

对该基因对应的所有样本的转换后的readscount计算均值

rowMeans(log(counts))

3.ln(counts)减去每一行的均值

log(counts)-rowMeans(log(counts))

4.筛选出step3得到的矩阵中不包含-Inf，NaN的基因行

(log(counts)-rowMeans(log(counts)))[is.finite(rowMeans(log(counts))),]

5.计算step4得到的矩阵中每列的中位数

median((log(counts)-rowMeans(log(counts)))[is.finite(rowMeans(log(counts))),])

6.计算每列的e^中位数，该数即为sizefactor

exp(median((log(counts)-rowMeans(log(counts)))[is.finite(rowMeans(log(counts))),]))

7.原始每一个readscount除以对应列的sizefactor即为归一化后的基因表达量

original readscount/size factors即为最终的每个样本每个基因对应的归一化的表达量。每列中的任意original readscount除以该列的sizefactor。

以上代码只具有表面含义，因为计算时各种数据类型的存在，实际在R中运算较为复杂。当然，以上步骤只是用DESeq2去estimating size factors的内部算法，实际可以直接利用sizeFactors(dds)去计算size factors。