前面提到基因組中的趣事（一）：這個基因編碼98種轉錄本，現在看看其它還有什么沒想到的？

序列最長和最短的基因

計算基因序列的長度，注意GTF中的位置是前閉后閉。

awk 'BEGIN{OFS=FS="\t"}{if($3=="gene") {len1=$5-$4+1; print $10, $14, $18, len1;}}' GRCh38.tab.gtf | sort -k4,4nr | sed '1 i\ID\tGene\tType\tLength' >Gene_length

查看最長和最短的3個基因

head -n 4 Gene_length; tail -n 3 Gene_length ID Gene Type Length ENSG00000078328 RBFOX1 protein_coding 2473539 ENSG00000174469 CNTNAP2 protein_coding 2304997 ENSG00000153707 PTPRD protein_coding 2298478 ENSG00000236597 IGHD7-27 IG_D_gene 11 ENSG00000237235 TRDD2 TR_D_gene 9 ENSG00000223997 TRDD1 TR_D_gene 8

可變剪接調控基因RBFOX1以2.7 million的長度超過之前文獻報道的最長基因CNTNAP2 (智力語言損傷相關基因)。RBFOX1編碼的蛋白倒不長，只有397個氨基酸，可見其內含子區特別長。

T細胞受體相關基因TRDD1作為最短的基因，長度只有8 nt,編碼的小肽序列包含2個氨基酸 EI。

直接用上面的數據繪制長度分布不太合適，拖尾很長。對數據根據長度區間重新做下分組 (當然還可以分的再細)，再進行繪制其長度分布。

awk 'BEGIN{OFS=FS="\t"}{if($3=="gene") {len1=$5-$4+1; if(len1<2000) type1="Less2000"; else if(len1<5000) type1="Less5000"; else if(len1<20000) type1="Less20000"; else type1="Other"; print type1, len1;}}' GRCh38.tab.gtf | sed '1 i\Type\tLength' >Gene_length.plot

數據放入ImageGP，設置統一的bin大小為100。跟我最初的印象還是不太一樣的，以前憑空以為1000-2000 nt的基因是最多的，實際看并不是，短基因還是更多的。

只看蛋白編碼基因的長度分布呢？蛋白編碼基因的長度分布比較均勻，太短和太長的都不多，集中在1000-5000 nt。

awk 'BEGIN{OFS=FS="\t"}{if($3=="gene" && $18=="protein_coding") {len1=$5-$4+1; if(len1<2000) type1="Less2000"; else if(len1<5000) type1="Less5000"; else if(len1<20000) type1="Less20000"; else type1="Other"; print type1, len1;}}' GRCh38.tab.gtf | sed '1 i\Type\tLength' >PC_Gene_length.plot

基因組中臨近基因最近和最遠的是多少 （不考慮正負鏈）

# 需要考慮的是跨染色體的情況 # 只輸出不重疊的基因或只重疊1個堿基的基因 awk 'BEGIN{OFS=FS="\t"; lastgene=""; lastchr=""}{if(lastchr=="") lastchr=$1; if($3=="gene") {gene=$14 ;start=$4; end=$5; if($1!=lastchr) {lastchr=$1; lastgene=""; } if(lastgene!="") {dist=start-lastend; if(dist>=0) print gene,lastgene,dist; } lastgene=gene; lastend=end} }' GRCh38.tab.gtf | sort -k3,3nr | sed '1 i\SecondGene\tFirstGene\tDist' >Gene_dist.txt

看下最近和最遠的基因是什么？CTBP2P1和CCNQP2相距最遠，距離有30個million。這是兩個pseudogene。

head Gene_dist.txt; tail Gene_dist.txt SecondGene FirstGene Dist CTBP2P1 CCNQP2 30228085 AC242852.1 LINC01691 21762656 BX088702.1 ABBA01045074.1 17301933 ANKRD26P1 PPP1R1AP2 10556825 ROCK1 RNU6-721P 5625962 BX322784.1 KRT8P17 4793117 EMB AC122694.1 4500222 AC128676.1 AC023141.13 4439110 AC069061.1 AC131274.2 4286863 FIBP CTSW 0 MTATP6P22 MTCO3P35 0 MT-CO3 MT-ATP6 0 MTCO3P10 AC099654.7 0 MTCO3P12 MTATP6P1 0 MTCO3P31 MTATP6P31 0 MTND4P26 MTND4LP14 0 MTND5P33 MTND6P33 0 MT-TY MT-TC 0 TRMT2A AC006547.2 0

距離最近的就是緊挨著了，主要是線粒體基因，串聯起來了。

# 前面做了排序，基因的順次就變了 # grep 'MT-' Gene_dist.txt # 重新算下，再捕獲下 awk 'BEGIN{OFS=FS="\t"; lastgene=""; lastchr=""}{if(lastchr=="") lastchr=$1; if($3=="gene") {gene=$14 ;start=$4; end=$5; if($1!=lastchr) {lastchr=$1; lastgene=""; } if(lastgene!="") {dist=start-lastend; if(dist>=0) print gene,lastgene,dist; } lastgene=gene; lastend=end} }' GRCh38.tab.gtf | grep 'MT-' MT-RNR1 MT-TF 1 MT-TV MT-RNR1 1 MT-RNR2 MT-TV 1 MT-TL1 MT-RNR2 1 MT-ND1 MT-TL1 3 MT-TI MT-ND1 1 MT-TM MT-TQ 2 MT-ND2 MT-TM 1 MT-TW MT-ND2 1 MT-TA MT-TW 8 MT-TN MT-TA 2 MT-TC MT-TN 32 MT-TY MT-TC 0 MT-CO1 MT-TY 13 MT-TS1 MT-CO1 1 MT-TD MT-TS1 4 MT-CO2 MT-TD 1 MT-TK MT-CO2 26 MT-ATP8 MT-TK 2 MT-CO3 MT-ATP6 0 MT-TG MT-CO3 1 MT-ND3 MT-TG 1 MT-TR MT-ND3 1 MT-ND4L MT-TR 1 MT-TH MT-ND4 1 MT-TS2 MT-TH 1 MT-TL2 MT-TS2 1 MT-ND5 MT-TL2 1 MT-ND6 MT-ND5 1 MT-TE MT-ND6 1 MT-CYB MT-TE 5 MT-TT MT-CYB 1 MT-TP MT-TT 3

包含外顯子最多的轉錄本

來一條代碼同時計算每個轉錄本外顯子的數目和每個外顯子的長度。

# 第三列為exon # exoncnt 外顯子數量 # exonlen 每個轉錄本每個外顯子長度 # 用到了二維數組。awk存儲二維數組時是用SUBSEP把兩個下標連起來存儲 # 索引取值時也需要先把兩個key切開再取 # 結果存入兩個文件transcript_exon_cnt.txt # 和transcript_exon_len.txt awk 'BEGIN{OFS=FS="\t"}{if($3=="exon") {trid=$20"\t"$14; exoncnt[trid]+=1; exonlen[trid, exoncnt[trid]]=$5-$4+1}}END{transcript_exon_cnt="transcript_exon_cnt.txt"; for(i in exoncnt) print i, exoncnt[i] >transcript_exon_cnt; transcript_exon_len="transcript_exon_len.txt"; for(i in exonlen) {split(i, subkey, SUBSEP); print subkey[1], subkey[2], exonlen[subkey[1], subkey[2]] > transcript_exon_len;}}' GRCh38.tab.gtf

排個序看下，妊娠綜合征相關lncRNA HELLPAR的外顯子長度最長，超20萬 nt。外顯子長度最長的蛋白編碼基因是NFIA，一個轉錄因子，其外顯子長度超4萬 nt。另外有33個基因各有一個長度為1 nt的外顯子。

HELLPAR ENST00000626826 1 205012 KCNQ1OT1 ENST00000597346 1 91667 AC003681.1 ENST00000624945 1 49287 NFIA ENST00000603233 1 44880 TSIX ENST00000604411 1 37027 AC245452.1 ENST00000458178 2 36531 FLRT2 ENST00000330753 2 33290 SMAD2 ENST00000262160 11 32994 PCDH9 ENST00000377861 2 27561 GRIN2B ENST00000609686 13 27303

包含外顯子數目最多的轉錄本是ENST00000589042，共有363個外顯子。其對應的基因是TTN，橫紋肌發育相關基因。外顯子數目排第二的基因NEB也是骨骼肌微絲發育相關基因。肌肉的復雜性也許跟這些基因有關系，前面提到的最長的基因、編碼轉錄本最多的基因也有一些是根肌肉發育相關的。

TTN ENST00000589042 363 TTN ENST00000591111 313 TTN ENST00000342992 312 TTN ENST00000615779 312 TTN ENST00000342175 191 TTN ENST00000359218 191 TTN ENST00000460472 191 NEB ENST00000618972 183 NEB ENST00000397345 182 NEB ENST00000427231 182

GTF怎么下載的？具體見推文NGS基礎 - 參考基因組和基因注釋文件和下圖。

創作挑戰賽新人創作獎勵來咯，堅持創作打卡瓜分現金大獎

總結

以上是生活随笔為你收集整理的基因组中的趣事（二）- 最长的基因2.7 million，最短的基因只有8 nt却能编码的全部內容，希望文章能夠幫你解決所遇到的問題。

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