“清晰而又吸引人——这无疑是学习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。
