网上的答案经常不靠谱，包括我的

通常情况下我会使用 featureCounts 得到表达矩阵是 raw counts, 但总是有人需要我转换成各种形式，比如 RPKM, FPKM and TPM

概念见：https://statquest.org/2015/07/09/rpkm-fpkm-and-tpm-clearly-explained/

起初，想偷一下懒，就搜索到了 这个答案：http://ny-shao.name/2016/11/18/a-short-script-to-calculate-rpkm-and-tpm-from-featurecounts-output.html

首先网络代码是错的

朋友提醒我，其实错了，因为很容易校验，明显 colSums(exprSet_tpm) 并不是一百万。（奇怪的是我居然没有对我的第一次代码进行同样的校验）

我第一次修改的代码仍然是错的

其实我老早就写过TPM公式，就是RPKM的百分比的百万倍扩大值，所以还是自己动手重新写了代码。

rpkm <- function(counts, lengths) {
  rate <- counts / lengths
  rate / sum(counts) * 1e9
}

exprSet_rpkm=rpkm(exprSet,len) 
exprSet_tpm=1e6*exprSet_rpkm/colSums(exprSet_rpkm)


不过，比较奇怪的是这个时候 colSums(exprSet_tpm) 是接近一百万，而不是精确的一百万，我当时还没有想清楚具体缘由，是不是R的计算小数点问题。

有关于的讨论：http://blog.nextgenetics.net/?e=51#body-anchor

有趣的是， 因为自己并不使用这个RPKM值，所以后面也没有继续校验代码，知道昨天学徒在使用我的数据时候很认真的发现了这个bugs并且指出来了。

不知道第一次发布这个教程，有多少人看了，如果真的有需求，理论上需要严格检查我的代码。

第二次修改

这次代码结合了我在单细胞课程的代码，方法一：

exprSet[1:4,1:4] 
len=meta[match(rownames(exprSet),meta$Geneid),'Length']
head(len)
rpkm <- function(counts, lengths) {
  # 首先对矩阵进行基因长度归一化
  # 矩阵除以向量是按照行分开，表达矩阵的行是基因，所以每个基因除以各自的基因长度
  rate <- counts / lengths


  # 然后对矩阵进行文库大小归一化, 就复杂一点    # 注意这里的两次转置。


  t(t(rate) / colSums(counts) )* 1e9
}

exprSet_rpkm=rpkm(exprSet,len)  


方法二：

lengths=len
head(lengths)
head(rownames(exprSet))
# http://www.biotrainee.com/thread-1791-1-1.html
exprSet[1:4,1:4]
total_count<- colSums(exprSet)
head(total_count)
head(lengths)
total_count[4]
lengths[1]
1*10^9/(lengths[1]*total_count[4])
rpkm <- t(do.call( rbind,
                   lapply(1:length(total_count),
                          function(i){
                            10^9*exprSet[,i]/lengths/total_count[i]
                          }) ))

rpkm[1:4,1:4]
exprSet_rpkm[1:4,1:4]


嗯，差不多可以肯定，最后定稿的代码是对的了，可以看到两个公式效果一样，而且我还比较了其中一个基因：

1*10^9/(lengths[1]*total_count[4])