有一些学生的Linux功底实在是太差了,所以我不得不重启六年前的《生信工程师》面试题给他们练习,有一个题目就是探索gtf。有意思的是学生们给我反馈了有几个基因居然既是lncRNA又是protein_coding。
首先去下载 gtf文件:
这个 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冲突,并不是这个基因有两副面孔哦!
使用下面的命令的组合:
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
在人类五六万基因里面就这么少量的几十个基因的冲突,已经是非常的幸运了,因为人类基因被研究的非常多,如果是其它物种,会更可怕,到处是各种各样的冲突
这里我选择去genecards数据库查询:
可以看到,之所以这个ELFN2基因既是lncRNA又是protein_coding,其实是因为数据库ID的匹配出了问题, 另外一个基因名字应该是:Lnc-ELFN2-2 ,而不是ELFN2基因,这个ELFN2基因仍然是一个正常的蛋白编码基因?
其实我都不知道这个问题算不算生物学问题,就是有没有真的一个基因它就既是lncRNA又是protein_coding,但并不是这种数据库ID匹配的失误造成的,而是它基因真实的特性呢?
上面的代码大量使用了Linux的shell命令技巧:
大家一定要把Linux的6个阶段跨越过去 ,一般来说,每个阶段都需要至少一天以上的学习:
联系客服