Perl: вставить символ в хеш массива

Question

Я вычисляю баллы лог-шансов последовательностей и возвращаю мотив (небольшой фрагмент последовательности), который дает максимальный балл. У меня есть код, который вычисляет максимальную оценку для каждой последовательности в моем файле, и у меня возникают проблемы с сохранением мотива, который дает эту оценку. Смотрите мой другой пост (ы) для формата файла, общего вычисления оценок log-odds и т. Д. Perl: Создание и манипулирование хешем массивов для оценок log-odds последовательностей ДНК . Мой код выглядит следующим образом:

use strict;
use warnings;
use List::Util 'max';
use Data::Dumper; 

#USER SPECIFICATIONS
#User specifies motif width
my $width = 3;

#User enters the filename that contains the sequence data
print "Please enter the filename of the fasta sequence data: ";
my $filename1 = <STDIN>;

#Remove newline from file
chomp $filename1;

#Open the file and store each dna seq in hash
my %id2seq = ();
my %HoA = ();
my %loscore = ();
my %maxscore = ();
my %maxmot = ();
my $id = '';
open (FILE, '<', $filename1) or die "Cannot open $filename1.",$!;
my $dna;
while (<FILE>)
{
    if($_ =~ /^>(.+)/)
    {
        $id = $1; #Stores 'Sequence 1' as the first $id, for example
    }
    else
    {
        $HoA{$id} = [ split(//) ]; #Splits the contents to allow for position reference later
        $id2seq{$id} .= $_; #Creates a hash with each seq associated to an id number
        $maxmot{$id} = (); #Creates empty hash to push motifs to
        foreach $id (keys %HoA)
        {
            for my $len (0..(length($HoA{$id})-$width-1))
            {
                push @{ $loscore{$id} }, 0;
            }
        }
        push @{ $maxscore{$id} }, -30; #Creates a HoA with each id number to have a maxscore (initial score -30)
    }
}
close FILE;

foreach $id (keys %id2seq)
{
    print "$id2seq{$id}\n\n";
}
print "\n\n";



#EXPECTATION
#Create log-odds table of motifs
my %logodds;
$logodds{'A'}[0] = 0.1;
$logodds{'A'}[1] = 0.2;
$logodds{'A'}[2] = 0.3;
$logodds{'C'}[0] = 0.2;
$logodds{'C'}[1] = 0.5;
$logodds{'C'}[2] = 0.2;
$logodds{'G'}[0] = 0.3;
$logodds{'G'}[1] = 0.2;
$logodds{'G'}[2] = 0.4;
$logodds{'T'}[0] = 0.4;
$logodds{'T'}[1] = 0.1;
$logodds{'T'}[2] = 0.1;

#MAXIMIZATION
#Determine location for each sequence that maximally
#aligns to the motif pattern

foreach $id (keys %HoA)
{   
    for my $pos1 (0..length($HoA{$id})-$width-1)    #Look through all positions the motif can start at
    {
        for my $pos2 ($pos1..$pos1+($width-1))  #Define the positions within the motif (0 to width-1)
        {           
            for my $base (qw( A C G T))
            {
                if ($HoA{$id}[$pos2] eq $base)  #If the character matches a base:
                {
                    for my $pos3 (0..$width-1)  #Used for position reference in %logodds
                    {
                        #Calculate the log-odds score at each location
                        $loscore{$id}[$pos2] += $logodds{$base}[$pos3];

                        #Calculate the maximum log-odds score for each sequence

                        #Find the motif that gives the maximum score for each sequence
                        $maxscore{$id} = max( @{ $loscore{$id} });
                        if ($loscore{$id}[$pos2] == $maxscore{$id})
                        {
                            push @{ $maxmot{$id} }, $HoA{$id}[$pos3]; #NOT SURE THAT THIS IS THE CORRECT WAY TO PUSH WHAT I WANT
                        }
                    }
                }
            }
        }
    }
}

print "\n\n";
print Dumper(\%maxmot);

Ожидаемый результат для %maxmot должен быть примерно таким:

'Sequence 11' => [ 'C', 'A', 'T'], 'Sequence 14' => ['C', 'T', 'G'], etc.

В ожидаемом выводе должно быть три базы, потому что $width = 3. Вывод, который я получаю, дает мне кратные каждой базы, не в каком-либо заметном порядке (или, по крайней мере, я не могу заметить образец):

'Sequence 11' => [ 'C', 'C', 'C', 'A', 'A', 'A', 'A', 'T', 'A', 'A', 'T', 'T', 'T'], 'Sequence 14' => ['C', 'C', 'T', 'T', 'C', 'C', 'T', 'T', 'T', 'T', 'T'], etc. Я считаю, что проблема коренится в шаге push @{ $maxmot{$id} }, $HoA{$id}[$pos3];, но я не совсем уверен, как это исправить. Вместо этого я попытался нажать $HoA{$id}[$pos2], но это, похоже, тоже не решает мою проблему. Как всегда, любая помощь приветствуется! Я могу уточнить, если нужно, я знаю, что мой вопрос немного запутан. Заранее спасибо.

Answer 1

push() не проблема. После запуска вашего кода становится очевидным, что условное $loscore{$id}[$pos2] == $maxscore{$id} true чаще, чем вы ожидаете.

Вот несколько вопросов, которые я хотел бы задать в обзоре кода:

• почему вы используете length() для массива? Он не возвращает длину массива.
• Разве for my $base (qw( A C G T)) { if ($HoA{$id}[$pos2] eq $base) {... не является неэффективным способом эквивалента my $base = $HoA{$id}[$pos2];?
• вычисление для каждого $pos2 выполняется $pos2 + 1 раз, то есть один раз для 0, дважды для 1, ... то есть более поздние позиции получают более высокий балл.
• один расчет для $loscore{$id}[$pos2] является суммой @{ $logodds{$base} }, то есть база в позиции $pos2 + $pos3 игнорируется для расчета
• вы пересчитываете $maxscore{$id} при работе над последовательностями и , затем используете изменяющееся значение в условном выражении
• (моё предположение) мотив должен иметь длину $width баз, но ваш код хранит только одну базу в %maxmot

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

Я позволил себе переписать ваш код, чтобы быть более кратким и понятным. Но вы должны быть в состоянии отследить строки в новом коде до исходного кода.

#!/usr/bin/perl
use strict;
use warnings;

use List::Util 'max';
use Data::Dumper;

my $width = 3;

my %HoA;
my %maxpos;
my %loscore;
my $id = '';
while (<DATA>) {
    if (/^>(.+)/) {
        $id = $1;
    } else {
        $HoA{$id}     = [ split(//) ];
        $maxpos{$id}  = @{ $HoA{$id} } - $width - 1;
        $loscore{$id} = [ (0) x ($maxpos{$id} + 1) ];
    }
}

my %logodds = (
    A => [0.1, 0.2, 0.3],
    C => [0.2, 0.5, 0.2],
    G => [0.3, 0.2, 0.4],
    T => [0.4, 0.1, 0.1],
);

#MAXIMIZATION
my %maxscore;
my %maxmot;

# Calculate the log-odds score at each location
foreach $id (keys %HoA) {
    for my $index (0..$maxpos{$id}) {
        for my $offset (0..$width-1) {
            # look at base in sequence $id at $offset after $index
            my $base = $HoA{$id}[$index + $offset];
            $loscore{$id}[$index] += $logodds{$base}[$offset];
        }
    }
}

# Calculate the maximum log-odds score for each sequence
foreach $id (keys %HoA) {
    $maxscore{$id} = max( @{ $loscore{$id} });
}

# Find the motif that gives the maximum score for each sequence
foreach $id (keys %HoA) {
    for my $index (0..$maxpos{$id}) {
        if ($loscore{$id}[$index] == $maxscore{$id}) {
            # motif of length $width
            my $motif = join('', @{ $HoA{$id} }[$index..$index + $width - 1]);
            $maxmot{$id}->{$motif}++;
        }
    }
}

print Data::Dumper->Dump([\%loscore, \%maxscore, \%maxmot],
                         [qw(*loscore *maxscore *maxmot)]);

exit 0;

__DATA__
>Sequence_1
TCAGAACCAGTTATAAATTTATCATTTCCTTCTCCACTCCT
>Sequence_2
CCCACGCAGCCGCCCTCCTCCCCGGTCACTGACTGGTCCTG
>Sequence_3
TCGACCCTCTGGAACCTATCAGGGACCACAGTCAGCCAGGCAAG

Тестовый прогон:

$ perl dummy.pl
%loscore = (
             'Sequence_1' => [
                               '1.2',
                               '0.8',
                               '0.6',
                               '0.8',
                               '0.5',
                               '0.8',
                               '1',
                               '0.8',
                               '0.4',
                               '0.5',
                               '0.8',
                               '0.7',
                               '0.5',
                               '0.9',
                               '0.6',
                               '0.4',
                               '0.3',
                               '0.6',
                               '0.8',
                               '0.7',
                               '0.4',
                               '1.2',
                               '0.5',
                               '0.3',
                               '0.6',
                               '0.7',
                               '1.1',
                               '0.8',
                               '0.4',
                               '0.7',
                               '1',
                               '0.5',
                               '1.1',
                               '1',
                               '0.6',
                               '0.7',
                               '0.5',
                               '1.1',
                               '0.8'
                             ],
             'Sequence_2' => [
                               '0.9',
                               '1',
                               '0.6',
                               '1',
                               '0.6',
                               '1.1',
                               '0.8',
                               '0.5',
                               '1',
                               '1.1',
                               '0.6',
                               '1',
                               '0.9',
                               '0.8',
                               '0.5',
                               '1.1',
                               '0.8',
                               '0.5',
                               '1.1',
                               '0.9',
                               '0.9',
                               '1.1',
                               '0.8',
                               '0.6',
                               '0.6',
                               '1.2',
                               '0.6',
                               '0.7',
                               '0.7',
                               '0.9',
                               '0.7',
                               '0.7',
                               '0.7',
                               '1',
                               '0.6',
                               '0.6',
                               '1.1',
                               '0.8',
                               '0.7'
                             ],
             'Sequence_3' => [
                               '1.3',
                               '0.7',
                               '0.7',
                               '0.8',
                               '0.9',
                               '0.8',
                               '0.5',
                               '1',
                               '0.7',
                               '1',
                               '0.8',
                               '0.8',
                               '0.5',
                               '0.8',
                               '0.8',
                               '0.6',
                               '0.7',
                               '0.4',
                               '1.2',
                               '0.8',
                               '0.7',
                               '0.9',
                               '0.8',
                               '0.7',
                               '0.8',
                               '1',
                               '0.6',
                               '0.9',
                               '0.8',
                               '0.4',
                               '0.6',
                               '1.2',
                               '0.8',
                               '0.5',
                               '1',
                               '1',
                               '0.8',
                               '0.7',
                               '0.7',
                               '1.1',
                               '0.7',
                               '0.7'
                             ]
           );
%maxscore = (
              'Sequence_1' => '1.2',
              'Sequence_3' => '1.3',
              'Sequence_2' => '1.2'
            );
%maxmot = (
            'Sequence_3' => {
                              'TCG' => 1
                            },
            'Sequence_2' => {
                              'TCA' => 1
                            },
            'Sequence_1' => {
                              'TCA' => 2
                            }
          );

Это выглядит намного ближе к ожидаемому результату. Но, конечно, я мог быть совершенно не уверен в своих догадках ...

---

Если я правильно понимаю вычисление logscore, тогда оценка за мотив является константой и, следовательно, может быть предварительно рассчитана. Что привело бы к следующему более простому подходу:

#!/usr/bin/perl
use strict;
use warnings;

use Data::Dumper;

my $width = 3;

my %logodds = (
    A => [0.1, 0.2, 0.3],
    C => [0.2, 0.5, 0.2],
    G => [0.3, 0.2, 0.4],
    T => [0.4, 0.1, 0.1],
);

# calculate log score for each motif combination
my $motif_score = {'' => 0}; # start with a 0-length motif
foreach my $offset (0..$width - 1) {
    my %scores;

    # for all motifs...
    foreach my $prefix (keys %{ $motif_score }) {
        my $base_score = $motif_score->{$prefix};

        # ... add another base to the motif
        for my $base (qw(A G C T)) {
            $scores{"${prefix}${base}"} = $base_score + $logodds{$base}[$offset];
        }
    }

    # store the scores for the new sequences
    $motif_score = \%scores;
}

#print Data::Dumper->Dump([$motif_score], [qw(motif_score)]);

my $id;
my %maxmot;
while (<DATA>) {
    if (/^>(.+)/) {
        $id = $1;
    } else {
        chomp(my $sequence = $_);
        my $max = -1;

        # run a window of length $width over the sequence
        foreach my $index (0..length($sequence) - $width - 1) {

            # extract the motif at $index from sequence
            my $motif = substr($sequence, $index, $width);
            my $score = $motif_score->{$motif};

            # update max score if the motif has a higher score
            if ($score > $max) {
                $max         = $score;
                $maxmot{$id} = $motif;
            }
        }
    }
}

print Data::Dumper->Dump([\%maxmot], [qw(*maxmot)]);

exit 0;

__DATA__
>Sequence_1
TCAGAACCAGTTATAAATTTATCATTTCCTTCTCCACTCCT
>Sequence_2
CCCACGCAGCCGCCCTCCTCCCCGGTCACTGACTGGTCCTG
>Sequence_3
TCGACCCTCTGGAACCTATCAGGGACCACAGTCAGCCAGGCAAG

Тестовый прогон:

$ perl dummy.pl
%maxmot = (
            'Sequence_2' => 'TCA',
            'Sequence_3' => 'TCG',
            'Sequence_1' => 'TCA'
          );

---

Поскольку лог-лог по каждому мотиву является константой, список мотивов, отсортированный по порядку лог-логоров, также будет постоянным. Учитывая этот список, нам нужно будет только найти первый мотив, который соответствует данной последовательности. Следовательно, мы можем применить к этой задаче высоко оптимизированный механизм регулярных выражений. В зависимости от вашего реального размера проблемы это, вероятно, будет более эффективным решением:

#!/usr/bin/perl
use warnings;
use strict;

use List::Util qw(first pairs);
use Data::Dumper;

my $width = 3;

my %logodds = (
    A => [0.1, 0.2, 0.3],
    C => [0.2, 0.5, 0.2],
    G => [0.3, 0.2, 0.4],
    T => [0.4, 0.1, 0.1],
);
my @bases = keys %logodds;

# calculate log score for each motif combination
my $motif_logscore = {'' => 0}; # start with a 0-length motif
foreach my $offset (0..$width - 1) {
    my %score;

    # for all motifs...
    foreach my $prefix (keys %{ $motif_logscore }) {
        my $base_score = $motif_logscore->{$prefix};

        # ... add another base to the motif
        for my $base (@bases) {
            $score{"${prefix}${base}"} = $base_score + $logodds{$base}[$offset];
        }
    }

    # update hash ref to new motif scores
    $motif_logscore = \%score;
}

#print Data::Dumper->Dump([$motif_logscore], [qw(motif_logscore)]);

my @motifs_sorted =
    # list of [<motif>, <regular expression>] array refs
    map    { [$_->[0], qr/$_->[0]/] }
    # sort in descending numerical score order
    sort   { $b->[1] cmp $a->[1] }
    # list of [<motif>, <score>] array refs
    pairs %{ $motif_logscore };

#print Data::Dumper->Dump([\@motifs_sorted], [qw(*motifs_sorted)]);

my $id;
my %maxmot;
while (<DATA>) {
    if (/^>(.+)/) {
        $id = $1;
    } else {
        my $sequence = $_;
        # find the first pair where the regex matches -> store motif
        $maxmot{$id} = (
            first { ($sequence =~ $_->[1])[0] } @motifs_sorted
        )->[0];
    }
}

print Data::Dumper->Dump([\%maxmot], [qw(*maxmot)]);

exit 0;

__DATA__
>Sequence_1
TCAGAACCAGTTATAAATTTATCATTTCCTTCTCCACTCCT
>Sequence_2
CCCACGCAGCCGCCCTCCTCCCCGGTCACTGACTGGTCCTG
>Sequence_3
TCGACCCTCTGGAACCTATCAGGGACCACAGTCAGCCAGGCAAG

Answer 2

Может быть, вам нужен не массив, а хеш?

Изменить толчок на

undef $maxmot{$id}{ $HoA{$id}[$pos3] };

Создает хеш (с неопределенными значениями, поэтому важны только ключи). В выводе я не вижу более 3 клавиш для каждой последовательности во входных данных из вашего связанного вопроса.