打开APP
userphoto
未登录

开通VIP,畅享免费电子书等14项超值服

开通VIP
为什么一个基因可以既是lncRNA又是protein_coding

有一些学生的Linux功底实在是太差了,所以我不得不重启六年前的《生信工程师》面试题给他们练习,有一个题目就是探索gtf。有意思的是学生们给我反馈了有几个基因居然既是lncRNA又是protein_coding。

首先去下载 gtf文件:

  • https://www.gencodegenes.org/human/
  • https://www.gencodegenes.org/mouse/

这个  gencode.v36.annotation.gtf.gz  文件也就是不到50M,所以很快就下载完毕,然后使用下面的代码格式化:

zcat gencode.v36.annotation.gtf.gz |  grep -v '_PAR_Y' |\
perl -alne  '{next unless  $F[1] eq "HAVANA";next unless $F[2] eq "gene";/gene_id \"(.*?)\.\d+\"; gene_type \"(.*?)\"; gene_name \"(.*?)\"/;print "$3\t$2\t$1\t$F[0]\t$F[3]\t$F[4]"}' \
> HAVANA_v36_human_gene_info

zcat gencode.v36.annotation.gtf.gz |  grep -v '_PAR_Y' |\
perl -alne  '{next unless  $F[1] ne "HAVANA";next unless $F[2] eq "gene";/gene_id \"(.*?)\.\d+\"; gene_type \"(.*?)\"; gene_name \"(.*?)\"/;print "$3\t$2\t$1\t$F[0]\t$F[3]\t$F[4]"}' \
> ENSEMBL_v36_human_gene_info

得到了两个基因信息文件, 简单的shell命令串起来统计一下,就知道 HAVANA和 ENSEMB 的区别了:

$ cut -f 2 ENSEMBL_v36_human_gene_info |sort |uniq -c|sort -nk1,1
      1 processed_pseudogene
      2 Mt_rRNA
      5 sRNA
      8 ribozyme
     18 pseudogene
     20 protein_coding
     22 Mt_tRNA
     47 rRNA
     49 scaRNA
    496 rRNA_pseudogene
    935 snoRNA
   1879 miRNA
   1899 snRNA
   2210 misc_RNA

上面的 ENSEMBL数据库里面都是一些比较少见的RNA基因,然后下面的 HAVANA数据库来源的基因要稍微正常一点 :

$ cut -f 2 HAVANA_v36_human_gene_info  |sort |uniq -c|sort -nk1,1
      1 IG_pseudogene
      1 misc_RNA
      1 scRNA
      1 translated_unprocessed_pseudogene
      1 vault_RNA
      2 snRNA
      2 translated_processed_pseudogene
      3 IG_J_pseudogene
      4 TR_D_gene
      4 TR_J_pseudogene
      6 TR_C_gene
      8 snoRNA
      9 IG_C_pseudogene
     14 IG_C_gene
     18 IG_J_gene
     33 TR_V_pseudogene
     37 IG_D_gene
     48 polymorphic_pseudogene
     79 TR_J_gene
     98 unitary_pseudogene
    106 TR_V_gene
    138 transcribed_unitary_pseudogene
    145 IG_V_gene
    187 IG_V_pseudogene
    500 transcribed_processed_pseudogene
    938 transcribed_unprocessed_pseudogene
   1057 TEC
   2612 unprocessed_pseudogene
  10159 processed_pseudogene
  16889 lncRNA
  19924 protein_coding

有意思的是HAVANA数据库来源的文件有53025,但是基因名字是有冗余的,只有52999行 :

wc -l *info
   7591 ENSEMBL_v36_human_gene_info
  53025 HAVANA_v36_human_gene_info
  60616 total
ubuntu 15:04:41 ~
$ cut -f 1 HAVANA_v36_human_gene_info |sort -u|wc
  52999   52999  464077

那我就简单探索一下是哪几个基因居然有冗余:


$ cut -f 1 HAVANA_v36_human_gene_info |sort |uniq -c |grep -v -w 1 
      2 ALG1L9P
      2 ANAPC1P2
      2 ARMCX5-GPRASP2
      2 BMS1P4
      2 CCDC39
      2 CLCA4-AS1
      2 DGCR5
      2 DNAJC9-AS1
      2 DUXAP8
      2 ELFN2
      2 GABARAPL3
      2 HERC2P7
      2 ITFG2-AS1
      2 LINC00484
      2 LINC00486
      2 MATR3
      2 PDE11A
      2 POLR2J3
      2 POLR2J4
      2 PRICKLE2-AS1
      2 RAET1E-AS1
      2 RGS5
      2 SFTA3
      2 SLFN12L
      2 SMIM40
      2 TMSB15B
ubuntu 15:05:06 ~
$ grep ELFN2 HAVANA_v36_human_gene_info
ELFN2 lncRNA ENSG00000243902 chr22 37339583 37427445
ELFN2 protein_coding ENSG00000166897 chr22 37367960 37427479

我随手查了一下,这个ELFN2基因,居然真的既是lncRNA又是protein_coding,可是我再认真看的时候,发现它其实是有两个完全不同的ensembl数据库ID,也就是说,其实是ID冲突,并不是这个基因有两副面孔哦!

检查全部的 ID 冲突的地方

使用下面的命令的组合:

cut -f 1 HAVANA_v36_human_gene_info |sort |uniq -c |grep -v -w 1 |awk '{print $2}' > id
$ grep -w -f id HAVANA_v36_human_gene_info 
## 得到的结果如下:
CLCA4-AS1 lncRNA ENSG00000236915 chr1 86571181 86696311
CLCA4-AS1 lncRNA ENSG00000261737 chr1 86703502 86704462
RGS5 protein_coding ENSG00000143248 chr1 163111121 163321791
RGS5-AS1 lncRNA ENSG00000232892 chr1 163161675 163213023
RGS5 lncRNA ENSG00000232995 chr1 163244505 163321894
LINC00486 lncRNA ENSG00000230876 chr2 32825359 32926693
LINC00486 lncRNA ENSG00000236854 chr2 32927085 32946149
ANAPC1P2 lncRNA ENSG00000285793 chr2 87030675 87076429
ANAPC1P2 unprocessed_pseudogene ENSG00000231259 chr2 87031815 87052992
PDE11A protein_coding ENSG00000128655 chr2 177623244 178072777
PDE11A protein_coding ENSG00000284741 chr2 177628069 178108339
PDE11A-AS1 lncRNA ENSG00000229941 chr2 177653419 177723289
PRICKLE2-AS1 lncRNA ENSG00000241111 chr3 64067964 64103131
PRICKLE2-AS1 lncRNA ENSG00000241572 chr3 64099273 64101122
CCDC39 protein_coding ENSG00000284862 chr3 180602858 180684942
CCDC39 lncRNA ENSG00000145075 chr3 180614008 180870933
CCDC39-AS1 lncRNA ENSG00000243187 chr3 180680084 180700449
MATR3 protein_coding ENSG00000280987 chr5 139273752 139331671
MATR3 protein_coding ENSG00000015479 chr5 139293674 139331677
SMIM40 protein_coding ENSG00000285064 chr6 33321386 33329286
SMIM40 protein_coding ENSG00000286920 chr6 33323628 33329279
RAET1E-AS1 lncRNA ENSG00000268592 chr6 149863494 149919507
RAET1E-AS1 lncRNA ENSG00000223701 chr6 149884431 149919508
POLR2J4 lncRNA ENSG00000214783 chr7 43940895 44019175
POLR2J4 transcribed_unprocessed_pseudogene ENSG00000272655 chr7 44013562 44019170
POLR2J3 protein_coding ENSG00000168255 chr7 102537918 102572653
POLR2J3 protein_coding ENSG00000285437 chr7 102562133 102572583
LINC00484 lncRNA ENSG00000229694 chr9 91118592 91182762
LINC00484 lncRNA ENSG00000235641 chr9 91159573 91166272
DNAJC9-AS1 lncRNA ENSG00000236756 chr10 73247360 73276984
DNAJC9-AS1 lncRNA ENSG00000227540 chr10 73252791 73254349
BMS1P4-AGAP5 lncRNA ENSG00000242288 chr10 73674295 73730466
BMS1P4 lncRNA ENSG00000271816 chr10 73699151 73730487
BMS1P4 transcribed_unprocessed_pseudogene ENSG00000242338 chr10 73715843 73730469
ALG1L9P lncRNA ENSG00000248671 chr11 71673885 71818238
ALG1L9P transcribed_unprocessed_pseudogene ENSG00000254978 chr11 71800541 71804640
ITFG2-AS1 lncRNA ENSG00000256150 chr12 2695765 2812902
ITFG2-AS1 lncRNA ENSG00000258325 chr12 2796877 2812902
SFTA3 protein_coding ENSG00000257520 chr14 36473207 36521149
SFTA3 lncRNA ENSG00000229415 chr14 36473288 36513829
HERC2P7 unprocessed_pseudogene ENSG00000281909 chr15 22480439 22484840
HERC2P7 transcribed_unprocessed_pseudogene ENSG00000274471 chr15 23309607 23314566
GABARAPL3 TEC ENSG00000279980 chr15 90347587 90349437
GABARAPL3 processed_pseudogene ENSG00000238244 chr15 90348844 90349197
SLFN12L protein_coding ENSG00000205045 chr17 35464249 35487857
SLFN12L lncRNA ENSG00000286065 chr17 35474904 35537861
DUXAP8 lncRNA ENSG00000206195 chr22 15784959 15829984
DUXAP8 transcribed_processed_pseudogene ENSG00000271672 chr22 15826566 15827187
DGCR5 lncRNA ENSG00000273032 chr22 18970439 19031242
DGCR5 transcribed_unprocessed_pseudogene ENSG00000237517 chr22 18985836 18994501
ELFN2 lncRNA ENSG00000243902 chr22 37339583 37427445
ELFN2 protein_coding ENSG00000166897 chr22 37367960 37427479
ARMCX5-GPRASP2 lncRNA ENSG00000271147 chrX 102599512 102714671
ARMCX5-GPRASP2 protein_coding ENSG00000286237 chrX 102712495 102753530
TMSB15B-AS1 lncRNA ENSG00000231728 chrX 103845151 103919548
TMSB15B protein_coding ENSG00000158427 chrX 103918896 103966712
TMSB15B protein_coding ENSG00000269226 chrX 104063871 104076212

在人类五六万基因里面就这么少量的几十个基因的冲突,已经是非常的幸运了,因为人类基因被研究的非常多,如果是其它物种,会更可怕,到处是各种各样的冲突

那么到底哪个ID更佳呢?

这里我选择去genecards数据库查询:

  • https://www.genecards.org/cgi-bin/carddisp.pl?gene=ELFN2 (数据库 :Ensembl: ENSG00000166897   )
  • https://www.genecards.org/Search/Keyword?queryString=ENSG00000243902 (数据库 :Lnc-ELFN2-2  )

可以看到,之所以这个ELFN2基因既是lncRNA又是protein_coding,其实是因为数据库ID的匹配出了问题, 另外一个基因名字应该是:Lnc-ELFN2-2 ,而不是ELFN2基因,这个ELFN2基因仍然是一个正常的蛋白编码基因?

一个生物学问题?

其实我都不知道这个问题算不算生物学问题,就是有没有真的一个基因它就既是lncRNA又是protein_coding,但并不是这种数据库ID匹配的失误造成的,而是它基因真实的特性呢?

上面的代码大量使用了Linux的shell命令技巧:

大家一定要把Linux的6个阶段跨越过去 ,一般来说,每个阶段都需要至少一天以上的学习:

  • 第1阶段:把linux系统玩得跟Windows或者MacOS那样的桌面操作系统一样顺畅,主要目的就是去可视化,熟悉黑白命令行界面,可以仅仅以键盘交互模式完成常规文件夹及文件管理工作。
  • 第2阶段:做到文本文件的表格化处理,类似于以键盘交互模式完成Excel表格的排序、计数、筛选、去冗余,查找,切割,替换,合并,补齐,熟练掌握awk,sed,grep这文本处理的三驾马车。
  • 第3阶段:元字符,通配符及shell中的各种扩展,从此linux操作不再神秘!
  • 第4阶段:高级目录管理:软硬链接,绝对路径和相对路径,环境变量。
  • 第5阶段:任务提交及批处理,脚本编写解放你的双手。
  • 第6阶段:软件安装及conda管理,让linux系统实用性放飞自我。
本站仅提供存储服务,所有内容均由用户发布,如发现有害或侵权内容,请点击举报
打开APP,阅读全文并永久保存 查看更多类似文章
猜你喜欢
类似文章
新版TCGAbiolinks包学习:表达矩阵提取(mRNA/lncRNA/counts/tpm/fpkm)
lncRNA的一些基础知识
Gene and transcript types
GTF文件格式简介
不同转录组流程结果到底该如何比较
TCGA lncRNA的提取 | 生信笔记
更多类似文章 >>
生活服务
热点新闻
分享 收藏 导长图 关注 下载文章
绑定账号成功
后续可登录账号畅享VIP特权!
如果VIP功能使用有故障,
可点击这里联系客服!

联系客服