python 学习之提取 GTF 文件转录本序列
测试开头


测试结尾
学而时习之,不亦说乎

1引言
我们之前根据 GTF 文件筛选并 导出了最长转录本对应的 gtf 文件, 我们可以进一步利用这个文件,结合下载的 基因组序列文件 来 提取最长转录本的序列信息 。
当然也可以提取所有 gtf 文件里的转录本的序列。
2加入狗窝
欢迎加入我的狗窝, 更多资料, 知识, 代码, 技能 等你发现:

3读取基因组文件
我们知道对于 人和小鼠等基因组文件是很大的, 我们利用一般 逐行读取储存为字典 的方法操作 效率很慢 :
# 将基因组读取为字典
genome = {}
with open('GRCm39.primary_assembly.genome.fa','r') as genoemefa:
for line in genoemefa:
if line.startswith('>'):
name = line.replace('n','')
key = name[1:]
genome[key] = ''
else:
genome[key] += line.replace('n','')
上面要等好长一段时间,我们前几期介绍了 pyfastx 软件,可以快速读取进来:
import pyfastx
# 读取基因组
genome = pyfastx.Fasta('Mus_musculus.GRCm39.dna.primary_assembly.fa.gz')
# 提取序列
genome.fetch('1',(1,10))
'NNNNNNNNNN'
4提取序列
基本思路是根据 exon 来识别, 正链则根据位置从基因组来提取 exon 序列,多个就连接起来, 负链则取反向互补序列进行连接储存,最终保存在字典里:
请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!