Bioperl 将 $seq->id 推送到数组

Question

我对 Perl 和 Bioperl 还很陌生，我正在尝试编写一个脚本来识别相同序列的实例。为了实现这一点，我设想了一个包含 2 个 infile 的脚本，第一个是 fasta 格式的多重对齐，第二个是附件文件，它将 fasta id 链接到其他相关信息。我的方法是通读 Bio::SeqIO 的多重比对，并将文件内容放在哈希中，其中序列是键，id 是值，或者在序列共享的情况下，id 数组是值.

我认为它应该是这样的：

"AATTTGTTGTTGTACC" => ('Seq1', 'Seq13'),

"TTTCTCTTTCCCAAAG" => 'Seq2',

目前，我认为我被卡住了，因为在序列共享的情况下尝试将第二个 id 推送到数组中时出错（即上面示例中的“Seq13”）。

这是我正在使用的测试多重对齐：

    >Seq1
    AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
    >Seq2
    TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    >Seq13
    AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA

下面是我目前写的代码：

#!/usr/bin/perl
use strict;
use warnings;
use Bio::Seq;
use Bio::SeqIO;
use Data::Dumper;

my $seqs = shift @ARGV or die "please provide a multiple alignment file and an accesory information file: $!\n";
my $info = shift @ARGV or die "please provide a multiple alignment file and an accesory information file: $!\n";

#open(INFO, '<', $info);

my $inseq = Bio::SeqIO->new(
    -file => $seqs,
    -format => "fasta",
    );

my %hts;

while (my $seq = $inseq->next_seq) {
#    print $seq->seq(), "\t", $seq->id, "\n";
    if (defined $hts{$seq->seq()}) {
    print "Sequence already in hash:\t$seq->id\n";
    push @{$hts{$seq->seq()}}, ${$seq->id};
    }
    else {
    $hts{$seq->seq()} = $seq->id;
    }
    print Dumper \%hts
}

所以我希望得到一些帮助

1) 我收到一个我不太理解的错误，但我认为它与 push 语句有关 --> 在 ht_sharing.pl 第 24 行第 3 行使用“strict refs”时，不能使用字符串（“Seq1”）作为 ARRAY 引用。

2) 当 if 循环外的 print 语句处于活动状态时，它会打印我认为应该的 id（即 Seq1），但在 if 循环内的 print 语句中，相同的调用 $seq->id 会产生一个引用（即 Bio::Seq=HASH(0x19e7210)->id）。为什么是这样？我不明白为什么打印 $seq->id 在同一个 while 循环中有不同的输出。

如果有人能提供澄清，我将不胜感激，当然，因为在最佳实践或解决问题的更好方法方面对这个 cmets 还是很陌生的人也很好。

干杯， 安娜

Answer 1

您的代码非常接近，但有几个小问题。第一个是你想使用语法if (exists $hash{$key}) { ... } 来查看一个键是否存在，defined 告诉你该值是否被定义。第二件事是您无缘无故地取消引用您的 $seq 对象。

当您在 Bio::SeqIO 对象上调用方法“next_seq”时，它会返回一个 Bio::Seq 对象。如果你在那个 Bio::Seq 对象上调用 'id' 方法，它会按预期返回 ID，所以不需要遵守任何东西。此外，无需显式导入 Bio::Seq（这只是注释，不是问题）。

其他cmets：

• 尝试将您的print Dumper %hts; 调用放在while (my $seq ...) 循环之后（即，在您完成所有seq 对象之后）。在浏览文件时转储​​哈希值在这里不是很丰富。
• 您真的需要保留与序列关联的 ID 吗？如果您只想检查重复序列，请在循环中尝试$hts{$seq->seq}++，然后查看排序值以查看是否有重复。这样会更快。
• 这可能是一个练习，但请考虑数据的大小。如果您尝试在数亿个序列上执行此操作（这些天很常见），您可能会等待很长时间，然后执行此操作时内存不足。我提到这一点是因为还有其他方法可以做到这一点，使用集群方法，但我不想劝阻你尝试 BioPerl 解决方案。