有一些學生的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系統實用性放飛自我。