基因长度之多少

2019-04-02 本文已影响0人因地制宜的生信达人

基因长度之多少

无论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

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

基因长度之多少

基因长度之多少

代码

差异

猜你喜欢

热点阅读