基于标识的两个表的比较和筛选

Question

我在编写一个脚本时遇到了一些小麻烦，该脚本根据身份遍历两个不同的表和筛选行。我刚刚意识到，到目前为止，这可能超出了我对perl的了解，所以我希望我能从你们那里得到一些方便的提示！

我有两个选项卡分隔的表，如下所示：

alleles.txt：

chr     pse.bp  bp      nalleles        maf     acc-1   acc-2   acc-3   acc-4    acc-5   acc-6    acc-7    acc-8   acc-9    acc-10    acc11    acc12    acc13 acc14    acc15
1       11      11      2       18      T       T       T       T       T       T       T       T       T       T       T       T       T       T       C       T
1       18      18      2       18      T       T       T       T       T       T       T       T       T       T       T       T       T       T       C       T
1       22      21.5    3       16      0       0       0       T       0       0       0       0       0       0       T       TCCTAAAT        0       0       0

hmp.txt：

rs#     alleles chrom   pos     strand  assembly#       center  protLSID        assayLSID       panelLSID       QCcode  M.10        M.101       M.103
NA      NA      1       10971981        NA      NA      NA      NA      NA      NA      NA      2       1       2
NA      NA      1       11716572        NA      NA      NA      NA      NA      NA      NA      1       1       1

我想用这样的方式编写一个脚本来处理这些数据：

我需要alleles.txt的信息，但我不需要每一行。我希望程序接受alleles.txt的第一行，如果它在hmp.txt中找到一行alleles.txt的第一列与hmp.txt的第三列相匹配，而alleles.txt的第三列与hmp.txt的第四列相匹配，则应该将这一行写入一个新文件中。对于allels.txt中的每一行，我都想这样做。我尝试用嵌套的while循环来完成这个任务，下面是我尝试运行的代码：

#!/usr/bin/perl
# maghap.pl
# converts pre-processed alleles layout into TASSEL-readable hapmap format.
# type ./maghap.pl hmp.txt alleles.txt to use this program.
use strict;#use warnings;

die "usage: ./maghap.pl RSB.lars.hmp.txt alleles.txt\n" unless (@ARGV == 2);
#open(my $hapmap, "<", "$ARGV[0]") or die "ERROR loading $ARGV[0]\n";
open(my $alleles, "<", "$ARGV[1]") or die "ERROR loading $ARGV[1]\n";
open(my $out, ">", "$ARGV[1].realsnps") or die "ERROR creating $ARGV[1].realsnps\n";

while (my $allelesline = <$alleles>) {
  #chomp;
  my @alleles_columns = split (/\t/, $allelesline);
  #print $out "@alleles_columns";
  #my $hit = 0;
  open(my $hapmap, "<", "$ARGV[0]") or die "ERROR loading $ARGV[0]\n";
  while (my $hapmapline = <$hapmap>) {
  #chomp;
    my @hapmap_columns = split(/\t/, $hapmapline);
    #print $out "@hapmap_columns";
    if ($alleles_columns[0] == $hapmap_columns[2]) {
      if ($alleles_columns[2] == $hapmap_columns[3]) {
        print $out "@alleles_columns";
        #print $out "@hapmap_columns";
        #$hit = 1;
        last;
      }
    }
  #print $out "@alleles_columns" if $hit;
  }
  close $hapmap;
}
#close $hapmap;
close $alleles;
close $out;

你可以从所有评论行中看到，我尝试了很多东西，但我现在似乎被困住了……到目前为止，该程序至少在运行，但由于某种原因，它找不到任何匹配(我检查过有匹配)。如果我关闭第二个If -条件(只在第一个if -条件中寻找匹配的东西)，它确实会找到很多匹配条件；但是，如果我关闭了第一个条件(只寻找匹配第二个条件的东西)，它就什么都找不到。我可能还应该提到，这两个文件都包含大约800.000行，所以为每一行遍历一个文件可能不是一个好的选择；但是，我想不出更好的选择了！我认为使用散列可能更容易，但我想不出一种聪明的方法将其编译成散列。如果你们对我有什么建议，我会很高兴的！

提前谢谢！

对虾

Answer 1

这适用于您提供的示例数据。以script.pl alleles.txt hmp.txt的形式运行

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

open my $AL,  '<', shift or die $!;
open my $HMP, '<', shift or die $!;

# Skip headers
<$AL>;
<$HMP>;

my ($chr_h, $pos_h) = (-1, -1);
while (<$AL>) {
    my ($chr_a, $pos_a) = (split /\t/)[0, 2];
    while ($chr_h < $chr_a and $pos_h < $pos_a) {
        ($chr_h, $pos_h) = (split /\t/, <$HMP>)[2, 3];
    }
    if ($chr_h == $chr_a and $pos_h == $pos_a) {
        print;
       ($chr_h, $pos_h) = (-1, -1);
    }
}

Answer 2

看来我必须把它放到一个答复中，因为我不能将代码放入注释=)

好的，我刚刚创建了这个迷你数据集，在其中我搜索了两行在两个文件中匹配的行，通过更少的方式：

alleles.txt：

rs#     alleles chrom   pos     strand  assembly#       center  protLSID        assayLSID       panelLSID       QCcode  M.10    M.101   M.103
NA      NA      1       10971981        NA      NA      NA      NA      NA      NA      NA      2       1       2
NA      NA      1       11716572        NA      NA      NA      NA      NA      NA      NA      1       1       1

hmp.txt：

chr     pse.bp  bp      nalleles        maf     bur-0   can-0   col-0   ct-1    edi-0   hi-0    kn-0    ler-0   mt-0    no-0    oy-0    po-0    rsch-4  sf-2    ts
1       11230382        10971981        3       14      GGTA    GGTA    GG      GG      GG      GG      GG      GG      GG      GG      GG      GGTA    GG      0
1       12050466        11716572        2       15      A       A       A       A       A       T       A       A       T       A       A       A       T       T

这两行在需要的字段中互相匹配，对吗？但是，如果我运行您给我的脚本，它只返回：

1       11230382        10971981        3       14      GGTA    GGTA    GG      GG      GG      GG      GG      GG      GG      GG      GG      GGTA    GG      0

这只是第一行。为了确保这不是原因，我应该提到，hmp.txt在末尾包含的列远远多于M.10、M.101、M.103。我只包含了前三列，因为实际的文件包含了大约1000列。顺便问一下，有没有办法也返回标题行？我知道我可以手动写，因为它不是那么多的文字，但只是为了它的地狱.

再次感谢你的帮助！

对虾