2

我正在尝试使用以下代码从数据库中提取序列：

use strict;
use Bio::SearchIO; 
use Bio::DB::Fasta;


my ($file, $id, $start, $end) = ("secondround_merged_expanded.fasta","C7136661:0-107",1,10);
my $db = Bio::DB::Fasta->new($file);
my $seq = $db->seq($id, $start, $end);
print $seq,"\n";

我试图提取的序列的标题是：C7136661:0-107，如文件中所示：

>C7047455:0-100
TATAATGCGAATATCGACATTCATTTGAACTGTTAAATCGGTAACATAAGCAGCACACCTGGGCAGATAGTAAAGGCATATGATAATAAGCTGGGGGCTA

当我将标题切换到更标准的东西（如test）时，代码工作正常。我认为 BioPerl 不喜欢非标准标题。有什么办法可以解决这个问题，这样我就不必重新编码 FASTA 文件了？

4

2 回答 2

3

默认情况下，Bio::DB::Fasta将使用紧跟>在标题行上的所有非空格字符来形成序列的键。在您的情况下，这看起来像C7047455:0-100，这与子序列的内置缩写相同。如此处$db->seq($id, $start, $stop)所述，您可以使用,而不是$db->seq("$id:$start-$stop")调用 , 所以调用$db->seq('C7136661:0-107')看起来像您要求的那样$db->seq('C7136661', 0, 107)，并且该密钥不存在。

我无法知道您的数据中有什么，但是如果仅使用标题的第一部分直到冒号作为键就足够了，那么您可以使用-makeid回调来修改键。然后，如果您仅C7136661用于检索序列，它将起作用。

此代码演示。请注意，您可能已经有一个.index缓存文件，您必须在看到任何行为变化之前将其删除。

use strict;
use warnings;

use Bio::DB::Fasta;

my ($file, $id, $start, $end) = qw(
  secondround_merged_expanded.fasta
  C7136661
  1 10
);

my $db = Bio::DB::Fasta->new($file, -makeid => \&makeid);

sub makeid {
  my ($head) = @_;
  $head =~ /^>([^:]+)/ or die qq(Invalid header "$head");
  $1;
}

my $seq = $db->seq($id, $start, $end);
print $seq, "\n";
于 2012-12-04T20:57:09.263 回答
0

我有与这篇文章相关的问题。我想知道是否有人尝试过当查询中的位置超出 fasta 位置的限制时会发生什么。所以可以说，fasta 包含 100 个碱基，而您查询包含位置 102，此方法是否会捕获错误。我在一些真实数据中尝试了这个，它似乎总是返回“1”，但是，我的 fasta 序列包含 0/1，所以很难理解这是否是一个错误代码/它正在返回错误基数的输出。

我尝试查看文档，但找不到任何东西。

于 2014-02-04T23:11:13.387 回答