counts2rpkm

Posted on by

为什么要做这个计算

大家都知道在真核生物里面，一个基因有多个转录本，每个转录组又是由不同的外显子组合而成，以前老旧的RNA-seq分析流程比较喜欢用RPKM值来量化表达量。大家很容易可以搜索到RPKM的公式，里面考虑到了基因长度。

而公式里面的基因长度并不是基因在染色体的起始终止坐标的差值，而通常是该基因的最长的转录本的外显子的长度之和。

所以计算起来就比较复杂，但是作为一个合格的生物信息学工程师，这种个性化分析需求要能hold住哈。

下面我就简单给出一个解决方案：

下载CCDS数据库文件

这里我们拿小鼠做例子

很容易谷歌得到CCDS的下载链接，下载查看文件格式如下：

wget ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_mouse/CCDS.current.txt 
head CCDS.current.txt 
#chromosome nc_accession    gene    gene_id ccds_id ccds_status cds_strand  cds_from    cds_to  cds_locations   match_type
1   NC_000067.6 Xkr4    497097  CCDS14803.1 Public  -   3216021 3671347 [3216021-3216967, 3421701-3421900, 3670551-3671347] Identical
1   NC_000067.6 Rp1 19888   CCDS14804.1 Public  -   4344599 4352824 [4344599-4350090, 4351909-4352080, 4352201-4352824] Identical
1   NC_000067.6 Sox17   20671   CCDS14805.1 Public  -   4491715 4493405 [4491715-4492667, 4493099-4493405]  Identical
​

第三列是基因名，第十列的中括号里面的是一个个外显子的起始终止坐标。这就是我们计算的依据。

perl单行命令处理

代码如下：

grep -v '^#' CCDS.current.txt | perl -alne '{/\[(.*?)\]/;$len=0;foreach(split/,/,$1){@tmp=split/-/;$len+=($tmp[1]-$tmp[0])};$h{$F[2]}=$len if $len >$h{$F[2]}} END{print "$_\t$h{$_}" foreach sort keys %h}' >mm10_ccds_length.txt
​

简单解释一下这个命令吧：

  • 首先把以#开头的注释信息去除 grep -v '^#'

  • 然后用正则匹配抓取中括号里面的信息 /\[(.*?)\]/

  • 接着对中括号里面的信息依据逗号进行分割成数组 split/,/,$1

  • 再接着对数组里面的外显子进行循环取长度并累加 @tmp=split/-/;$len+=($tmp[1]-$tmp[0])

  • 然后判断该转录本的外显子长度之和是否大于记录到该转录本对应的基因记录的长度，如果大于，就重新记录下来，因为我们要取最长转录本。 $h{$F[2]}=$len if $len >$h{$F[2]}

  • 最后把记录好的基因对应的最长转录本的外显子长度之和打印出来即可 print "$_\t$h{$_}" foreach sort keys %h}'

得到的结果如下：

tail mm10_ccds_length.txt 
Zwilch  1773
Zwint   752
Zxdb    2621
Zxdc    2567
Zyg11a  2107
Zyg11b  2221
Zyx 1686
Zzef1   8621
Zzz3    2722
a   393

小鼠里面居然有一个基因名字就是a,真搞笑。

是不是很简单呀，赶快练习一下吧！

也欢迎大家去我们论坛留言讨论。http://www.biotrainee.com/thread-1791-1-1.html

RPKM的公式难点就是这个基因的长度的求取，后面的没什么好讲的了。

大家自己点击下面的阅读原文即可去查看基因的counts矩阵如何转换为RPKM矩阵！

Comments are closed.