前面提到基因组中的趣事（一）：这个基因编码98种转录本，现在看看其它还有什么没想到的？
计算基因序列的长度，注意GTF中的位置是前闭后闭。
awk 'BEGIN{OFS=FS=" "}{if($3=="gene") {len1=$5-$4+1; print $10, $14, $18, len1;}}' GRCh38.tab.gtf | sort -k4,4nr | sed '1 i\ID Gene Type Length' >Gene_length查看最长和最短的3个基因
head -n 4 Gene_length; tail -n 3 Gene_lengthID Gene Type LengthENSG00000078328 RBFOX1 protein_coding 2473539ENSG00000174469 CNTNAP2 protein_coding 2304997ENSG00000153707 PTPRD protein_coding 2298478ENSG00000236597 IGHD7-27 IG_D_gene 11ENSG00000237235 TRDD2 TR_D_gene 9ENSG00000223997 TRDD1 TR_D_gene 8可变剪接调控基因RBFOX1以2.7 million的长度超过之前文献报道的最长基因CNTNAP2 (智力语言损伤相关基因)。RBFOX1编码的蛋白倒不长，只有397个氨基酸，可见其内含子区特别长。
T细胞受体相关基因TRDD1作为最短的基因，长度只有8 nt,编码的小肽序列包含2个氨基酸 EI。
直接用上面的数据绘制长度分布不太合适，拖尾很长。对数据根据长度区间重新做下分组 (当然还可以分的再细)，再进行绘制其长度分布。
awk 'BEGIN{OFS=FS=" "}{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 Length' >Gene_length.plot数据放入ImageGP，设置统一的bin大小为100。跟我最初的印象还是不太一样的，以前凭空以为1000-2000 nt的基因是最多的，实际看并不是，短基因还是更多的。
只看蛋白编码基因的长度分布呢？蛋白编码基因的长度分布比较均匀，太短和太长的都不多，集中在1000-5000 nt。
awk 'BEGIN{OFS=FS=" "}{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 Length' >PC_Gene_length.plot# 需要考虑的是跨染色体的情况# 只输出不重叠的基因或只重叠1个碱基的基因awk 'BEGIN{OFS=FS=" "; 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 FirstGene Dist' >Gene_dist.txt看下最近和最远的基因是什么？CTBP2P1和CCNQP2相距最远，距离有30个million。这是两个pseudogene。
head Gene_dist.txt; tail Gene_dist.txtSecondGene FirstGene DistCTBP2P1 CCNQP2 30228085AC242852.1 LINC01691 21762656BX088702.1 ABBA01045074.1 17301933ANKRD26P1 PPP1R1AP2 10556825ROCK1 RNU6-721P 5625962BX322784.1 KRT8P17 4793117EMB AC122694.1 4500222AC128676.1 AC023141.13 4439110AC069061.1 AC131274.2 4286863FIBP CTSW 0MTATP6P22 MTCO3P35 0MT-CO3 MT-ATP6 0MTCO3P10 AC099654.7 0MTCO3P12 MTATP6P1 0MTCO3P31 MTATP6P31 0MTND4P26 MTND4LP14 0MTND5P33 MTND6P33 0MT-TY MT-TC 0TRMT2A AC006547.2 0距离最近的就是紧挨着了，主要是线粒体基因，串联起来了。
# 前面做了排序，基因的顺次就变了# grep 'MT-' Gene_dist.txt# 重新算下，再捕获下awk 'BEGIN{OFS=FS=" "; 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 1MT-TV MT-RNR1 1MT-RNR2 MT-TV 1MT-TL1 MT-RNR2 1MT-ND1 MT-TL1 3MT-TI MT-ND1 1MT-TM MT-TQ 2MT-ND2 MT-TM 1MT-TW MT-ND2 1MT-TA MT-TW 8MT-TN MT-TA 2MT-TC MT-TN 32MT-TY MT-TC 0MT-CO1 MT-TY 13MT-TS1 MT-CO1 1MT-TD MT-TS1 4MT-CO2 MT-TD 1MT-TK MT-CO2 26MT-ATP8 MT-TK 2MT-CO3 MT-ATP6 0MT-TG MT-CO3 1MT-ND3 MT-TG 1MT-TR MT-ND3 1MT-ND4L MT-TR 1MT-TH MT-ND4 1MT-TS2 MT-TH 1MT-TL2 MT-TS2 1MT-ND5 MT-TL2 1MT-ND6 MT-ND5 1MT-TE MT-ND6 1MT-CYB MT-TE 5MT-TT MT-CYB 1MT-TP MT-TT 3来一条代码同时计算每个转录本外显子的数目和每个外显子的长度。
# 第三列为exon# exoncnt 外显子数量# exonlen 每个转录本每个外显子长度# 用到了二维数组。awk存储二维数组时是用SUBSEP把两个下标连起来存储# 索引取值时也需要先把两个key切开再取# 结果存入两个文件transcript_exon_cnt.txt# 和transcript_exon_len.txtawk 'BEGIN{OFS=FS=" "}{if($3=="exon") {trid=$20" "$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 205012KCNQ1OT1 ENST00000597346 1 91667AC003681.1 ENST00000624945 1 49287NFIA ENST00000603233 1 44880TSIX ENST00000604411 1 37027AC245452.1 ENST00000458178 2 36531FLRT2 ENST00000330753 2 33290SMAD2 ENST00000262160 11 32994PCDH9 ENST00000377861 2 27561GRIN2B ENST00000609686 13 27303包含外显子数目最多的转录本是ENST00000589042，共有363个外显子。其对应的基因是TTN，横纹肌发育相关基因。外显子数目排第二的基因NEB也是骨骼肌微丝发育相关基因。肌肉的复杂性也许跟这些基因有关系，前面提到的最长的基因、编码转录本最多的基因也有一些是根肌肉发育相关的。
TTN ENST00000589042 363TTN ENST00000591111 313TTN ENST00000342992 312TTN ENST00000615779 312TTN ENST00000342175 191TTN ENST00000359218 191TTN ENST00000460472 191NEB ENST00000618972 183NEB ENST00000397345 182NEB ENST00000427231 182GTF怎么下载的？具体见推文NGS基础 - 参考基因组和基因注释文件和下图。
| 留言与评论（共有 0 条评论） “” |
