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”，“生信星球”的支持！
