利用gff提取某个基因的最长转录本（Python实现）

Python018

利用gff提取某个基因的最长转录本（Python实现）,第1张

from Bio import SeqIO

from Bio.Seq import Seq

from Bio.SeqRecord import SeqRecord

import gffutils

import re

#######################

#从cds获得每条转录本的长度，从gff获得每个基因包含哪些转录本

fout = 'output/longest.fa'

cds = 'data/rna.fna'

gff = 'data/genomic.gff'

fout = 'output/longest.fa'

seq_index = SeqIO.index(cds ,'fasta')

gffdb = gffutils.create_db(gff,'gff.db'.force=True,

                          merge_strategy='creat_unique')

#以基因为key，转录本为value的字典

gene2longestcds = {}

#对基因进行迭代

for g in gffdb.all_features(featuretype='gene'):

  g_id= g.id

  for m in gffdb.children(g, featuretype='mRNA'):##只关注mRNA信息

    m_id = m.id.replace('rna-' , '')

    m_len = len(seq_index[m_id].seq)

    #一个基因中最长转录本的名字

    if g_id not in gene2longestcds:

      gene2longestcds[g_id] = m_id

    elif m_len >len(seq_index[gene2longestcds[g_id]].seq):

      gene2longestcds[g_id] = m_id

sr_list = [v for k,v in seq_index.items() if k in gene2longestcds.values()]

##输出

SeqIO.write(sr_list, fout, 'fasta')

re模块是Python中的正则表达式调用模块，在python中，通过将正则表达式内嵌集成re模块，程序员们可以直接调用来实现正则匹配。

正则表达式的大致匹配过程是：

re模块所支持的方法有如下：

其中，pattern为匹配模式，由re.compile生成，例如： pattern = re.compile(r'hello')

参数flag是匹配模式，取值可以使用按位或运算符’|’表示同时生效，如 re.I | re.M。可选值有：

注：以上七个方法中的flags同样是代表匹配模式的意思，如果在pattern生成时已经指明了flags，那么在下面的方法中就不需要传入这个参数了。

1. re.match(pattern, string[, flags])

2. re.search(pattern, string[, flags])

3. re.split(pattern, string[, maxsplit])

4. re.findall(pattern, string[, flags])

5. re.finditer(pattern, string[, flags])

6. re.sub(pattern, repl, string[, count])

7. re.subn(pattern, repl, string[, count])

描述

注意：该方法只能删除开头或是结尾的字符，不能删除中间部分的字符。

语法

参数

返回值

实例

以上实例输出结果如下：

参考

https://www.cnblogs.com/chengege/p/11190782.html

https://www.runoob.com/python/att-string-strip.html