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

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了利用gff提取某个基因的最长转录本(Python实现)相关的知识,希望对你有一定的参考价值。

参考技术A 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')

scikit-bio 从 gff3 文件中提取基因组特征

【中文标题】scikit-bio 从 gff3 文件中提取基因组特征【英文标题】:scikit-bio extract genomic features from gff3 file 【发布时间】:2016-07-11 07:59:43 【问题描述】:

scikit-bio 是否可以从基因组 fasta 文件中提取存储在 gff3 格式文件中的基因组特征?

例子:


基因组.fasta

>sequence1
ATGGAGAGAGAGAGAGAGAGGGGGCAGCATACGCATCGACATACGACATACATCAGATACGACATACTACTACTATGA

annotation.gff3

#gff-version 3
sequence1   source  gene    1   78  .   +   .   ID=gene1
sequence1   source  mRNA    1   78  .   +   .   ID=transcript1;parent=gene1
sequence1   source  CDS 1   6   .   +   0   ID=CDS1;parent=transcript1
sequence1   source  CDS 73  78  .   +   0   ID=CDS2;parent=transcript1

mRNA 特征 (transcript1) 的所需序列将是两个子 CDS 特征的串联。所以在这种情况下,这将是'ATGGAGCTATGA'

【问题讨论】:

从 scikit-bio 0.5.0 开始,不支持读取 gff3 文件。如果这是您希望添加到项目中的功能,请考虑在问题跟踪器上提交功能请求:github.com/biocore/scikit-bio/issues 【参考方案1】:

此功能已添加到 scikit-bio,但 bioconda 中可用的版本尚不支持(2017-12-15)。 gff3 的格式文件存在于Github repository。

您可以使用以下方法克隆 repo 并在本地安装它:

$ git clone https://github.com/biocore/scikit-bio.git
$ cd scikit-bio
$ python setup.py install

按照文件中给出的示例,以下代码应该可以工作:

import io
from skbio.metadata import IntervalMetadata
from skbio.io import read

gff = io.StringIO(open("annotations.gff3", "r").read())
im = read(gff, format='gff3', into=IntervalMetadata, seq_id="sequence1")

print(im)

对我来说,这会引发FormatIdentificationWarning,但条目报告正确:

4 interval features
-------------------
Interval(interval_metadata=<140154121000104>, bounds=[(0, 78)], fuzzy=[(False, False)], metadata='source': 'source', 'type': 'gene', 'score': '.', 'strand': '+', 'ID': 'gene1')
Interval(interval_metadata=<140154121000104>, bounds=[(0, 78)], fuzzy=[(False, False)], metadata='source': 'source', 'type': 'mRNA', 'score': '.', 'strand': '+', 'ID': 'transcript1', 'parent': 'gene1')
Interval(interval_metadata=<140154121000104>, bounds=[(0, 6)], fuzzy=[(False, False)], metadata='source': 'source', 'type': 'CDS', 'score': '.', 'strand': '+', 'phase': 0, 'ID': 'CDS1', 'parent': 'transcript1')
Interval(interval_metadata=<140154121000104>, bounds=[(72, 78)], fuzzy=[(False, False)], metadata='source': 'source', 'type': 'CDS', 'score': '.', 'strand': '+', 'phase': 0, 'ID': 'CDS2', 'parent': 'transcript1')

在代码中的示例中,GFF3 和 FASTA 文件在用于读取功能的输入字符串中连接。也许这可以解决这个问题。此外,我不是 100% 确定如何使用返回的间隔来提取特征。

【讨论】:

以上是关于利用gff提取某个基因的最长转录本(Python实现)的主要内容,如果未能解决你的问题,请参考以下文章

Corset轻松搞定无参转录组差异基因(转载)

基因组注释文件(GTF/GFF)格式介绍

scikit-bio 从 gff3 文件中提取基因组特征

转录组分析的正确姿势

NCBI获取基因序列以及不同转录本序列

基因组注释文件(二) gff 和 gtf文件格式说明