前面提到基因组中的趣事（一）：这个基因编码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基础 - 参考基因组和基因注释文件和下图。