前面提到基因組中的趣事（一）：這個(gè)基因編碼98種轉(zhuǎn)錄本，現(xiàn)在看看其它還有什么沒(méi)想到的？

序列最長(zhǎng)和最短的基因

計(jì)算基因序列的長(zhǎng)度，注意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

查看最長(zhǎng)和最短的3個(gè)基因

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

可變剪接調(diào)控基因RBFOX1以2.7 million的長(zhǎng)度超過(guò)之前文獻(xiàn)報(bào)道的最長(zhǎng)基因CNTNAP2 (智力語(yǔ)言損傷相關(guān)基因)。RBFOX1編碼的蛋白倒不長(zhǎng)，只有397個(gè)氨基酸，可見(jiàn)其內(nèi)含子區(qū)特別長(zhǎng)。

T細(xì)胞受體相關(guān)基因TRDD1作為最短的基因，長(zhǎng)度只有8 nt,編碼的小肽序列包含2個(gè)氨基酸 EI。

直接用上面的數(shù)據(jù)繪制長(zhǎng)度分布不太合適，拖尾很長(zhǎng)。對(duì)數(shù)據(jù)根據(jù)長(zhǎng)度區(qū)間重新做下分組 (當(dāng)然還可以分的再細(xì))，再進(jìn)行繪制其長(zhǎng)度分布。

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

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

只看蛋白編碼基因的長(zhǎng)度分布呢？蛋白編碼基因的長(zhǎng)度分布比較均勻，太短和太長(zhǎng)的都不多，集中在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

基因組中臨近基因最近和最遠(yuǎn)的是多少 （不考慮正負(fù)鏈）

# 需要考慮的是跨染色體的情況 # 只輸出不重疊的基因或只重疊1個(gè)堿基的基因 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

看下最近和最遠(yuǎn)的基因是什么？CTBP2P1和CCNQP2相距最遠(yuǎn)，距離有30個(gè)million。這是兩個(gè)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

距離最近的就是緊挨著了，主要是線粒體基因，串聯(lián)起來(lái)了。

# 前面做了排序，基因的順次就變了 # 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

包含外顯子最多的轉(zhuǎn)錄本

來(lái)一條代碼同時(shí)計(jì)算每個(gè)轉(zhuǎn)錄本外顯子的數(shù)目和每個(gè)外顯子的長(zhǎng)度。

# 第三列為exon # exoncnt 外顯子數(shù)量 # exonlen 每個(gè)轉(zhuǎn)錄本每個(gè)外顯子長(zhǎng)度 # 用到了二維數(shù)組。awk存儲(chǔ)二維數(shù)組時(shí)是用SUBSEP把兩個(gè)下標(biāo)連起來(lái)存儲(chǔ) # 索引取值時(shí)也需要先把兩個(gè)key切開(kāi)再取 # 結(jié)果存入兩個(gè)文件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

排個(gè)序看下，妊娠綜合征相關(guān)lncRNA HELLPAR的外顯子長(zhǎng)度最長(zhǎng)，超20萬(wàn) nt。外顯子長(zhǎng)度最長(zhǎng)的蛋白編碼基因是NFIA，一個(gè)轉(zhuǎn)錄因子，其外顯子長(zhǎng)度超4萬(wàn) nt。另外有33個(gè)基因各有一個(gè)長(zhǎng)度為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

包含外顯子數(shù)目最多的轉(zhuǎn)錄本是ENST00000589042，共有363個(gè)外顯子。其對(duì)應(yīng)的基因是TTN，橫紋肌發(fā)育相關(guān)基因。外顯子數(shù)目排第二的基因NEB也是骨骼肌微絲發(fā)育相關(guān)基因。肌肉的復(fù)雜性也許跟這些基因有關(guān)系，前面提到的最長(zhǎng)的基因、編碼轉(zhuǎn)錄本最多的基因也有一些是根肌肉發(fā)育相關(guān)的。

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怎么下載的？具體見(jiàn)推文NGS基礎(chǔ) - 參考基因組和基因注釋文件和下圖。

創(chuàng)作挑戰(zhàn)賽新人創(chuàng)作獎(jiǎng)勵(lì)來(lái)咯，堅(jiān)持創(chuàng)作打卡瓜分現(xiàn)金大獎(jiǎng)

總結(jié)

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

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。