基因组中的趣事（二）

2023-08-16 06:41| 来源: 网络整理| 查看: 265

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

【本文地址】

公司简介

联系我们