基因的counts矩阵转换为RPKM矩阵

为什么要做这个计算

大家都知道在真核生物里面，一个基因有多个转录本，每个转录组又是由不同的外显子组合而成，以前老旧的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的公式难点就是这个基因的长度的求取，后面的没什么好讲的了。

赞 (0)

相关推荐