单细胞转录组学习笔记-12-RPKM概念及计算方法

2019-07-12 本文已影响35人刘小泽

刘小泽写于19.7.9+12-第二单元第十讲：RPKM概念及计算方法
笔记目的：根据生信技能树的单细胞转录组课程探索smart-seq2技术相关的分析技术
课程链接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53

RPKM须知

核心就在于基因长度的计算

参考这个网站，做的还是很清爽的：http://www.metagenomics.wiki/pdf/definition/rpkm-calculation

它的公式看上去很简单，就是对基因长度（也就是公式里的K）以及文库大小(也就是M) 进行标准化，其中文库大小很好理解，就是一个样本全部测序reads之和。但是这里有一个之前也没有注意到的问题，基因长度是什么？怎么计算出来的？

几个问题：它是简单的终止坐标减去起始坐标吗？它是外显子长度之和吗？那外显子之间有重合怎么处理？它是某一条转录本的长度吗？还是说基因长度是一个转录本的平均值？

这个问题是要好好总结一下，但首先还是要明确一些基因相关名词：

基因与转录本：一个基因可以转录成多个转录本
转录本与外显子：真核生物的每个转录本一般是由一个或多个外显子组成
转录本与cDNA：转录本逆转录得到cDNA
外显子与编码区CDS、非编码区UTR：可以翻译成蛋白的外显子区域是CDS区域，不能翻译的外显子开头、结尾是UTR区域
CDS与开放阅读框ORF：ORF是从起始密码子(ATG, but not always)到终止密码子(TAA, TAG, TGA)的DNA序列，另外又有两种链的方向，因此总共有6种阅读框，它也会包含内含子(这导致了真核生物的CDS与ORF不一致；另外在原核生物中它们是一样的)
http://www.scfbio-iitd.res.in/tutorial/orf.html

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

外显子和内含子是DNA层面的gene feature，密码子(codon)是RNA层面上的feature；外显子和内含子都存在于双链DNA（dsDNA）的蛋白编码基因区域中，一般通过sense strand，也即是5’=》3‘查看
然后dsDNA变成hnRNA (不均一核RNA：heterogeneous nuclear RNA)，和dsDNA的 5'-3' 方向一致，另外将T碱基换成了U。这里虽然hnRNA中也含有dsDNA的exon、intron序列，但称呼要换："exon transcripts"和"intron transcripts"
接着hnRNA进行了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

关于基因长度

既然基因包含这么多"组件"，那么求它的长度也会有几种方法：

选择最长的转录本
多个转录本的均值
非冗余外显子长度之和
非冗余CDS之和

注意到这里的"非冗余"，就是存在一个基因的多个外显子之间存在重叠（比如基因A的1号外显子较短，2号外显子长，1号包含在2号中），单纯的相加会重复计算

计算基因长度

第一种：非冗余外显子之和

载入原始表达矩阵

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

选择哪个包来计算基因长度呢？

不知道就搜索一下：

然后找到TxDb这个包的介绍，看Bioconductor这本电子书的41页 3.4.19节：https://bioconductor.github.io/BiocWorkshops/BioC2018.pdf

其中写道：

起步

因为处理的是小鼠数据，所以要加载相应小鼠的包：

library("TxDb.Mmusculus.UCSC.mm10.knownGene")
# 加载包以后，看看帮助文档，发现这么一句话(使用它只需要call它的名字)
## show the db object that is loaded by calling it's name

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

点进去发现，存在好多函数，因为我们想采用非冗余外显子加和的方法计算基因长度，因此选择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
# 看到结果是一个GRanges object，那么不懂就再搜索什么是GRanges
?GRanges
# The GRanges class is a container for the genomic locations and their associated annotations. 就是存储基因组坐标和相关注释信息的"容器"

因为取出的exon和gene都是坐标信息(可以看到上面的21062393-21075496就是坐标区间)，那么二者既然都是坐标而且不一致，就会有交叉，如何找到这个交叉？

又引入了一个新的函数=》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因为写在前面，于是它就作为queryHits；可以看到：它的第一个元素和genes_txdb的6729存在交叉，取出来看一下结果
> 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的第一个元素是chr1 4807893-4807982，对应上了基因的第6729个元素：chr1 4807893-4846735，于是这个exon就属于这个基因

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

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

t1=as.data.frame(t1) #更直观地来查看t1，结果将25w个外显子的坐标

为了给t1一个对应的基因ID，需要使用t2

t1$geneid=mcols(t2)[,1]
# mcols作用是提取 a DataFrame object containing the metadata columns
# 看下面👇这个，mcols(t2)[,1]提取的也就是虚线右侧的一列"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

再看一眼结果：

关键一步

好，看到基因18777对应着10个exon，那么如何求这个18777的长度呢？是简单把10个exon的第四列加起来吗？

不是的！请看：第8和第9个exon，它们是不是有重叠？或者说exon9是包含exon8的。那么应该用sum-overlap的方法，sum好求，但overlap呢？

可以这样：先不考虑具体长度值，例如现在不直接让exon1的长度为90，而是输出4807893 - 4807982这90个数值，目的就是使用unique()函数对这些数值去重复，最后一个length就求出来了

这个思路有了，而且这个过程一定是个循环，那么第一步就是：将对应到同一个基因ID的外显子都放一起

# 利用split(x,f)，需要一个x，一个f参数，其中x是向量或数据框，f是分组的因子。拆分完返回列表
split(t1,as.factor(t1$geneid))

然后第二步是对列表的每个元素取start到end的全部数值，也就是第二列到第三列，它返回的是一个列表

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

最后一步就是对列表去重、求长度。

去重有两种方法：
一是求总长度，然后去overlap；二是先去overlap，再求总长度
这里采用迂回式的第二种，也更容易操作
就是检测y[2]:y[3]，发现有重复的数字计算一遍即可，因为这一个数字就代表一个碱基位点

length(unique(unlist(tmp)))

综上，得到一个循环嵌套：

g_l = lapply(split(t1,t1$geneid),function(x){
    tmp=apply(x,1,function(y){
      y[2]:y[3]
    })
    length(unique(unlist(tmp)))
  })
# 再变成一个数据框
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
# 结果将gene_id对应到symbol
library(org.Mm.eg.db)
s2g=toTable(org.Mm.egSYMBOL)
g_l=merge(g_l,s2g,by='gene_id')

第二种：根据最长转录本

这个就稍微方便一点，有一个求转录本长度的函数：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
# 发现有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这一列，有重复的ID，另外最后一列的length值也不同，说明这里一个基因有多个不同长度的转录本

那么就对gene_id排序(为了让同样id的基因排在一起)，对tx_len排序(为了找最长的转录本)

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这个基因就会选取第一个转录本15232，因为它的长度为3075
# 其实这样降序排列完，就可以直接去重，保留最长的那一个(就是第一个)了
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

有了基因长度，计算RPKM

上面得到了有长度信息的基因，首先要和原始矩阵的行名取交集

ng=intersect(rownames(a),g_l$symbol)
# 得到的新表达矩阵=》有长度信息的基因表达矩阵
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"

# 简单探索
> 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

接下来就是对新表达矩阵进行操作：

先求文库大小：

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

然后需要对每个基因表达量值除以对应的基因长度（单位是Kb），再除以总文库（单位是Mb）大小

# 如果说i代表一行表达量，exprSet[,i]就是表达矩阵的第一列，即22731个基因的表达量，lengths就是22731个基因的长度，它和exprSet[,i]是一一对应的。total_count[i]就是第一个样本的文库大小。最后需要乘以10^9来抵消单位的影响
10^9*exprSet[,i]/lengths/total_count[i]

再做一个循环：

实现了从第一个样本到最后一个求得RPKM，结果返回一个包含768个元素的列表，每个样本求得的RPKM值占一行

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

接着对每个元素进行rbind操作，拼接到一起，但是拼接的结果是：768行，22731列

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

最后要转置一下：

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

回顾总结

下面就以刚刚计算出来的第一行，第四列的7.347796这个值为例，看看到底RPKM是怎么算出来的，算是一个复习

首先看下表达矩阵：

> 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样本中的原始表达量为1
# 然后总的文库大小是121297
> total_count[4]
SS2_15_0048_A4 
        121297 
# 它的第一个基因0610005C13Rik长度是1122
> lengths[1]
[1] 1122
# 因此计算公式就是
# 10^9*exprSet[,i]/lengths/total_count[i]
> 10^9*1/1122/121297
[1] 7.347796

从这个计算结果也可以看到：RPKM值为7，它的count值才为1，这样会过滤掉很多存在RPKM表达量的基因，因此过滤基因设定count值为0就好

注意：不要认为count值为1了，RPKM就是7左右。这个要取决于文库大小，如果文库很小，那么count为1时，RPKM也能达到几万

最后和文章提供的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

可以看到有一点差别，但差别不大，这是因为计算的基因长度方法是有差别的