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 序列,多个就连接起来, 负链则取反向互补序列进行连接储存,最终保存在字典里: