还是使用CCDS数据库的基因坐标信息来计算基因长度吧

生信技能树

发布于 2022-07-26 09:52:54

8130

发布于 2022-07-26 09:52:54

文章被收录于专栏：生信技能树

我的一个GitHub项目里面的代码：https://github.com/jmzeng1314/scRNA_smart_seq2/blob/master/RNA-seq/step7-counts2rpkm.R ，有探索过3种方法获取基因长度，然后发现 同样的基因在不同数据库记录的位置信息差距好离谱 所以不得不弃用 TxDb.Hsapiens.UCSC.hg38.knownGene 包。

还是使用CCDS记录文件吧，CCDS 数据库旨在确定一组核心的人类和小鼠蛋白质编码区域，这些区域具有一致的注释和高质量。人类数据更新到了2018 年，包括了 33397 个 CCDS IDs，共 19033 个 Gene 。

在数据库：ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/ 可以下载，然后需要在Linux或者Mac环境下面使用 bedtools 软件加上perl代码，完成下面的操作。

首先是下载和查看 CCDS.current.txt 文件，有两个不同协议都可以下载，代码如下所示：

wget -c  ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_human/CCDS.current.txt  -O CCDS.current.ftp.txt
wget -c https://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_human/CCDS.current.txt -O CCDS.current.http.txt

不同地区网络不一样，建议反复确认哦。简单的看了看，两个协议下载的都是同一个文件：


$ wc -l CCDS.current.*
   35139 CCDS.current.ftp.txt
   35139 CCDS.current.http.txt

$ md5sum CCDS.current.*
72b9980c5fca7fc8b3127bb5381815a6  CCDS.current.ftp.txt
72b9980c5fca7fc8b3127bb5381815a6  CCDS.current.http.txt

cp  CCDS.current.http.txt CCDS.current.txt

一定要检查md5哦，保证文件下载完整。

该文件总共有 11 列，每一列的信息分别是

第1列：染色体
第2列：nc_accession ID
第3列：基因名，即 Symbol ID
第4列：gene_id，即 Entrez ID
第5列：ccds_id，后面的小数点为版本号
第6列：ccds_status 的注释状态，可以进行统计一下：

$ cat  CCDS.current.txt |sed '1d'| cut -f 6 | sort | uniq -c
  32582 Public
    249 Reviewed, update pending
    134 Reviewed, withdrawal pending
    331 Under review, update
    138 Under review, withdrawal
   1693 Withdrawn
     11 Withdrawn, inconsistent annotation

其中部分状态为Withdrawn 是由于缺乏完整的数据支撑。一般来说，后续分析可以直接过滤掉它。

第7列：正负链
第8列：CDS 起始坐标，需要注意的是，该坐标是 0-base 的标注方法
第9列：CDS 终止坐标，同上
第10列：CDS 中每个 exon 的具体坐标信息
第11列：match_type，匹配程度，指的是与 Ensembl and NCBI 等的注释是否一致

$ cat CCDS.current.txt | sed '1d' |cut -f 11 | sort |uniq -c
  35061 Identical
     63 None
     14 Partial

需要注意的是：Exon 包含 UTR 和 CDS 区域，而 CDS 仅仅指蛋白编码区域。

接下来我们简单的统计一下这个文件：



cat CCDS.current.txt |grep -v Withdrawn|wc -l 
   33435 
cat CCDS.current.txt |grep -v Withdrawn|cut -f3|sort |uniq |wc -l 
   19032 

cat CCDS.current.txt |grep -v Withdrawn|cut -f3|sort |uniq -c|sort -k1,1nr |head
    25 CACNA1G
     21 ST3GAL3
     20 CREM
     19 CACNA1C
     18 TCF4
     17 VEGFA
     15 CTNND1
     15 NR1I3
     14 CSF2RA
     14 DMKN

可以看到，这里面基本上包含了人类的接近2万个蛋白质编码基因，但是让我诧异的是蛮多基因都是有重复的记录。

# 需要安装 bedtools 软件
cat CCDS.current.txt |grep -v Withdrawn |perl -alne '{/\[(.*?)\]/;next unless $1;$gene=$F[2];$exons=$1;$exons=~s/\s//g;$exons=~s/-/\t/g;print "$F[0]\t$_\t$gene" foreach split/,/,$exons;}'|sort -u |bedtools sort -i >exon_probe.hg38.gene.bed
# 前面的15614个基因的非冗余外显子数量是 152616 个，都是在exon_probe.hg38.gene.bed 文件里面

cat exon_probe.hg38.gene.bed |perl -alne '{$s+=abs($F[2]-$F[1])}END{print $s}'
## 计算得到 WES 全长是 34673268 ， 约 35Mb，符合我们的认知哦
# 如果你对上面的perl语法不熟悉，看下面的awk语法也可以：
cat exon_probe.hg38.gene.bed  | awk 'BEGIN{cds_total=0} {cds_total=cds_total+($3-$2)} END{print cds_total}'

这里其实统计的是 CDS 总长度：前面提到了 Exon 包含 UTR 和 CDS 区域，而 CDS 仅仅指蛋白编码区域。所以 CDS 数据库不记录 UTR 信息的坐标。只能统计 CDS 区域的总长度，无法统计全外显子 Exon 的总长度，约 35 Mb。

可以看到比较长和比较短的基因分别是：

cat exon_probe.hg38.gene.bed|perl -alne '{$s{$F[3]}+=$F[2]-$F[1]}END{print "$_\t$s{$_}" foreach keys %s}' > gene_length.human.txt
 
 
sort -k2,2nr gene_length.human.txt |head
TTN     114068
MUC16   43440
OBSCN   27932
SYNE1   27005
NEB     25704
DST     21315
CCDC168 21242
SYNE2   21127
FSIP2   20701
ADGRV1  18831

sort -k2,2nr gene_length.human.txt |tail 
SLN     95
MTRNR2L4        86
RPL41   75
MTRNR2L10       74
MTRNR2L1        74
MTRNR2L3        74
MTRNR2L5        74
MTRNR2L6        74
MTRNR2L7        74
MTRNR2L8        74

我查了查这样的 MTRNR2L1 基因是 MT-RNR2 Like 1 (Pseudogene) ，不知道为什么这个 CCDS数据库要记录这些东西，奇奇怪怪的。

而 RPL41, a small ribosomal peptide 确实是很短，看到了一个文章提到：The RPL41 peptide (NH 2 ‐MRAKWRKKRMRRLKRKRRKMRQRSK‐OH) was synthesized by GenScript (Nanjing, China), 就是26个氨基酸，对应的是78个碱基，不知道为什么这个CCDS数据库记录它是75bp的长度。

写在文末

我在《生信技能树》，《生信菜鸟团》，《单细胞天地》的大量推文教程里面共享的代码都是复制粘贴即可使用的，有任何疑问欢迎留言讨论，也可以发邮件给我，详细描述你遇到的困难的前因后果给我，我的邮箱地址是 jmzeng1314@163.com

如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：

We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2022-03-08，如有侵权请联系 cloudcommunity@tencent.com 删除

数据库