8

我开始深入研究 Perl，但是在编写“Perl-ly”代码而不是用 Perl 编写 C 时遇到了麻烦。如何更改以下代码以使用更多 Perl 习语，我应该如何学习这些习语？

只是对它在做什么的解释：此例程是对齐 DNA 或氨基酸序列的模块的一部分（如果您关心此类事情，请使用 Needelman-Wunch）。它创建了两个二维数组，一个用于存储两个序列中每个位置的分数，另一个用于跟踪路径，以便以后可以重新创建得分最高的对齐方式。它工作得很好，但我知道我做的事情不是很简洁明了。

编辑：这是为了分配。我完成了它，但想稍微清理一下我的代码。如果您有兴趣，可以在课程网站上找到有关实现该算法的详细信息。

sub create_matrix {
    my $self = shift;
    #empty array reference
    my $matrix = $self->{score_matrix};
    #empty array ref
    my $path_matrix = $self->{path_matrix};
    #$seq1 and $seq2 are strings set previously
    my $num_of_rows = length($self->{seq1}) + 1;
    my $num_of_columns = length($self->{seq2}) + 1;

    #create the 2d array of scores
    for (my $i = 0; $i < $num_of_rows; $i++) {
        push(@$matrix, []);
        push(@$path_matrix, []);
        $$matrix[$i][0] = $i * $self->{gap_cost};
        $$path_matrix[$i][0] = 1;
    }

    #fill out the first row
    for (my $i = 0; $i < $num_of_columns; $i++) {
        $$matrix[0][$i] = $i * $self->{gap_cost};
        $$path_matrix[0][$i] = -1;
    }
    #flag to signal end of traceback
    $$path_matrix[0][0] = 2;
    #double for loop to fill out each row
    for (my $row = 1; $row < $num_of_rows; $row++) {
        for (my $column = 1; $column < $num_of_columns; $column++) {
            my $seq1_gap = $$matrix[$row-1][$column] + $self->{gap_cost};
            my $seq2_gap = $$matrix[$row][$column-1] + $self->{gap_cost};
            my $match_mismatch = $$matrix[$row-1][$column-1] + $self->get_match_score(substr($self->{seq1}, $row-1, 1), substr($self->{seq2}, $column-1, 1));
            $$matrix[$row][$column] = max($seq1_gap, $seq2_gap, $match_mismatch);

            #set the path matrix
            #if it was a gap in seq1, -1, if was a (mis)match 0 if was a gap in seq2 1
            if ($$matrix[$row][$column] == $seq1_gap) {
                $$path_matrix[$row][$column] = -1;
            }
            elsif ($$matrix[$row][$column] == $match_mismatch) {
                $$path_matrix[$row][$column] = 0;
            }
            elsif ($$matrix[$row][$column] == $seq2_gap) {
                $$path_matrix[$row][$column] = 1;
            }
        }
    }
}
4

6 回答 6

9

您收到了一些关于语法的建议，但我也建议采用更模块化的方法，如果没有其他原因代码可读性。如果您能够在担心低级细节之前了解大局，那么加快代码速度会容易得多。

您的主要方法可能如下所示。

sub create_matrix {
    my $self = shift;
    $self->create_2d_array_of_scores;
    $self->fill_out_first_row;
    $self->fill_out_other_rows;
}

你也会有几个像这样的小方法：

n_of_rows
n_of_cols
create_2d_array_of_scores
fill_out_first_row
fill_out_other_rows

并且您可以通过定义更小的方法（getter、setter 等）来更进一步。那时，您的中级方法create_2d_array_of_scores根本不会直接触及底层数据结构。

sub matrix      { shift->{score_matrix} }
sub gap_cost    { shift->{gap_cost}     }

sub set_matrix_value {
    my ($self, $r, $c, $val) = @_;
    $self->matrix->[$r][$c] = $val;
}

# Etc.
于 2009-10-23T17:27:00.020 回答
8

一个简单的改变是使用for这样的循环：

for my $i (0 .. $num_of_rows){
    # Do stuff.
}

有关详细信息，请参阅有关foreach 循环和范围运算符的 Perl 文档。

于 2009-10-23T16:32:27.620 回答
7

我还有其他一些评论，但这是第一个观察结果：

my $num_of_rows = length($self->{seq1}) + 1;
my $num_of_columns = length($self->{seq2}) + 1;

字符串也是如此$self->{seq1}，$self->{seq2}您可以继续使用substr. 我更愿意将它们存储为字符数组：

$self->{seq1} = [ split //, $seq1 ];

这是我写它的方式：

sub create_matrix {
    my $self = shift;

    my $matrix      = $self->{score_matrix};
    my $path_matrix = $self->{path_matrix};

    my $rows = @{ $self->{seq1} };
    my $cols = @{ $self->{seq2} };

    for my $row (0 .. $rows) {
        $matrix->[$row]->[0] =  $row * $self->{gap_cost};
        $path_matrix->[$row]->[0] = 1;
    }

    my $gap_cost = $self->{gap_cost};

    $matrix->[0] = [ map { $_ * $gap_cost } 0 .. $cols ];
    $path_matrix->[0] = [ (-1) x ($cols + 1) ];

    $path_matrix->[0]->[0] = 2;

    for my $row (1 .. $rows) {
        for my $col (1 .. $cols) {
            my $gap1 = $matrix->[$row - 1]->[$col] + $gap_cost;
            my $gap2 = $matrix->[$row]->[$col - 1] + $gap_cost;
            my $match_mismatch =
                $matrix->[$row - 1]->[$col - 1] +
                $self->get_match_score(
                    $self->{seq1}->[$row - 1],
                    $self->{seq2}->[$col - 1]
                );

            my $max = $matrix->[$row]->[$col] =
                max($gap1, $gap2, $match_mismatch);

            $path_matrix->[$row]->[$col] = $max == $gap1
                    ? -1
                    : $max == $gap2
                    ? 1
                    : 0;
            }
        }
    }
于 2009-10-23T16:33:11.677 回答
7

而不是像这样取消引用您的二维数组：

$$path_matrix[0][0] = 2;

做这个：

$path_matrix->[0][0] = 2;

此外，你正在做很多 if/then/else 语句来匹配特定的子序列：这可以更好地写成given语句（perl5.10 相当于 C's switch）。在perldoc perlsyn阅读它：

given ($matrix->[$row][$column])
{
    when ($seq1_gap)       { $path_matrix->[$row][$column] = -1; }
    when ($match_mismatch) { $path_matrix->[$row][$column] = 0; }
    when ($seq2_gap)       { $path_matrix->[$row][$column] = 1; }
}
于 2009-10-23T16:36:25.147 回答
5

您的大部分代码都在处理二维数组。如果您想对数组做很多事情，我认为最大的改进将是改用PDL ，尤其是在考虑效率的情况下。它是一个提供出色数组支持的 Perl 模块。为了提高效率，底层例程在 C 中实现，因此速度也很快。

于 2009-10-23T16:34:16.573 回答
0

我总是建议查看CPAN以获取以前的解决方案或如何在 Perl 中做事的示例。你看过Algorithm::NeedlemanWunsch吗？

该模块的文档包括一个匹配 DNA 序列的示例。这是一个使用来自wikipedia的相似度矩阵的示例。

#!/usr/bin/perl -w
use strict;
use warnings;
use Inline::Files;                 #multiple virtual files inside code
use Algorithm::NeedlemanWunsch;    # refer CPAN - good style guide

# Read DNA sequences
my @a = read_DNA_seq("DNA_SEQ_A");
my @b = read_DNA_seq("DNA_SEQ_B");

# Read Similarity Matrix (held as a Hash of Hashes)
my %SM = read_Sim_Matrix();

# Define scoring based on "Similarity Matrix" %SM
sub score_sub {
    if ( !@_ ) {
        return -3;                 # gap penalty same as wikipedia)
    }
    return $SM{ $_[0] }{ $_[1] };    # Similarity Value matrix
}

my $matcher = Algorithm::NeedlemanWunsch->new( \&score_sub, -3 );
my $score = $matcher->align( \@a, \@b, { align => \&check_align, } );

print "\nThe maximum score is $score\n";

sub check_align {
    my ( $i, $j ) = @_;              # @a[i], @b[j]
    print "seqA pos: $i, seqB pos: $j\t base \'$a[$i]\'\n";
}

sub read_DNA_seq {
    my $source = shift;
    my @data;
    while (<$source>) {
        push @data, /[ACGT-]{1}/g;
    }
    return @data;
}

sub read_Sim_Matrix {

    #Read DNA similarity matrix (scores per Wikipedia)
    my ( @AoA, %HoH );
    while (<SIMILARITY_MATRIX>) {
        push @AoA, [/(\S+)+/g];
    }

    for ( my $row = 1 ; $row < 5 ; $row++ ) {
        for ( my $col = 1 ; $col < 5 ; $col++ ) {
            $HoH{ $AoA[0][$col] }{ $AoA[$row][0] } = $AoA[$row][$col];
        }
    }
    return %HoH;
}

__DNA_SEQ_A__
A T G T A G T G T A T A G T
A C A T G C A
__DNA_SEQ_B__
A T G T A G T A C A T G C A
__SIMILARITY_MATRIX__
-  A  G  C  T
A  10  -1  -3  -4
G  -1  7  -5  -3
C  -3  -5  9  0
T  -4  -3  0  8

这是一些示例输出：

seqA pos: 7, seqB pos: 2  base 'G'
seqA pos: 6, seqB pos: 1  base 'T'
seqA pos: 4, seqB pos: 0  base 'A'

The maximum score is 100
于 2009-10-25T11:32:24.240 回答