基因长度之多少

Posted on by

基因长度之多少

无论RPKM或FPKM多么的遭人诟病，它的真实需求还是存在， 那么我们该如何合理的定义基因的长度呢？

这个时候需要理解：基因，转录本(transcripts,isoform，mRNA序列)、EXON区域，cDNA序列、UTR区域，ORF序列、CDS序列 这些概念，一个基因可以转录为多个转录本，真核生物里面每个转录本通常是由一个或者多个EXON组成，能翻译为蛋白的EXON区域是CDS区域，不能翻译的那些EXON的开头和结尾是UTR区域，翻译区域合起来是ORF序列，而转录本逆转录就是cDNA序列。

有了这些概念，我们就能理解目前主流定义基因长度的几种方式。

  • 挑选基因的最长转录本
  • 选取多个转录本长度的平均值
  • 非冗余外显子(EXON)长度之和
  • 非冗余 CDS（Coding DNA Sequence） 长度之和

代码

这几种情况的代码如下：

library("TxDb.Mmusculus.UCSC.mm10.knownGene")
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
## 下面是定义基因长度为 非冗余exon长度之和
if(F){
 exon_txdb=exons(txdb)
 genes_txdb=genes(txdb)

 o = findOverlaps(exon_txdb,genes_txdb)
 o
 t1=exon_txdb[queryHits(o)]
 t2=genes_txdb[subjectHits(o)]
 t1=as.data.frame(t1)
 t1$geneid=mcols(t2)[,1]
 # 如果觉得速度不够，就参考R语言实现并行计算
 # http://www.bio-info-trainee.com/956.html
 g_l = lapply(split(t1,t1$geneid),function(x){
 # x=split(t1,t1$geneid)[[1]]
 head(x)
 tmp=apply(x,1,function(y){
 y[2]:y[3]
 })
 length(unique(unlist(tmp)))
 })
 head(g_l)
 g_l=data.frame(gene_id=names(g_l),length=as.numeric(g_l))

 save(g_l,file = 'step7-g_l.Rdata')
}
load(file = 'step7-g_l.Rdata')
## 下面是定义基因长度为 最长转录本长度
if(F){
 # 拿到所有基因的所有转录本长度信息
 t_l=transcriptLengths(txdb)
 head(t_l)
 t_l=na.omit(t_l)
 # 这里把同样的基因的多个转录本长度排序了
 t_l=t_l[order(t_l$gene_id,t_l$tx_len,decreasing = T),]
 str(t_l)
 # 冗余的那些转录本都是长度比较小的，直接去除即可
 t_l=t_l[!duplicated(t_l$gene_id),]
 head(t_l)
 g_l=t_l[,c(3,5)]

}
## 如果需要取同一个基因的多个转录本的长度平均值，上面代码稍微变化一下即可。
head(g_l)
library(org.Mm.eg.db)
s2g=toTable(org.Mm.egSYMBOL)
head(s2g)
g_l=merge(g_l,s2g,by='gene_id')

差异

参考counts2rpkm，定义基因长度为非冗余CDS之和 http://www.bio-info-trainee.com/3298.html

举例：

CCDS文件里面：

X NC_000086.7 Zxdb 668166 CCDS41065.1 Public + 94724674 94727295 [94724674-94727295] Identical

可以看到不同的定义，得到的基因长度定义还是差异很大的。

​

Comments are closed.