從ENSEMBL的注釋來看，人基因組中包含60,676個注釋的基因，19968個蛋白編碼基因。這些基因長度不同、位置不同、轉錄出的轉錄本不同，下面我們用幾篇推文一步步去了解下基因組中的基因都有哪些令我們驚訝的地方。

GFF全稱為general feature format，這種格式主要是用來注釋基因組中的基因信息。在推文NGS基礎 - GTF/GFF文件格式解讀和轉換我們對這個格式做了詳細解釋。

基本結構如下：

其最后一列為屬性列，包含的屬性信息可多可少，以ENSEMBL提供的人的GTF為例，包括基因的名字、ID和編碼信息等。

那么有了這個文件 (GRCh38.gtf)，我們能做些什么呢？

人GTF中注釋了多少種基因類型？

首先對GTF文件做個小處理，所有的雙引號"都替換為\t。

再利用下面的代碼組合確定每一列具體對應什么信息，省卻了人工去數的麻煩 (代碼解釋見Linux學習 - SED操作，awk的姊妹篇)。

# sed = 給每行加行號 # N 先讀進去一行，再讀一行，對2行同時進行操作 sed 's/"/\t/g' GRCh38.gtf | head -n 1 | tr '\t' '\n' | sed = | sed 'N;s/\n/\t/' 1 chr1 2 havana 3 gene 4 11869 5 14409 6 . 7 + 8 . 9 gene_id 10 ENSG00000223972 11 ; gene_version 12 5 13 ; gene_name 14 DDX11L1 15 ; gene_source 16 havana 17 ; gene_biotype 18 transcribed_unprocessed_pseudogene 19 ;

轉換這一步比較耗時，直接轉一份存起來

sed 's/"/\t/g' GRCh38.gtf >GRCh38.tab.gtf

提取并計數基因的類型

# 根據第三列選擇基因行 # 第18列為基因類型，進行計數 awk 'BEGIN{OFS=FS="\t"}{if($3=="gene") a[$18]+=1}END{ for(i in a) print i,a[i]}' GRCh38.tab.gtf | sort -k2,2nr | sed '1 i\Gene_type\tCount' | head

最多的是蛋白編碼基因近20000個，其次是lncRNA，pusudogene的數目也很多，這是有點出乎意料的。（其它類型解釋見：https://m.ensembl.org/info/genome/genebuild/biotypes.html 和 https://www.gencodegenes.org/pages/biotypes.html）

Gene_type Count protein_coding 19968 lncRNA 16880 processed_pseudogene 10168 unprocessed_pseudogene 2627 misc_RNA 2220

繪個圖吧，數據往高顏值免費在線繪圖工具 ImageGP一放，選擇Wide format，點擊提交，圖就出來了。

蛋白編碼基因最多能產生多少已經注釋的轉錄本？

# 根據第三列選擇轉錄本行 # 根據類型選擇蛋白編碼的轉錄本 # 不知道哪一列是什么信息，用下面這句 # sed -n '2p' GRCh38.tab.gtf | tr '\t' '\n' | sed = | sed 'N;s/\n/\t/'awk 'BEGIN{OFS=FS="\t"}{if($3=="transcript" && $22=="protein_coding") a[$10"\t"$18]+=1 }END{for(i in a) print i,a[i], "Gene"}' GRCh38.tab.gtf | sort -k3,3nr | sed '1 i\Gene\tSymbol\tTranscript_count\tGroup' >Protein_coding_gene_transcript_count

編碼轉錄本最多的是肢帶型肌營養不良相關基因SGCE，可以編碼98條轉錄本，這應該很重要，在不同的情況下使用不同的剪接形式。

# 最后一列只是用來作圖的 Gene Symbol Transcript_count Group ENSG00000127990 SGCE 98 Gene ENSG00000197912 SPG7 96 Gene ENSG00000145362 ANK2 94 Gene ENSG00000156113 KCNMA1 93 Gene ENSG00000196628 TCF4 91 Gene ENSG00000126091 ST3GAL3 85 Gene ENSG00000007372 PAX6 82 Gene ENSG00000205336 ADGRG1 79 Gene ENSG00000165795 NDRG2 78 Gene

目前沒看到哪個網站可以一起查看多個基因的功能，ImageGP可以嘗試支持一下。現在還是用命令來查找下吧，看上去也沒什么特別的，轉錄因子、G蛋白偶聯受體、鈣信號通路。PAX6是控制眼睛和其它感官發育的。SPG7是跨線粒體內膜的3A基因。ANK2在心肌細胞特異高表達。

grep -f <(head Protein_coding_gene_transcript_count | tail -n+2 | cut -f 1) Human.anno.withalias.ortholog2.txt | cut -f 1-4

把數據拷貝到ImageGP，畫一個(bin=1)的直方圖 (也就是線圖了，省去了排序和計數了)，可以看到單個轉錄本的基因還是最多的。

一個顏色不好看，奇數和偶數顏色分開展示一下 （重新生成數據）

awk 'BEGIN{OFS=FS="\t"}{if($3=="transcript" && $22=="protein_coding") a[$10"\t"$18]+=1 }END{for(i in a) {cnt=a[i]; type1="odd"; if(cnt%2==0) type1="even"; print cnt, type1}}' GRCh38.tab.gtf | sed '1 i\Transcript_count\tGroup' >Protein_coding_gene_transcript_count2.plot

多于50個轉錄本的基因不多，合并下，省得看著拖尾 （直接在線修改參數）。

一個基因編碼多種不同類型的轉錄本

以轉錄本最多的基因SGCE (肢帶型肌營養不良相關基因)為例，其轉錄出4種不同類型的轉錄本。(processed_transcript: 總稱沒有開放閱讀框的轉錄本)

grep 'ENSG00000127990' GRCh38.tab.gtf | grep -P '\ttranscript\t' | cut -f 28 | sort | uniq -c23 nonsense_mediated_decay5 processed_transcript45 protein_coding25 retained_intron

今天就先探索到這，下期繼續 (實際是沒時間寫這么多了，分開寫吧，邊寫大家邊一起交流)。

這些基因的名字太有才了，研究一下都可以發10分文章

Excel改變了你的基因名，30% 相關Nature文章受影響，NCBI也受波及

猜一猜最長的基因和最短的基因分別是多少？

總結

以上是生活随笔為你收集整理的基因组中的趣事（一）：这个基因编码98种转录本的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。