是否可以在 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'
。
最佳答案
此功能已添加到 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% 确定如何使用返回的间隔来提取特征。
关于python - scikit-bio 从 gff3 文件中提取基因组特征,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/38301880/