perl - 循环遍历数组以使用 Perl 比较两个值

Question

我有一个大型数据集(29 列 x 19000 行)，我希望能够比较每一行的值并打印描述性输出。

具体来说，我想查询 A 列 (@WTcall) 中的值，它实际上是一个通过/失败语句。如果数据失败，我想打印“失败语句”并转到下一行，但如果数据通过，我想继续描述数据。

接下来的问题是判断X列(@positive)和Y列(@negative)的数据属于哪一类:

(例如:

如果 X 列和 Y 列 >= 0.6，则打印“ABC”

如果 X 列和 Y 列 < 0.6，则打印“CBA”

如果 X 列 >= 0.6 但 Y 列 < 0.6 打印“DEF”

如果 X 列 < 0.6 但 Y 列 >= 0.6 打印“FED”

否则打印“缺失数据”。 )

我已经包含了我在下面编写的代码以及示例数据的一个子集。

我在发布之前运行的测试在代码中被注释掉了。简而言之，如果我注释掉“if 和 elsif 语句”列表，打印“@WTcall\t@positive\t@negative\n”并通过 head 命令将其通过管道传输 - 我的变量似乎提取了正确的信息。

实际比较中出现了问题，因为每一行都使用“甲基化\t甲基化\n”描述进行分类。我不清楚这是为什么。我知道我有一些条目，其中 @WTcall 列应该匹配 $BadPosition(通过/失败检查)。此外，如果我再次注释掉“if 语句”，打印“@WTcall\n$BadPosition”并将其通过 sort 和 uniq 进行管道传输 - 我只得到一个“No_WT_Concensus”的值，因此应该没有拼写错误或匹配问题这些值。

我确信这个问题很明显，而且就在我面前，所以我非常感谢任何帮助。

谢谢。

代码:

#!/usr/bin/perl

use strict;
use warnings;
use autodie;
die "Usage: $0 Filename\n" if @ARGV != 1;

my $file = shift;
my @line;
my @positive; 
my @negative; 
my @WTcall;
my $BadPosition = 'No_WT_Concensus';
my $i;

open my $infh, '<', $file;

while (<$infh>) {
    chomp;
    @line = split(/\t/,$_);
    $WTcall[0]=$line[0];
    $positive[0]=$line[14];
    $negative[0]=$line[29];

    #foreach $1 (<$infh>) {
    foreach $1 (@WTcall) {
        if (@WTcall eq $BadPosition){
            print "No_WT_Concensus\tNo_WT_Concensus\n";
        } elsif (@positive >= 0.6 && @negative >= 0.6){
            print "Methylated\tMethylated\n";
        } elsif (@positive <= 0.6 && @negative <= 0.6){
            print "Under-methylated\tUnder-methylated\n";
        } elsif(@positive >= 0.6 && @negative <=0.6){
            print "Hemimethylated (m6A)\tHemimethylated (A)\n";
        } elsif(@positive <= 0.6 && @negative >= 0.6){
            print "Hemimethylated (A)\tHemimethylated (m6A)\n";
        } else{
            print "Missing_Site\tMissing_Site\n";
        }
        #print "@WTcall\n$BadPosition\n";
        #print "@WTcall\t@positive\t@negative\n"
        #print "@negative\n";
    }
}

close $infh;

示例数据:

Methylated  coding gene 619 thrA    NC_000913.3 pos 3   1   0.9535  1   NC_000913.3 619 +   18  0.8889  Methylated  coding gene 620 thrA    NC_000913.3 neg 3   0.9429  0.9756  0.9714  NC_000913.3 620 -   14  1
No_WT_Concensus coding gene 195410  ispU    NC_000913.3 pos 2   0.5789  0.766   0.6071  NC_000913.3 195410  +   39  0.5897  Methylated  coding gene 195411  ispU    NC_000913.3 neg 3   0.75    0.9074  0.9306  NC_000913.3 195411  -   21  0.8095
Under-methylated    pseudogene  3632965 yhiL    NC_000913.3 pos 0   0.0323  0.1429  0.0962  NC_000913.3 3632965 +   22  0.0909  Under-methylated    pseudogene  3632966 yhiL    NC_000913.3 neg 0   0.1463  0.175   0.1429  NC_000913.3 3632966 -   23  0
Methylated  intergenic  164636  hrpB-mrcB   NC_000913.3 pos 3   0.7381  0.7647  0.7273  NC_000913.3 164636  +   25  0.8 Methylated  intergenic  164637  hrpB-mrcB   NC_000913.3 neg 3   0.7 0.7931  0.7213  NC_000913.3 164637  -   25  0.4
Methylated  intergenic  269287  ykfA-perR   NC_000913.3 pos 3   0.875   0.8833  0.931   NC_000913.3 269287  +   22  0.8182  Methylated  intergenic  269288  ykfA-perR   NC_000913.3 neg 3   0.8077  0.6866  0.6491  NC_000913.3 269288  -   17  0.5294
Methylated  coding gene 4397856 mutL    NC_000913.3 pos 3   0.9245  0.9831  0.9661  blank   blank   blank   blank   blank   Methylated  coding gene 4397857 mutL    NC_000913.3 neg 3   0.9783  0.9808  0.9683  NC_000913.3 4397857 -   1   0
Methylated  coding gene 4397969 mutL    NC_000913.3 pos 3   0.9643  0.9524  1   blank   blank   blank   blank   blank   Methylated  coding gene 4397970 mutL    NC_000913.3 neg 3   1   1   1   blank   blank   blank   blank   blank
Methylated  coding gene 2761    thrA    NC_000913.3 pos 3   0.9259  0.8654  0.9242  NC_000913.3 2761    +   31  1   Methylated  coding gene 2762    thrA    NC_000913.3 neg 3   0.913   0.9636  0.9767  NC_000913.3 2762    -   29  0.9655
Methylated  coding gene 3073    thrB    NC_000913.3 pos 3   0.9677  0.8983  1   NC_000913.3 3073    +   29  1   Methylated  coding gene 3074    thrB    NC_000913.3 neg 3   1   0.9038  0.9778  NC_000913.3 3074    -   31  1

Answer 1

这是您的要求，经过编辑以显示并行结构:

• X >= 0.6 和 Y >= 0.6 然后是“ABC”

• X < 0.6 和 Y < 0.6 然后是“CBA”

• X >= 0.6 但 Y < 0.6 那么“DEF”

• X < 0.6 但 Y >= 0.6 那么“FED”

一些要求是< 0.6 , 但在你的代码中你有 <= 0.6 .

您有两件事要测试，您应该首先寻找一种结构，每件事只测试一次。这是表达这一点的伪代码:

if X >= 0.6
    if Y >= 0.6
        "ABC"
    else
        "DEF"
else
    if Y >= 0.6
        "FED"
    else
        "CBA"

一旦您知道某个值是否大于或等于某个值，您也知道它是否小于，因此您不必再次测试。测试只是 fork ；如果您不选择一个分支，则必须选择另一个分支。

您的代码有点不同，因为它同时匹配 $x >= 0.6或 $x <= 0.6 , $y 也一样.这意味着如果 $x和 $y是 0.6 ，那么任何分支都可以匹配，您将获得链中的第一个分支。这似乎是代码中的一个错误，因为它不是您在文本中描述的内容。

去掉结构

我不得不做很多项目，这些项目有很长的此类选择列表，其中许多项目有数百个要测试的东西。

现在，诀窍是将其转换为 Perl。请记住，子例程返回最后计算的表达式，所以这是可行的:

my @x_y = (
    [ 0.1, 0.7 ],
    [ 0.1, 0.1 ],
    [ 0.7, 0.1 ],
    [ 0.7, 0.7 ]
    );

foreach my $x_y ( @x_y ) {
    printf "X: %.1f Y: %.1f --> %s\n", @$x_y, get_value( @$x_y );
    }

sub get_value {
    my( $x, $y ) = @_;

    if( $x >= 0.6 ) { $y >= 0.6 ? 'ABC' : 'DEF' }
    else            { $y >= 0.6 ? 'FED' : 'CBA' }
    }

如果我不传递一个值，我什至可能会参数化枢轴值并给它一个值:

sub get_value {
    my( $x, $y, $pivot ) = @_;
    $pivot //= 0.6;  # default value

    if( $x >= $pivot ) { $y >= $pivot ? 'ABC' : 'DEF' }
    else               { $y >= $pivot ? 'FED' : 'CBA' }
    }

使用(实验性的)子例程签名，这会更简洁一些，因为我可以设置默认值:

use v5.22;
use feature qw(signatures);
no warnings qw(experimental::signatures);

sub get_value ( $x, $y, $pivot = 0.6 ){
    if( $x >= $pivot ) { $y >= $pivot ? 'ABC' : 'DEF' }
    else               { $y >= $pivot ? 'FED' : 'CBA' }
    }

但是事情会变得更有趣。我已经为 $y 重复了该测试，但我可以保存比较结果:

use v5.22;
use feature qw(signatures);
no warnings qw(experimental::signatures);

sub get_value ( $x, $y, $pivot = 0.6 ){
    my $boolean = ($y >= $pivot);
    if( $x >= $pivot ) { $boolean ? 'ABC' : 'DEF' }
    else               { $boolean ? 'FED' : 'CBA' }
    }

但是，我到底在做什么？我正在尝试选择一个值。我在代码中将其表示为决策树。如果我可以翻转它呢？我可以用 $x 做同样的保存结果技巧:

use v5.22;
use feature qw(signatures);
no warnings qw(experimental::signatures);

sub get_value ( $x, $y, $pivot = 0.6 ){
    my $y_boolean = ($y >= $pivot);
    my $x_boolean = ($x >= $pivot);

    if( $x_boolean ) { $y_boolean ? 'ABC' : 'DEF' }
    else             { $y_boolean ? 'FED' : 'CBA' }
    }

现在我遇到了 bool 组合 (0,0)、(0,1)、(1,0) 和 (1,1) 的情况。我将使用它们作为数组索引。 state生成一个持久变量，这样我就不需要在每次调用子例程时都重新定义它。 Perl v5.28 allows you to initialize arrays and hashes in state ，对于早期版本，您只需使用引用:

use v5.28;
use feature qw(signatures);
no warnings qw(experimental::signatures);

sub get_value ( $x, $y, $pivot = 0.6 ) {
    state @table = (
        [ qw(CBA FED) ],
        [ qw(DEF ABC) ]
        );

    my $y_boolean = ($y >= $pivot);
    my $x_boolean = ($x >= $pivot);

    $table[$x_boolean][$y_boolean];
    }

或者，稍微紧凑一些，我可以将比较放在数组索引中:

use v5.28;
use feature qw(signatures);
no warnings qw(experimental::signatures);

sub get_value ( $x, $y, $pivot = 0.6 ){
    state @table = (
        [ qw(CBA FED) ],
        [ qw(DEF ABC) ]
        );
    $table[$x >= $pivot][$y >= $pivot];
    }

现在值与机制分离了——我在 Mastering Perl 中花了很多时间.这也可以是一个参数，尽管现在它必须是一个数组引用，因为数组参数不能有默认值:

use v5.28;
use feature qw(signatures);
no warnings qw(experimental::signatures);

my @x_y = (
    [ 0.1, 0.7 ],
    [ 0.1, 0.1 ],
    [ 0.7, 0.1 ],
    [ 0.7, 0.7 ]
    );

my $pivot = 0.6;
my @table = ( [ qw(CBA FED) ], [ qw(DEF ABC) ] );

foreach my $x_y ( @x_y ) {
    printf "X: %.1f Y: %.1f --> %s\n", @$x_y,
        get_value( @$x_y, $pivot, \@table );
    }


sub get_value ( $x, $y, $pivot = 0.6,
    @table = ([ qw(DEF FED) ], [ qw(ABC CBA) ]) )
    {
    $table[$x >= $pivot][$y >= $pivot];
    }

这更容易管理。您可以在机制不变的情况下调整枢轴和返回的值。

更进一步，将值完全从程序中取出，并将它们放入配置文件中。然后，同一程序可以处理许多其他情况，而无需您对其进行编辑。

针对您的情况

回到你的代码。 choroba shows you this ，它解决了一些问题，但留下了 <=机智的问题:

while (<$infh>) {
    chomp;
    my @columns = split /\t/;
    my ($wt_call, $positive, $negative) = @columns[0, 14, 29];

    if ($wt_call eq $BadPosition) {
        print "No_WT_Concensus\tNo_WT_Concensus\n";
    } elsif ($positive >= 0.6 && $negative >= 0.6) {
        print "Methylated\tMethylated\n";
    } elsif ($positive <= 0.6 && $negative <= 0.6) {
        print "Under-methylated\tUnder-methylated\n";
    } elsif ($positive >= 0.6 && $negative <=0.6) {
        print "Hemimethylated (m6A)\tHemimethylated (A)\n";
    } elsif ($positive <= 0.6 && $negative >= 0.6) {
        print "Hemimethylated (A)\tHemimethylated (m6A)\n";
    } else {
        print "Missing_Site\tMissing_Site\n";
    }
}

清理了一下，你有一个while控制处理每一行的部分，以及一个处理控制值的部分的子程序。

use v5.28;
use feature qw(signatures);
no warnings qw(experimental::signatures);

while( <$infh> ) {
    chomp;
    my( $wt_call, $positive, $negative ) = (split /\t/)[0,14,29];
    if( $wt_call eq ... ) {
        print "No_WT_Concensus\tNo_WT_Concensus\n";
        next;
        }
    say get_value( $positive, $negative );
    }


sub get_value ( $x, $y, $pivot = 0.6 ){
    state @table = (
        [ qw(CBA FED) ],
        [ qw(DEF ABC) ]
        );
    $table[$x >= $pivot][$y >= $pivot];
    }

请注意 else condition 是无法达到的，因为你已经涵盖了所有的情况，所以我省略了。如果在另一种情况下你想查看空字段(null 与 0)，你会在 get_value 之前捕获它.一种方法是查看字段的长度。如果它是 0(没有字符)，则不要计算它。您可以有 0、1 或 2 个空字段，这些可能是不同的情况:

while( <$infh> ) {
    chomp;
    my( $wt_call, $positive, $negative ) = (split /\t/)[0,14,29];
    if( $wt_call eq ... ) {
        print "No_WT_Concensus\tNo_WT_Concensus\n";
        next;
        }
    # what if this is 1?
    unless( 2 == grep { length } ($positive, $negative) ) {
        print "No_WT_Concensus\tNo_WT_Concensus\n";
        next;
        }
    say get_value( $positive, $negative );
    }