0

私は Perl と Bioperl にかなり慣れていません。同一のシーケンスのインスタンスを識別するスクリプトを作成しようとしています。これを実現するために、私は 2 つの infiles を使用するスクリプトを考えています。1 つ目は fasta 形式の複数のアラインメントで、2 つ目は fasta id を他の関連情報にリンクするアクセサリ ファイルです。私のアプローチは、Bio::SeqIO を使用して複数のアラインメントを読み取り、シーケンスがキーで ID が値であるハッシュにファイルの内容を配置することです。シーケンス共有の場合は ID の配列が値です。 .

私はそれが次のように見えるべきだと思います：

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

"TTTCTCTTTCCCAAAG" => 'Seq2',

現時点では、シーケンス共有の場合に 2 番目の 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) よくわからないエラーが表示されますが、プッシュ ステートメントに関連していると思われます --> ht_sharing で "strict refs" が使用されている間、文字列 ("Seq1") を ARRAY ref として使用できません。 pl 24 行目、3 行目。

2）ifループの外側のprintステートメントがアクティブな場合、私が信じているようにIDを出力します（つまり、Seq1）が、ifループ内のprintステートメントでは、同じ呼び出し$ seq-> idが代わりに参照を生成します（つまり、Bio ::Seq=HASH(0x19e7210)->id)。どうしてこれなの？$seq->id を印刷すると、同じ while ループ内で異なる出力が得られる理由がわかりません。

誰かが明確化を提供できれば本当に感謝しています。もちろん、ベストプラクティスや問題にアプローチするためのより良い方法について、このコメントにまだ慣れていない人も素晴らしいです。

乾杯、アナ

4

1 に答える 1

1

あなたのコードはかなり似ていますが、小さな問題がいくつかあります。最初のものは、構文を使用しif (exists $hash{$key}) { ... }てキーが存在するdefinedことを確認し、値が定義されているかどうかを示します。$seq2 つ目は、理由もなくオブジェクトを逆参照していることです。

Bio::SeqIO オブジェクトでメソッド 'next_seq' を呼び出すと、Bio::Seq オブジェクトが返されます。その Bio::Seq オブジェクトで 'id' メソッドを呼び出すと、期待どおりに ID が返されるため、何も参照する必要はありません。また、Bio::Seq を明示的にインポートする必要はありません (これは単なるコメントであり、問​​題ではありません)。

他のコメント：

  • print Dumper %hts;呼び出しをwhile (my $seq ...)ループの後 (つまり、すべての seq オブジェクトを処理した後)に置くようにしてください。ファイルを調べているときにハッシュをダンプすることは、ここではあまり有益ではありません。
  • シーケンスに関連付けられた ID を保持する必要がありますか? 重複したシーケンスを確認したいだけの場合$hts{$seq->seq}++は、ループを試して、並べ替えられた値を見て、重複があるかどうかを確認してください。その方が速くなります。
  • これは演習かもしれませんが、データのサイズを考慮してください。数億のシーケンスでこれを実行しようとすると (最近では非常に一般的です)、長い時間待機してからメモリ不足になる可能性があります。クラスタリング手法を使用してこれを行う方法は他にもあるため、これについて言及しますが、BioPerl ソリューションを試すことを思いとどまらせたくありません。
于 2013-11-11T23:56:02.480 に答える