参考文章: 基因长度之多少
示例数据下载：rawCounts & rpkmNormalized
RPKM及CPM的定义
顺便说一下，RPKM和FPKM早就被诟病了，它们并不是一个有意义的统计量，本质上是不可解释的，具体可参考为什么说FPKM和RPKM都错了？。然而，无论RPKM或FPKM多么的遭人诟病，它的真实需求还是存在， 其中关键的部分就是算
geneLength，我们该如何计算呢？
这个时候需要理解：基因，转录本(transcripts,isoform，mRNA序列)、EXON区域，cDNA序列、UTR区域，ORF序列、CDS序列 这些概念，一个基因可以转录为多个转录本，真核生物里面每个转录本通常是由一个或者多个EXON组成，能翻译为蛋白的EXON区域是CDS区域，不能翻译的那些EXON的开头和结尾是UTR区域，翻译区域合起来是ORF序列，而转录本逆转录就是cDNA序列
有了这些概念，我们就能理解目前主流定义基因长度的几种方式：
- 挑选基因的最长转录本
- 选取多个转录本长度的平均值
- 非冗余外显子(EXON)长度之和
- 非冗余 CDS（Coding DNA Sequence） 长度之和
非冗余exon长度
我们先来看基因长度为 非冗余exon长度之和的计算方式：
library("TxDb.Mmusculus.UCSC.mm10.knownGene")
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
## 下面是定义基因长度为 非冗余exon长度之和
exon_txdb <- exons(txdb) # 得到外显子
genes_txdb <- genes(txdb) # 得到基因
o <- findOverlaps(exon_txdb, genes_txdb) # 寻找重叠部分，即基因外显子对应情况
t1 <- exon_txdb[queryHits(o)] # 筛选外显子
t2 <- genes_txdb[subjectHits(o)] # 筛选基因
t1 <- as.data.frame(t1)
t1$geneid <- mcols(t2)[ ,1] # 得到基因外显子对应情况，同时包含外显子的长度范围
# 如果觉得速度不够，就参考R语言实现并行计算
# http://www.bio-info-trainee.com/956.html
# 函数split()可以按照分组因子，把向量，矩阵和数据框进行适当的分组
g_l <- lapply(split(t1, t1$geneid), function(x){
tmp <- apply(x, 1, function(y){ # tmp是一个列表，
y[2]:y[3] # 每个元素为geneid对应的其中一个外显子的长度向量
})
length(unique(unlist(tmp))) # 剔除外显子中重复的部分，剩余的长度即为基因的长度
})
head(g_l) # lapply函数返回的值是列表形式 ,apply则为向量、矩阵或列表
g_l <- data.frame(gene_id=names(g_l), length=as.numeric(g_l)) # 封装成数据框，便于操作
save(g_l, file = 'g_l.Rdata')
# 保存数据，以后通过load(file = 'g_l.Rdata')读取
随后，我们通过将geneid换成symbol，筛选过后即可算RPKM值，具体如下：
library(org.Mm.eg.db) # 根据需求来选，如org.Mm.eg.db等
s2g <- toTable(org.Mm.egSYMBOL)
head(s2g)
g_l <- merge(g_l, s2g, by='gene_id') #把g_l,s2g两个数据框以'gene_id'为连接进行拼接
ng <- intersect(rownames(a), g_l$symbol) #取a数据框的行名与g_l数据框的symbol列的交集
a <- read.table('GSE111229_Mammary_Tumor_fibroblasts_768samples_rawCounts.txt.gz',
header = T , sep = '\t') # 读入原始表达矩阵
exprSet <- a[ng, ] # 筛选出表达矩阵
lengths <- g_l[match(ng, g_l$symbol), 2] # 以ng的顺序取基因长度
# http://www.biotrainee.com/thread-1791-1-1.html
total_count <- colSums(exprSet) # 对exprSet求基因文库的大小
rpkm <- t(do.call( rbind, # 将下列的结果按行合并
lapply(1:length(total_count), # 对每列细胞做同样操作
function(i){
10^9*exprSet[ , i]/lengths/total_count[i] # RPKM的定义
}) ))
rpkm[1:4, 1:4] # 注意，此时并没有行名和列名，后续可自行添加
# 下面可以比较一下 自己根据counts值算出来的RPKM和作者提供的RPKM区别。
b <- read.table('GSE111229_Mammary_Tumor_fibroblasts_768samples_rpkmNormalized.txt.gz',
header = T , sep = '\t')
# 检测数据，观察差异
b[1:4, 1:4]
rpkm_paper <- b[ng, ]
rpkm_paper[1:4, 1:4]
rpkm[1:4, 1:4]
最长转录本长度
不同于非冗余exon长度那样计算麻烦，有专门的函数可以得到最长转录本长度，代码如下所示：
library("TxDb.Mmusculus.UCSC.mm10.knownGene")
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
t_l <- transcriptLengths(txdb) # 基因id与转录本的对应情况，有转录本长度的信息
t_l <- na.omit(t_l) # 忽略其中的缺失值
t_l <- t_l[order(t_l$gene_id,t_l$tx_len,decreasing = T), ] # 将t_l按gene_id和tx_len递减排序
t_l <- t_l[!duplicated(t_l$gene_id), ] # 保留最长的转录本（一个geneid对应多个转录本）
head(t_l)
g_l <- t_l[ , c(3,5)] # 取第3、5列，与非冗余exon长度的g_l一致
后续操作请参考非冗余exon长度的第二部分代码
多个转录本长度的平均值
可参考最长转录本长度，只不过我们不取最大值，而是取平均值
非冗余CDS之和
CDS暂时还未考虑，后续会添加的
