參考文章: 基因長度之多少
示例數(shù)據(jù)下載：rawCounts & rpkmNormalized
RPKM及CPM的定義
順便說一下遣疯，RPKM和FPKM早就被詬病了唱蒸，它們并不是一個有意義的統(tǒng)計量闯团，本質(zhì)上是不可解釋的篮绿，具體可參考為什么說FPKM和RPKM都錯了？。然而，無論RPKM或FPKM多么的遭人詬病祸挪，它的真實需求還是存在， 其中關(guān)鍵的部分就是算
geneLength贞间，我們該如何計算呢贿条？
這個時候需要理解：基因，轉(zhuǎn)錄本(transcripts,isoform增热，mRNA序列)整以、EXON區(qū)域，cDNA序列峻仇、UTR區(qū)域公黑，ORF序列、CDS序列 這些概念摄咆，一個基因可以轉(zhuǎn)錄為多個轉(zhuǎn)錄本凡蚜，真核生物里面每個轉(zhuǎn)錄本通常是由一個或者多個EXON組成，能翻譯為蛋白的EXON區(qū)域是CDS區(qū)域吭从，不能翻譯的那些EXON的開頭和結(jié)尾是UTR區(qū)域朝蜘，翻譯區(qū)域合起來是ORF序列，而轉(zhuǎn)錄本逆轉(zhuǎn)錄就是cDNA序列
有了這些概念涩金，我們就能理解目前主流定義基因長度的幾種方式：
- 挑選基因的最長轉(zhuǎn)錄本
- 選取多個轉(zhuǎn)錄本長度的平均值
- 非冗余外顯子(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語言實現(xiàn)并行計算
# http://www.bio-info-trainee.com/956.html
# 函數(shù)split()可以按照分組因子副渴，把向量，矩陣和數(shù)據(jù)框進行適當?shù)姆纸M
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函數(shù)返回的值是列表形式 ,apply則為向量、矩陣或列表
g_l <- data.frame(gene_id=names(g_l), length=as.numeric(g_l)) # 封裝成數(shù)據(jù)框讼载，便于操作
save(g_l, file = 'g_l.Rdata')
# 保存數(shù)據(jù)轿秧，以后通過load(file = 'g_l.Rdata')讀取
隨后，我們通過將geneid換成symbol咨堤，篩選過后即可算RPKM值菇篡，具體如下：
library(org.Mm.eg.db) # 根據(jù)需求來選，如org.Mm.eg.db等
s2g <- toTable(org.Mm.egSYMBOL)
head(s2g)
g_l <- merge(g_l, s2g, by='gene_id') #把g_l,s2g兩個數(shù)據(jù)框以'gene_id'為連接進行拼接
ng <- intersect(rownames(a), g_l$symbol) #取a數(shù)據(jù)框的行名與g_l數(shù)據(jù)框的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, # 將下列的結(jié)果按行合并
lapply(1:length(total_count), # 對每列細胞做同樣操作
function(i){
10^9*exprSet[ , i]/lengths/total_count[i] # RPKM的定義
}) ))
rpkm[1:4, 1:4] # 注意一喘，此時并沒有行名和列名驱还，后續(xù)可自行添加
# 下面可以比較一下 自己根據(jù)counts值算出來的RPKM和作者提供的RPKM區(qū)別。
b <- read.table('GSE111229_Mammary_Tumor_fibroblasts_768samples_rpkmNormalized.txt.gz',
header = T , sep = '\t')
# 檢測數(shù)據(jù)凸克，觀察差異
b[1:4, 1:4]
rpkm_paper <- b[ng, ]
rpkm_paper[1:4, 1:4]
rpkm[1:4, 1:4]
最長轉(zhuǎn)錄本長度
不同于非冗余exon長度那樣計算麻煩议蟆，有專門的函數(shù)可以得到最長轉(zhuǎn)錄本長度，代碼如下所示：
library("TxDb.Mmusculus.UCSC.mm10.knownGene")
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
t_l <- transcriptLengths(txdb) # 基因id與轉(zhuǎn)錄本的對應情況萎战，有轉(zhuǎn)錄本長度的信息
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), ] # 保留最長的轉(zhuǎn)錄本（一個geneid對應多個轉(zhuǎn)錄本）
head(t_l)
g_l <- t_l[ , c(3,5)] # 取第3咐容、5列，與非冗余exon長度的g_l一致
后續(xù)操作請參考非冗余exon長度的第二部分代碼
多個轉(zhuǎn)錄本長度的平均值
可參考最長轉(zhuǎn)錄本長度蚂维，只不過我們不取最大值戳粒，而是取平均值
非冗余CDS之和
CDS暫時還未考慮，后續(xù)會添加的
