序列最长和最短的基因
计算基因序列的长度，注意GTF中的位置是前闭后闭。
代码语言：javascript复制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 iIDtGenetTypetLength' >Gene_length查看最长和最短的3个基因
代码语言：javascript复制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。
直接用上面的数据绘制长度分布不太合适，拖尾很长。对数据根据长度区间重新做下分组 (当然还可以分的再细)，再进行绘制其长度分布。
代码语言：javascript复制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 iTypetLength' >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 iTypetLength' >PC_Gene_length.plot
基因组中临近基因最近和最远的是多少 （不考虑正负链）
代码语言：javascript复制# 需要考虑的是跨染色体的情况
# 只输出不重叠的基因或只重叠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 iSecondGenetFirstGenetDist' >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距离最近的就是紧挨着了，主要是线粒体基因，串联起来了。
代码语言：javascript复制# 前面做了排序，基因的顺次就变了
# 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包含外显子最多的转录本
来一条代码同时计算每个转录本外显子的数目和每个外显子的长度。
代码语言：javascript复制# 第三列为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 182GTF怎么下载的？具体见推文NGS基础 - 参考基因组和基因注释文件和下图。
