Bioperl push $ seq-> id to array

Question

Я довольно новичок в Perl и Bioperl, я пытаюсь написать скрипт, который будет идентифицировать экземпляры одинаковых последовательностей. Чтобы достичь этого, я представляю сценарий, который принимает 2 infiles, первый - множественное выравнивание в формате fasta, а второй - вспомогательный файл, который связывает идентификаторы fasta с другой соответствующей информацией. Мой подход состоит в том, чтобы прочитать множественное выравнивание с Bio :: SeqIO и поместить содержимое файла в хэш, где последовательность является ключом, а id - значением, или массив идентификаторов - это значение в случае совместного использования последовательности ,

Я думаю, что это должно выглядеть примерно так:

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

"TTTCTCTTTCCCAAAG" => 'SEQ2',

В в тот момент, я считаю, что я застреваю из-за ошибки при попытке нажать второй идентификатор на массив в случае совместного использования последовательностей (т. е. «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 -> Невозможно использовать строку («Seq1») в качестве массива ref в то время как «строгие ссылки» используются в строке ht_sharing.pl 24, строка 3.

2) Когда оператор печати за пределами цикла if активен, он печатает идентификатор, как я полагаю, должен (т. Seq1), но в выражении print внутри цикла if один и тот же вызов $ seq-> id вместо этого создает ссылку (то есть Bio :: Seq = HASH (0x19e7210) -> id). Почему это? Я не понимаю, почему печать $ seq-> id имеет разные выходы в одном и том же цикле while.

Я был бы очень благодарен, если кто-нибудь может предложить разъяснения и, конечно, как кто-то еще совершенно новый для этого замечания по лучшей практике или лучший способ подойти к проблеме, также замечательны.

Приветствия, Ana

Answer 1

Ваш код довольно близко, но есть несколько незначительных вопросов. Первый заключается в том, что вы хотите использовать синтаксис if (exists $hash{$key}) { ... }, чтобы увидеть, что существует ключ, defined сообщает, определено ли значение. Во-вторых, вы разыскиваете свой объект $seq без всякой причины.

Когда вы вызываете метод 'next_seq' в объекте Bio :: SeqIO, он возвращает объект Bio :: Seq. Если вы вызываете метод «id» в этом объекте Bio :: Seq, он возвращает идентификатор, как ожидалось, поэтому нет необходимости ничего уважать. Кроме того, нет необходимости напрямую импортировать Bio :: Seq (это просто комментарий, а не проблема).

Другие комментарии:

• Попробуйте поставить print Dumper %hts; вызов после while (my $seq ...) цикла (то есть, после того, как вы прошли через все объекты Seq). Сбрасывание хеша, когда вы просматриваете файл, не очень информативно.
• Вам действительно нужно сохранить идентификатор, связанный с последовательностью? Если вы просто хотите проверить повторяющиеся последовательности, попробуйте $hts{$seq->seq}++ в вашем цикле и загляните в отсортированные значения, чтобы увидеть, есть ли у вас дубликаты. Это будет быстрее.
• Это может быть упражнение, но учитывайте размер ваших данных. Если вы пытаетесь сделать это на сотнях миллионов последовательностей (очень распространенных в наши дни), вы можете долго ждать, а затем исчерпать память. Я упоминаю, что, поскольку есть другие способы сделать это, используя методы кластеризации, но я не хочу отговаривать вас за попытку решения BioPerl.