基因組中的趣事（一）：這個基因編碼98種轉(zhuǎn)錄本
從ENSEMBL的注釋來看，人基因組中包含60,676個注釋的基因，19968個蛋白編碼基因。這些基因長度不同、位置不同、轉(zhuǎn)錄出的轉(zhuǎn)錄本不同，下面我們用幾篇推文一步步去了解下基因組中的基因都有哪些令我們驚訝的地方。
GFF全稱為general feature format，這種格式主要是用來注釋基因組中的基因信息。在推文NGS基礎(chǔ) - GTF/GFF文件格式解讀和轉(zhuǎn)換我們對這個格式做了詳細(xì)解釋。
基本結(jié)構(gòu)如下：
其最后一列為屬性列，包含的屬性信息可多可少，以ENSEMBL提供的人的GTF為例，包括基因的名字、ID和編碼信息等。
那么有了這個文件 (GRCh38.gtf)，我們能做些什么呢？
人GTF中注釋了多少種基因類型？
首先對GTF文件做個小處理，所有的雙引號"都替換為\t。
再利用下面的代碼組合確定每一列具體對應(yīng)什么信息，省卻了人工去數(shù)的麻煩 (代碼解釋見Linux學(xué)習(xí) - 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 ;轉(zhuǎn)換這一步比較耗時，直接轉(zhuǎn)一份存起來
sed 's/"/\t/g' GRCh38.gtf >GRCh38.tab.gtf提取并計數(shù)基因的類型
# 根據(jù)第三列選擇基因行
# 第18列為基因類型，進行計數(shù)
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的數(shù)目也很多，這是有點出乎意料的。（其它類型解釋見：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繪個圖吧，數(shù)據(jù)往高顏值免費在線繪圖工具 ImageGP一放，選擇Wide format，點擊提交，圖就出來了。
蛋白編碼基因最多能產(chǎn)生多少已經(jīng)注釋的轉(zhuǎn)錄本？
# 根據(jù)第三列選擇轉(zhuǎn)錄本行
# 根據(jù)類型選擇蛋白編碼的轉(zhuǎn)錄本
# 不知道哪一列是什么信息，用下面這句
# 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編碼轉(zhuǎn)錄本最多的是肢帶型肌營養(yǎng)不良相關(guān)基因SGCE，可以編碼98條轉(zhuǎn)錄本，這應(yīng)該很重要，在不同的情況下使用不同的剪接形式。
# 最后一列只是用來作圖的
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目前沒看到哪個網(wǎng)站可以一起查看多個基因的功能，ImageGP可以嘗試支持一下?，F(xiàn)在還是用命令來查找下吧，看上去也沒什么特別的，轉(zhuǎn)錄因子、G蛋白偶聯(lián)受體、鈣信號通路。PAX6是控制眼睛和其它感官發(fā)育的。SPG7是跨線粒體內(nèi)膜的3A基因。ANK2在心肌細(xì)胞特異高表達。
grep -f <(head Protein_coding_gene_transcript_count | tail -n+2 | cut -f 1) Human.anno.withalias.ortholog2.txt | cut -f 1-4
把數(shù)據(jù)拷貝到ImageGP，畫一個(bin=1)的直方圖 (也就是線圖了，省去了排序和計數(shù)了)，可以看到單個轉(zhuǎn)錄本的基因還是最多的。
一個顏色不好看，奇數(shù)和偶數(shù)顏色分開展示一下 （重新生成數(shù)據(jù)）
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個轉(zhuǎn)錄本的基因不多，合并下，省得看著拖尾 （直接在線修改參數(shù)）。
一個基因編碼多種不同類型的轉(zhuǎn)錄本
以轉(zhuǎn)錄本最多的基因SGCE (肢帶型肌營養(yǎng)不良相關(guān)基因)為例，其轉(zhuǎn)錄出4種不同類型的轉(zhuǎn)錄本。(processed_transcript: 總稱沒有開放閱讀框的轉(zhuǎn)錄本)
grep 'ENSG00000127990' GRCh38.tab.gtf | grep -P '\ttranscript\t' | cut -f 28 | sort | uniq -c
23 nonsense_mediated_decay
5 processed_transcript
45 protein_coding
25 retained_intron今天就先探索到這，下期繼續(xù) (實際是沒時間寫這么多了，分開寫吧，邊寫大家邊一起交流)。
Excel改變了你的基因名，30% 相關(guān)Nature文章受影響，NCBI也受波及
猜一猜最長的基因和最短的基因分別是多少？