在做表达矩阵的counts值作为RPKM的时候发现的这个知识点细节问题， 因为矩阵需要每一个样本除以它各自的文库大小，然后呢，每个基因又需要除以各自的基因长度。
所以呢，我们的表达矩阵，其实是需要除以两个长度不一的向量，而且方向不一样，一个是按照行来除以，一个是按照列来除以，我最后写的代码是： rpkm <- function(counts, lengths) {
# 首先对矩阵进行基因长度归一化
# 矩阵除以向量是按照行分开，表达矩阵的行是基因，所以每个基因除以各自的基因长度
rate <- counts / lengths
# 然后对矩阵进行文库大小归一化
t(t(rate) / colSums(counts) )* 1e9
}
对很多朋友来说，是需要解释一下的。 很明显 counts 是表达矩阵，lengths 是不同基因长度向量，而 colSums(counts) 是不同样本的长度向量。 一个简单的例子这里还是生成随机数： counts=1:10
dim(counts)=c(2,5)
lengths=c(1:2)
lib=1:5
counts/lengths
counts/lib
t(t(counts)/lib)
结果如下： 
可以看到，矩阵除以向量，是按行的顺序来的，如果需要列，就得先转置，再转回来。
|
