7

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

例子：

基因组.fasta

>sequence1
ATGGAGAGAGAGAGAGAGAGGGGGCAGCATACGCATCGACATACGACATACATCAGATACGACATACTACTACTATGA

注释.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'.

4

1 回答 1

1

此功能已添加到 scikit-bio，但 bioconda 中可用的版本尚不支持它（2017-12-15）。gff3 的格式文件存在于Github 存储库中。

您可以使用以下命令克隆 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)

对我来说，这会引发 a 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% 确定如何使用返回的间隔来提取特征。

于 2017-12-15T16:02:15.470 回答