python 提取基因非冗余 exon-5UTR-CDS-3UTR 长度
天空黑暗到一定程度,星辰便会熠熠生辉

1引言
之前有根据 featureCount 函数的原理来计算基因的 非冗余外显子的长度, 基于这个思想,我们还可以计算出 蛋白编码基因 的 5'UTR/CDS/3'UTR
三个区域的长度,用来做其他的分析。
前面刚学过 numpy 库,现学现用,我们使用 numpy 来储存每个特征的位置区间
来加快运行速度。
2操作
储存区间
我们用字典的值
来储存每个特征的区间位置, 测试 GTF 注释文件为 ensembl 数据库的。
查看所有特征类型,因为不同数据库命名不太一样:
$ awk '{print $3}' Homo_sapiens.GRCh38.101.gtf |sort -u
CDS
Selenocysteine
exon
five_prime_utr
gene
start_codon
stop_codon
three_prime_utr
transcript
储存不同特征的区间:
# 打开测试 gtf 文件
with open('Homo_sapiens.GRCh38.101.gtf','r') as gtf:
# 信息保存在字典里
info = {}
info_5UTR = {}
info_CDS = {}
info_3UTR = {}
for line in gtf:
# 跳过注释行
if line.startswith('#'):
continue
# 分割
fields = line.split()
# 类型
type = fields[2]
if type == 'exon':
# 名称
gene_name = fields[19].replace('"','').replace(';','')
gene_id = fields[9].replace('"','').replace(';','')
gene_type = fields[23].replace('"','').replace(';','')
# 连接名称
key = '|'.join([gene_name,gene_id,gene_type])
# 计算多个外显子长度
start = int(fields[3])
end = int(fields[4]) + 1
tmpfield = list(range(start,end))
# 储存所有exon位置信息
if info.get(key) is None:
info[key] = [range(start,end)]
else:
tmpfield = info[key]
tmpfield.append(range(start,end))
info[key] = tmpfield
else:
pass
if type == 'five_prime_utr':
# 名称
gene_name = fields[17].replace('"','').replace(';','')
gene_id = fields[9].replace('"','').replace(';','')
gene_type = fields[21].replace('"','').replace(';','')
# 连接名称
key = '|'.join([gene_name,gene_id,gene_type])
# 计算多个外显子长度
start = int(fields[3])
end = int(fields[4]) + 1
tmpfield = list(range(start,end))
# 储存所有5UTR位置信息
if info_5UTR.get(key) is None:
info_5UTR[key] = [range(start,end)]
else:
tmpfield = info_5UTR[key]
tmpfield.append(range(start,end))
info_5UTR[key] = tmpfield
else:
pass
if type == 'CDS':
# 名称
gene_name = fields[19].replace('"','').replace(';','')
gene_id = fields[9].replace('"','').replace(';','')
gene_type = fields[23].replace('"','').replace(';','')
# 连接名称
key = '|'.join([gene_name,gene_id,gene_type])
# 计算多个外显子长度
start = int(fields[3])
end = int(fields[4]) + 1
tmpfield = list(range(start,end))
# 储存所有CDS位置信息
if info_CDS.get(key) is None:
info_CDS[key] = [range(start,end)]
else:
tmpfield = info_CDS[key]
tmpfield.append(range(start,end))
info_CDS[key] = tmpfield
else:
pass
if type == 'three_prime_utr':
# 名称
gene_name = fields[17].replace('"','').replace(';','')
gene_id = fields[9].replace('"','').replace(';','')
gene_type = fields[21].replace('"','').replace(';','')
# 连接名称
key = '|'.join([gene_name,gene_id,gene_type])
# 计算多个外显子长度
start = int(fields[3])
end = int(fields[4]) + 1
tmpfield = list(range(start,end))
# 储存所有3UTR位置信息
if info_3UTR.get(key) is None:
info_3UTR[key] = [range(start,end)]
else:
tmpfield = info_3UTR[key]
tmpfield.append(range(start,end))
info_3UTR[key] = tmpfield
else:
pass
查看结果:
info['DDX11L1|ENSG00000223972|transcribed_unprocessed_pseudogene']
[range(11869, 12228),
range(12613, 12722),
range(13221, 14410),
range(12010, 12058),
range(12179, 12228),
range(12613, 12698),
range(12975, 13053),
range(13221, 13375),
range(13453, 13671)]
info_5UTR['OR4F5|ENSG00000186092|protein_coding']
[range(65419, 65434), range(65520, 65565), range(69055, 69091)]
info_CDS['OR4F5|ENSG00000186092|protein_coding']
[range(65565, 65574), range(69037, 70006), range(69091, 70006)]
info_3UTR['OR4F5|ENSG00000186092|protein_coding']
[range(70009, 71586), range(70009, 70109)]
区间取并集
我们把每个特征的区间使用 list 展开,在使用 numpy 的 去重 ,然后 len 计算长度即可。
先定义一个函数,方便对每个特征计算非冗余长度:
import numpy as np
# func
def getNonRedundentLength(range_info):
# save in dict
length_info = {}
# get union length
for key,val in range_info.items():
# innitiate np.array
arr = np.array([1])
for slice_range in val:
arr = np.append(arr,list(slice_range))
length_info[key] = len(np.unique(arr))-1
# reture result dict
return length_info
然后计算每个特征:
# 计算每个特征长度
exon_length = getNonRedundentLength(info)
UTR5_length = getNonRedundentLength(info_5UTR)
CDS_length = getNonRedundentLength(info_CDS)
UTR3_length = getNonRedundentLength(info_3UTR)
合并导出保存
将字典转为数据框:
import pandas as pd
# transorm into dataframe and merge by id
df_tran = pd.DataFrame.from_dict(exon_length,orient='index',columns=['translength'])
df_tran['ID'] = df_tran.index
df_5utr = pd.DataFrame.from_dict(UTR5_length,orient='index',columns=['utr5length'])
df_5utr['ID'] = df_5utr.index
df_cds = pd.DataFrame.from_dict(CDS_length,orient='index',columns=['cdslength'])
df_cds['ID'] = df_cds.index
df_3utr = pd.DataFrame.from_dict(UTR3_length,orient='index',columns=['urt3length'])
df_3utr['ID'] = df_3utr.index
df_tran

合并表格:
全部合并,包含所有基因:
# 按id合并表格
data_info = df_cds.merge(df_3utr.merge(df_tran.merge(df_5utr,on='ID',how='outer'),on='ID',how='outer'),on='ID',how='outer')
data_info

只保留含有 CDS 区的基因:
# 按id合并表格
data_info = df_cds.merge(df_3utr.merge(df_tran.merge(df_5utr,on='ID',how='outer'),on='ID',how='outer'),on='ID',how='left')
data_info

根据自己选择哪个结果。
保存:
# 添加基因名列
data_info['gene_name'] = [i.split(sep='|')[0] for i in data_info['ID']]
# 列排序
data_info = data_info[['ID','gene_name','translength','utr5length','cdslength','urt3length']]
# 缺失值替换为0
data_info.fillna(0,inplace = True)
# 保存
data_info.to_csv(r'All_transcripts_length_info.csv', index=False)
data_info

可以看到非编码基因的编码特征区域的长度都是 0。
3结尾
有几点需要注意的:
CDS 区包括起始密码子 ATG,但不包括终止密码子。 3UTR 区不包括终止密码子。 ensembl 数据库 stop codon 的三个碱基属于单独一个特征,其长度不包括在以上里面。
但是:
gencode 数据库的 GTF 文件 stop codon 包括在 UTR 区域里面,但其注释文件没有区分 UTR 类型。 gencode 数据库的 蛋白编码转录本序列
文件中, stop codon 包括在 CDS 区域里面。
感兴趣的大家可以自己查看研究一下。

欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群
哦,数据代码已上传至QQ群,欢迎加入下载。
群二维码:

老俊俊微信:
知识星球:
所以今天你学习了吗?
欢迎小伙伴留言评论!
今天的分享就到这里了,敬请期待下一篇!
最后欢迎大家分享转发,您的点赞是对我的鼓励和肯定!
如果觉得对您帮助很大,赏杯快乐水喝喝吧!
往期回顾
◀python bioinfokit 计算 FPKM/RPKM/TPM 值
◀...
请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!