生信分析、meta分析、数据挖掘
TCGA、GEO、SEER、Oncomine

基因长度之多少

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

举例:

  • https://www.uniprot.org/uniprot/A2CE44

  • http://www.informatics.jax.org/marker/MGI:3694898

  • http://www.informatics.jax.org/sequence/marker/MGI:3694898?provider=RefSeq

  • https://asia.ensembl.org/Mus_musculus/Gene/Summary?g=ENSMUSG00000073062

  • https://www.ncbi.nlm.nih.gov/CCDS/CcdsBrowse.cgi?REQUEST=CCDS&DATA=CCDS41065.1

CCDS文件里面:

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

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

      ■


生信基础知识大全系列:生信基础知识100讲   

史上最强的生信自学环境准备课来啦!! 7次改版,11节课程,14K的讲稿,30个夜晚打磨,100页PPT的课程。   

如果需要组装自己的服务器;代办生物信息学服务器

如果需要帮忙下载海外数据(GEO/TCGA/GTEx等等),点我?

如果需要线下辅导及培训,看招学徒 

如果需要个人电脑:个人计算机推荐

如果需要置办生物信息学书籍,看:生信人必备书单

如果需要实习岗位:实习职位发布

如果需要售后:点我

如果需要入门资料大全:点我

赞(0) 打赏
未经允许不得转载:医学SCI科研之家 » 基因长度之多少
分享到: 更多 (0)

评论 抢沙发

评论前必须登录!

 

meta分析、生信分析

meta、生信交流群综合科研交流群