劉小澤寫于19.7.9+12-第二單元第十講：RPKM概念及計(jì)算方法
筆記目的：根據(jù)生信技能樹的單細(xì)胞轉(zhuǎn)錄組課程探索smart-seq2技術(shù)相關(guān)的分析技術(shù)
課程鏈接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53

RPKM須知

核心就在于基因長(zhǎng)度的計(jì)算

參考這個(gè)網(wǎng)站，做的還是很清爽的：http://www.metagenomics.wiki/pdf/definition/rpkm-calculation

它的公式看上去很簡(jiǎn)單，就是對(duì)基因長(zhǎng)度（也就是公式里的K）以及文庫(kù)大小(也就是M) 進(jìn)行標(biāo)準(zhǔn)化，其中文庫(kù)大小很好理解，就是一個(gè)樣本全部測(cè)序reads之和。但是這里有一個(gè)之前也沒有注意到的問題，基因長(zhǎng)度是什么？怎么計(jì)算出來的？

幾個(gè)問題：它是簡(jiǎn)單的終止坐標(biāo)減去起始坐標(biāo)嗎？它是外顯子長(zhǎng)度之和嗎？那外顯子之間有重合怎么處理？它是某一條轉(zhuǎn)錄本的長(zhǎng)度嗎？還是說基因長(zhǎng)度是一個(gè)轉(zhuǎn)錄本的平均值？

這個(gè)問題是要好好總結(jié)一下，但首先還是要明確一些基因相關(guān)名詞：

• 基因與轉(zhuǎn)錄本：一個(gè)基因可以轉(zhuǎn)錄成多個(gè)轉(zhuǎn)錄本
• 轉(zhuǎn)錄本與外顯子：真核生物的每個(gè)轉(zhuǎn)錄本一般是由一個(gè)或多個(gè)外顯子組成
• 轉(zhuǎn)錄本與cDNA：轉(zhuǎn)錄本逆轉(zhuǎn)錄得到cDNA
• 外顯子與編碼區(qū)CDS、非編碼區(qū)UTR：可以翻譯成蛋白的外顯子區(qū)域是CDS區(qū)域，不能翻譯的外顯子開頭、結(jié)尾是UTR區(qū)域
• CDS與開放閱讀框ORF：ORF是從起始密碼子(ATG, but not always)到終止密碼子(TAA, TAG, TGA)的DNA序列，另外又有兩種鏈的方向，因此總共有6種閱讀框，它也會(huì)包含內(nèi)含子(這導(dǎo)致了真核生物的CDS與ORF不一致；另外在原核生物中它們是一樣的)
  http://www.scfbio-iitd.res.in/tutorial/orf.html

中心法則方向："DNA makes RNA makes Protein"

• 外顯子和內(nèi)含子是DNA層面的gene feature，密碼子(codon)是RNA層面上的feature；外顯子和內(nèi)含子都存在于雙鏈DNA（dsDNA）的蛋白編碼基因區(qū)域中，一般通過sense strand，也即是5’=》3‘查看

• 然后dsDNA變成hnRNA (不均一核RNA：heterogeneous nuclear RNA)，和dsDNA的 5'-3' 方向一致，另外將T堿基換成了U。這里雖然hnRNA中也含有dsDNA的exon、intron序列，但稱呼要換："exon transcripts"和"intron transcripts"

• 接著hnRNA進(jìn)行了excision ('splicing out') 操作，切掉了intron transcripts，將剩余的exon transcripts拼接起來，得到mRNA。mRNA的序列是和DNA的exon序列一致的，它也可以視作由dsDNA的sense strand序列將T變U。密碼子在DNA中叫"condons"，在mRNA中叫"triplets"

  In bioinformatics, the 64 triplets are sometimes presented as a "translation table"** that can be used directly with the Sense Strand sequence to infer the protein sequence.

https://www.mun.ca/biology/scarr/Exons_Introns_Codons.html

關(guān)于基因長(zhǎng)度

既然基因包含這么多"組件"，那么求它的長(zhǎng)度也會(huì)有幾種方法：

• 選擇最長(zhǎng)的轉(zhuǎn)錄本
• 多個(gè)轉(zhuǎn)錄本的均值
• 非冗余外顯子長(zhǎng)度之和
• 非冗余CDS之和

注意到這里的"非冗余"，就是存在一個(gè)基因的多個(gè)外顯子之間存在重疊（比如基因A的1號(hào)外顯子較短，2號(hào)外顯子長(zhǎng)，1號(hào)包含在2號(hào)中），單純的相加會(huì)重復(fù)計(jì)算

計(jì)算基因長(zhǎng)度

第一種：非冗余外顯子之和

載入原始表達(dá)矩陣

rm(list = ls()) 
options(stringsAsFactors = F)
load(file = '../input.Rdata')
a[1:4,1:4]
head(df)

選擇哪個(gè)包來計(jì)算基因長(zhǎng)度呢？

不知道就搜索一下：

然后找到TxDb這個(gè)包的介紹，看Bioconductor這本電子書的41頁(yè) 3.4.19節(jié)：https://bioconductor.github.io/BiocWorkshops/BioC2018.pdf

起步

因?yàn)樘幚淼氖切∈髷?shù)據(jù)，所以要加載相應(yīng)小鼠的包：

library("TxDb.Mmusculus.UCSC.mm10.knownGene")
# 加載包以后，看看幫助文檔，發(fā)現(xiàn)這么一句話(使用它只需要call它的名字)
## show the db object that is loaded by calling it's name

txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
# 依然不知道怎么用，那么就再看幫助文檔的"see also"信息

點(diǎn)進(jìn)去發(fā)現(xiàn)，存在好多函數(shù)，因?yàn)槲覀兿氩捎梅侨哂嗤怙@子加和的方法計(jì)算基因長(zhǎng)度，因此選擇exon和genes就好了

# 取出exon和gene
exon_txdb=exons(txdb)
genes_txdb=genes(txdb)
>   genes_txdb
GRanges object with 24368 ranges and 1 metadata column:
            seqnames              ranges strand |     gene_id
               <Rle>           <IRanges>  <Rle> | <character>
  100009600     chr9   21062393-21075496      - |   100009600
  100009609     chr7   84940169-84964009      - |   100009609
# 看到結(jié)果是一個(gè)GRanges object，那么不懂就再搜索什么是GRanges
?GRanges
# The GRanges class is a container for the genomic locations and their associated annotations. 就是存儲(chǔ)基因組坐標(biāo)和相關(guān)注釋信息的"容器"

因?yàn)槿〕龅膃xon和gene都是坐標(biāo)信息(可以看到上面的21062393-21075496就是坐標(biāo)區(qū)間)，那么二者既然都是坐標(biāo)而且不一致，就會(huì)有交叉，如何找到這個(gè)交叉？

又引入了一個(gè)新的函數(shù)=》findOverlaps

o = findOverlaps(exon_txdb,genes_txdb)
>   o
Hits object with 252602 hits and 0 metadata columns:
           queryHits subjectHits
           <integer>   <integer>
       [1]         1        6729
       [2]         2        6729
       [3]         3        6729
       [4]         4        6729
       [5]         5        6729
       ...       ...         ...
# 其中exon_txdb因?yàn)閷懺谇懊?，于是它就作為queryHits；可以看到：它的第一個(gè)元素和genes_txdb的6729存在交叉，取出來看一下結(jié)果
> genes_txdb[6729]
GRanges object with 1 range and 1 metadata column:
        seqnames          ranges strand |     gene_id
           <Rle>       <IRanges>  <Rle> | <character>
  18777     chr1 4807893-4846735      + |       18777
> exon_txdb[1]
GRanges object with 1 range and 1 metadata column:
      seqnames          ranges strand |   exon_id
         <Rle>       <IRanges>  <Rle> | <integer>
  [1]     chr1 4807893-4807982      + |         1
# exon的第一個(gè)元素是chr1 4807893-4807982，對(duì)應(yīng)上了基因的第6729個(gè)元素：chr1 4807893-4846735，于是這個(gè)exon就屬于這個(gè)基因

既然得到了所有的overlap，那就分別提取出來exon和gene的信息

t1=exon_txdb[queryHits(o)]
t2=genes_txdb[subjectHits(o)]

t1=as.data.frame(t1) #更直觀地來查看t1，結(jié)果將25w個(gè)外顯子的坐標(biāo)

為了給t1一個(gè)對(duì)應(yīng)的基因ID，需要使用t2

t1$geneid=mcols(t2)[,1]
# mcols作用是提取 a DataFrame object containing the metadata columns
# 看下面??這個(gè)，mcols(t2)[,1]提取的也就是虛線右側(cè)的一列"gene_id"
> genes_txdb[6729]
GRanges object with 1 range and 1 metadata column:
        seqnames          ranges strand |     gene_id
           <Rle>       <IRanges>  <Rle> | <character>
  18777     chr1 4807893-4846735      + |       18777

再看一眼結(jié)果：

關(guān)鍵一步

好，看到基因18777對(duì)應(yīng)著10個(gè)exon，那么如何求這個(gè)18777的長(zhǎng)度呢？是簡(jiǎn)單把10個(gè)exon的第四列加起來嗎？

不是的！請(qǐng)看：第8和第9個(gè)exon，它們是不是有重疊？或者說exon9是包含exon8的。那么應(yīng)該用sum-overlap的方法，sum好求，但overlap呢？

可以這樣：先不考慮具體長(zhǎng)度值，例如現(xiàn)在不直接讓exon1的長(zhǎng)度為90，而是輸出4807893 - 4807982這90個(gè)數(shù)值，目的就是使用unique()函數(shù)對(duì)這些數(shù)值去重復(fù)，最后一個(gè)length就求出來了

這個(gè)思路有了，而且這個(gè)過程一定是個(gè)循環(huán)，那么第一步就是：將對(duì)應(yīng)到同一個(gè)基因ID的外顯子都放一起

# 利用split(x,f)，需要一個(gè)x，一個(gè)f參數(shù)，其中x是向量或數(shù)據(jù)框，f是分組的因子。拆分完返回列表
split(t1,as.factor(t1$geneid))

然后第二步是對(duì)列表的每個(gè)元素取start到end的全部數(shù)值，也就是第二列到第三列，它返回的是一個(gè)列表

apply(x,1,function(y){y[2]:y[3]})

最后一步就是對(duì)列表去重、求長(zhǎng)度。

去重有兩種方法：
一是求總長(zhǎng)度，然后去overlap；二是先去overlap，再求總長(zhǎng)度
這里采用迂回式的第二種，也更容易操作
就是檢測(cè)y[2]:y[3]，發(fā)現(xiàn)有重復(fù)的數(shù)字計(jì)算一遍即可，因?yàn)檫@一個(gè)數(shù)字就代表一個(gè)堿基位點(diǎn)

length(unique(unlist(tmp)))

綜上，得到一個(gè)循環(huán)嵌套：

g_l = lapply(split(t1,t1$geneid),function(x){
    tmp=apply(x,1,function(y){
      y[2]:y[3]
    })
    length(unique(unlist(tmp)))
  })
# 再變成一個(gè)數(shù)據(jù)框
g_l=data.frame(gene_id=names(g_l),length=as.numeric(g_l))
> head(g_l)
    gene_id length
1 100009600   4352
2 100009609   2538
3 100009614    564
4 100009664   2398
5    100012   1854
6    100017   2736
# 結(jié)果將gene_id對(duì)應(yīng)到symbol
library(org.Mm.eg.db)
s2g=toTable(org.Mm.egSYMBOL)
g_l=merge(g_l,s2g,by='gene_id')

第二種：根據(jù)最長(zhǎng)轉(zhuǎn)錄本

這個(gè)就稍微方便一點(diǎn)，有一個(gè)求轉(zhuǎn)錄本長(zhǎng)度的函數(shù)：transcriptLengths

t_l=transcriptLengths(txdb)
>   head(t_l)
  tx_id    tx_name gene_id nexon tx_len
1     1 uc007afg.1   18777     8   2355
2     2 uc007afh.1   18777     9   2433
3     3 uc007afi.2   21399    10   2671
4     4 uc011wht.1   21399    10   2668
5     5 uc011whu.1   21399    10   2564
6     6 uc057aty.1    <NA>     1   2719
# 發(fā)現(xiàn)有NA，去掉即可
t_l=na.omit(t_l)
>   head(t_l)
  tx_id    tx_name gene_id nexon tx_len
1     1 uc007afg.1   18777     8   2355
2     2 uc007afh.1   18777     9   2433
3     3 uc007afi.2   21399    10   2671
4     4 uc011wht.1   21399    10   2668
5     5 uc011whu.1   21399    10   2564
7     7 uc007afm.2  108664     6   2396

然后看到gene_id這一列，有重復(fù)的ID，另外最后一列的length值也不同，說明這里一個(gè)基因有多個(gè)不同長(zhǎng)度的轉(zhuǎn)錄本

那么就對(duì)gene_id排序(為了讓同樣id的基因排在一起)，對(duì)tx_len排序(為了找最長(zhǎng)的轉(zhuǎn)錄本)

t_l=t_l[order(t_l$gene_id,t_l$tx_len,decreasing = T),]
> head(t_l)
      tx_id    tx_name gene_id nexon tx_len
15232 15232 uc008vih.2   99982    20   3075
15231 15231 uc008vig.2   99982    19   3015
9376   9376 uc008pkr.1   99929     6   4154
11784 11784 uc008rsk.3   99899     8   2909
11479 11479 uc008rat.3   99890     2   3065
11480 11480 uc008rau.2   99890     1   2475
# 上面??可以看到，99982這個(gè)基因就會(huì)選取第一個(gè)轉(zhuǎn)錄本15232，因?yàn)樗拈L(zhǎng)度為3075
# 其實(shí)這樣降序排列完，就可以直接去重，保留最長(zhǎng)的那一個(gè)(就是第一個(gè))了
t_l=t_l[!duplicated(t_l$gene_id),]
>   head(t_l)
      tx_id    tx_name gene_id nexon tx_len
15232 15232 uc008vih.2   99982    20   3075
9376   9376 uc008pkr.1   99929     6   4154
11784 11784 uc008rsk.3   99899     8   2909
11479 11479 uc008rat.3   99890     2   3065
10866 10866 uc008pqg.3   99889    10   3248
11552 11552 uc008rdv.2   99887     8   6271

有了基因長(zhǎng)度，計(jì)算RPKM

上面得到了有長(zhǎng)度信息的基因，首先要和原始矩陣的行名取交集

ng=intersect(rownames(a),g_l$symbol)
# 得到的新表達(dá)矩陣=》有長(zhǎng)度信息的基因表達(dá)矩陣
exprSet=a[ng,]
lengths=g_l[match(ng,g_l$symbol),2] #
> head(lengths)
[1] 1122  795  619 4556 1743 1471
> head(rownames(exprSet))
[1] "0610005C13Rik" "0610009B22Rik" "0610009L18Rik" "0610010F05Rik"
[5] "0610010K14Rik" "0610012G03Rik"

# 簡(jiǎn)單探索
> exprSet[1:3,1:3]
              SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5
0610005C13Rik              0              0              0
0610009B22Rik              0              0              0
0610009L18Rik              0              0              0
> dim(exprSet)
[1] 22731   768

接下來就是對(duì)新表達(dá)矩陣進(jìn)行操作：

先求文庫(kù)大小：

total_count<- colSums(exprSet)
> head(total_count)
SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5 SS2_15_0048_A4 
         95099          92537         161260         121297 
SS2_15_0048_A1 SS2_15_0048_A2 
        263927         280547 

然后需要對(duì)每個(gè)基因表達(dá)量值除以對(duì)應(yīng)的基因長(zhǎng)度（單位是Kb），再除以總文庫(kù)（單位是Mb）大小

# 如果說i代表一行表達(dá)量，exprSet[,i]就是表達(dá)矩陣的第一列，即22731個(gè)基因的表達(dá)量，lengths就是22731個(gè)基因的長(zhǎng)度，它和exprSet[,i]是一一對(duì)應(yīng)的。total_count[i]就是第一個(gè)樣本的文庫(kù)大小。最后需要乘以10^9來抵消單位的影響
10^9*exprSet[,i]/lengths/total_count[i]

再做一個(gè)循環(huán)：

實(shí)現(xiàn)了從第一個(gè)樣本到最后一個(gè)求得RPKM，結(jié)果返回一個(gè)包含768個(gè)元素的列表，每個(gè)樣本求得的RPKM值占一行

lapply(1:length(total_count),
                          function(i){
  10^9*exprSet[,i]/lengths/total_count[i]
})

接著對(duì)每個(gè)元素進(jìn)行rbind操作，拼接到一起，但是拼接的結(jié)果是：768行，22731列

do.call(rbind,
        lapply(1:length(total_count),
        function(i){
  10^9*exprSet[,i]/lengths/total_count[i]
}))

最后要轉(zhuǎn)置一下：

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]
     [,1] [,2] [,3]      [,4]
[1,]    0    0    0  7.347796
[2,]    0    0    0  0.000000
[3,]    0    0    0  0.000000
[4,]    0    0    0 19.904850

回顧總結(jié)

下面就以剛剛計(jì)算出來的第一行，第四列的7.347796這個(gè)值為例，看看到底R(shí)PKM是怎么算出來的，算是一個(gè)復(fù)習(xí)

首先看下表達(dá)矩陣：

> exprSet[1:4,1:4]
              SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5
0610005C13Rik              0              0              0
0610009B22Rik              0              0              0
0610009L18Rik              0              0              0
0610010F05Rik              0              0              0
              SS2_15_0048_A4
0610005C13Rik              1
0610009B22Rik              0
0610009L18Rik              0
0610010F05Rik             11

# 第一行第四列，也就是0610005C13Rik基因在SS2_15_0048_A4樣本中的原始表達(dá)量為1
# 然后總的文庫(kù)大小是121297
> total_count[4]
SS2_15_0048_A4 
        121297 
# 它的第一個(gè)基因0610005C13Rik長(zhǎng)度是1122
> lengths[1]
[1] 1122
# 因此計(jì)算公式就是
# 10^9*exprSet[,i]/lengths/total_count[i]
> 10^9*1/1122/121297
[1] 7.347796

從這個(gè)計(jì)算結(jié)果也可以看到：RPKM值為7，它的count值才為1，這樣會(huì)過濾掉很多存在RPKM表達(dá)量的基因，因此過濾基因設(shè)定count值為0就好

注意：不要認(rèn)為count值為1了，RPKM就是7左右。這個(gè)要取決于文庫(kù)大小，如果文庫(kù)很小，那么count為1時(shí)，RPKM也能達(dá)到幾萬(wàn)

最后和文章提供的RPKM矩陣比較下：

a=read.table('../GSE111229_Mammary_Tumor_fibroblasts_768samples_rpkmNormalized.txt.gz',
             header = T ,sep = '\t')
a[1:4,1:4]
rpkm_paper=a[ng,] 
# 文章做的
> rpkm_paper[1:4,1:4]
              SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5
0610005C13Rik              0              0              0
0610009B22Rik              0              0              0
0610009L18Rik              0              0              0
0610010F05Rik              0              0              0
              SS2_15_0048_A4
0610005C13Rik       6.966712
0610009B22Rik       0.000000
0610009L18Rik       0.000000
0610010F05Rik      19.843349
# 我們做的
> rpkm[1:4,1:4]
     [,1] [,2] [,3]      [,4]
[1,]    0    0    0  7.347796
[2,]    0    0    0  0.000000
[3,]    0    0    0  0.000000
[4,]    0    0    0 19.904850

可以看到有一點(diǎn)差別，但差別不大，這是因?yàn)橛?jì)算的基因長(zhǎng)度方法是有差別的