同样的基因在不同数据库差距好离谱

Posted on by

同样的基因在不同数据库差距好离谱

之前使用过 TxDb.Hsapiens.UCSC.hg38.knownGene 包来获取基因的坐标及长度，代码如下：

library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
genes_txdb=genes(txdb)
tmp=as.data.frame(genes_txdb)

从来没有想过，这个包居然有可能是错误的，我在找寻外显子序列最长的基因的时候发现了下面这些基因诡异的都很长，而且长度还很相似。

gene_id length
102724594 589655
7307 589655
875 582547
1409 513284
102724652 513026
102723553 284605

NCBI记录

首先我去NCBI搜索了7307 这个基因：https://www.ncbi.nlm.nih.gov/gene/7307
发现其GRCh38.p12 (GCF_000001405.38) 21 NC_000021.9 (43092956..43108291, complement)是这样的记录，意味着，其就15kb的长度，不可能其外显子序列高达589kb，所以中间肯定有什么误会。
我继续查看坐标记录是：chr21从6484623到43107578
然后我继续查看另外一个基因，https://www.ncbi.nlm.nih.gov/gene/102724594 在NCBI记录是21号染色体的NC_000021.9 (6484623..6499969, complement) 也不可能是589kb，在TxDb.Hsapiens.UCSC.hg38.knownGene包里面也被记录错了。

GENCODE记录

https://www.gencodegenes.org/ 算是人类研究比较权威的了，下载最新版gtf文件，搜索

wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_29/gencode.v29.annotation.gtf.gz

发现其坐标与NCBI记录相差不多啦。

U2AF1L5 protein_coding ENSG00000275895 chr21 6484623 6499261
U2AF1 protein_coding ENSG00000160201 chr21 43092956 43107587

为什么bioconductor包错了呢？

现在两个源都反对bioconductor的 TxDb.Hsapiens.UCSC.hg38.knownGene 包，所以它大概率应该是错了，可是为什么它错了呢？
好纠结。

高手解答

王铁柱的留言

这是因为这些基因在原本ucsc的数据库里就多个map position，比如875在txdb的原始区域就是两个，一个是6450000附近，一个是43068k附近，直接求width就是最大值减最小值，得到的其实是多个map之间的距离。至于为什么ucsc的refGene会出现多个区域，就是基因判定的问题了，ncbi认为cbsl和cbs就是后面那个位置上的一个基因，gencode注释的是前后两个不同基因

黄磊的留言

基因长度的定义一直有争议，因为一个基因可能有多个长短不一的transcript. UCSC有时把属于一个基因的距离很远的transcripts连在一起，有可能导致基因长度异常之大。genes()函数本身有些问题，所以一般是用最长的transcript的长度作为所属基因的长度

>library(TxDb.Hsapiens.UCSC.hg19.knownGene)
>txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
>tx_by_gene <- transcriptsBy(txdb, by="gene")
>gene_lens <- max(width(tx_by_gene))
gene_lens['7307']
​

Comments are closed.