1

我有 2 个包含基因组数据的 .vcf 文件，我想删除第一个文件中也存在于第二个文件中的行。到目前为止，我的代码似乎只迭代一次，删除第一个命中，然后停止搜索。任何帮助将不胜感激，因为我无法弄清楚问题出在哪里。对不起，任何错误代码...

我选择了数组而不是哈希，因为必须维护文件的初始顺序，我认为使用数组可以更好地实现这一点......

代码：

#!/usr/bin/perl

use strict;
use warnings;

## bring files to program

MAIN: {

my ($chrC, $posC, $junkC, $baserefC, $allrestC);

my (@ref_arrayH, @ref_arrayC);

my ($chrH, $posH, $baserefH);
my $entriesH;

# open the het file first
open (HET, $het) or die "Could not open file $het - $!";
while (<HET>) {
    if (defined) {
        chomp;
        if (/^#/) { next; }
        if ( /(^.+?)\s(\d+?)\s(.+?)\s([A-Za-z\.]+?)\s([A-Za-z\.]+?)\s(.+?)\s(.+?)\s(.+)/m ) {   # is a VCF file
            my @line_arrayH = split(/\t/, $_);              
            push (@ref_arrayH, \@line_arrayH);      # the "reference" of each line_array is now in each element of @ref_array
            $entriesH = scalar @ref_arrayH;             # count the number of entries in the het file
        }
    }
}
close HET;

#   print $entriesH,"\n";

open (CNS, $cns) or die "Could not open file $cns - $!";

foreach my $refH ( @ref_arrayH ) {
    $chrH = $refH -> [0];
    $posH = $refH -> [1];
    $baserefH = $refH -> [3];

    foreach my $line (<CNS>) {
        chomp $line;
        if ($line =~ /^#/) { next; }
        if ($line =~ /(^.+?)\s(\d+?)\s(.+?)\s([A-Za-z\.]+?)\s([A-Za-z\.]+?)\s(.+?)\s(.+?)\s(.+)/m ) {   # is a VCF file
            ($chrC, $posC, $junkC, $baserefC, $allrestC) = split(/\t/,$line);
                if ( $chrC eq $chrH and $posC == $posH and $baserefC eq $baserefH ) { next }
                else { print "$line\n"; }
        }
    }
}
 #  close CNS;

 }

中枢神经系统文件：

chrI    1084    .   A   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1085    .   C   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1086    .   A   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1087    .   C   T   3.55    .   DP=1;AF1=1;AC1=2;DP4=0,0,1,0;MQ=52;FQ=-30   GT:PL:GQ    0/1:31,3,0:4
chrI    1088    .   T   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1089    .   A   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1090    .   C   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1091    .   T   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1099    .   A   .   32.8    .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30.2 PL  0
chrI    1100    .   G   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1101    .   G   .   12.3    .   DP=1;AF1=1;AC1=2;DP4=0,0,1,0;MQ=52;FQ=-30.1 PL  18
chrI    1102    .   G   .   32.9    .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30.1 PL  0
chrI    1103    .   A   .   5.45    .   DP=1;AF1=1;AC1=2;DP4=0,0,1,0;MQ=52;FQ=-30   PL  26
chrI    1104    .   C   T   3.55    .   DP=1;AF1=1;AC1=2;DP4=0,0,1,0;MQ=52;FQ=-30   GT:PL:GQ    0/1:31,3,0:4
chrI    1105    .   T   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0

HET 文件：

chrI    1087    .   C   T   3.55    .   DP=1;AF1=1;AC1=2;DP4=0,0,1,0;MQ=52;FQ=-30   GT:PL:GQ    0/1:31,3,0:4
chrI    1104    .   C   T   3.55    .   DP=1;AF1=1;AC1=2;DP4=0,0,1,0;MQ=52;FQ=-30   GT:PL:GQ    0/1:31,3,0:4

我希望输出是这样的

chrI    1084    .   A   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1085    .   C   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1086    .   A   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1088    .   T   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1089    .   A   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1090    .   C   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1091    .   T   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1099    .   A   .   32.8    .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30.2 PL  0
chrI    1100    .   G   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1101    .   G   .   12.3    .   DP=1;AF1=1;AC1=2;DP4=0,0,1,0;MQ=52;FQ=-30.1 PL  18
chrI    1102    .   G   .   32.9    .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30.1 PL  0
chrI    1103    .   A   .   5.45    .   DP=1;AF1=1;AC1=2;DP4=0,0,1,0;MQ=52;FQ=-30   PL  26
chrI    1105    .   T   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0

但给了我这个：

chrI    1084    .   A   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1085    .   C   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1086    .   A   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1088    .   T   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1089    .   A   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1090    .   C   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1091    .   T   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1099    .   A   .   32.8    .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30.2 PL  0
chrI    1100    .   G   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0
chrI    1101    .   G   .   12.3    .   DP=1;AF1=1;AC1=2;DP4=0,0,1,0;MQ=52;FQ=-30.1 PL  18
chrI    1102    .   G   .   32.9    .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30.1 PL  0
chrI    1103    .   A   .   5.45    .   DP=1;AF1=1;AC1=2;DP4=0,0,1,0;MQ=52;FQ=-30   PL  26
chrI    1104    .   C   T   3.55    .   DP=1;AF1=1;AC1=2;DP4=0,0,1,0;MQ=52;FQ=-30   GT:PL:GQ    0/1:31,3,0:4
chrI    1105    .   T   .   33  .   DP=1;AF1=0;AC1=0;DP4=1,0,0,0;MQ=52;FQ=-30   PL  0

为什么这个嵌套循环不能正常工作？如果我想保留这种数组数组结构，为什么只在第一次进行迭代？

改变foreach循环会更好吗

foreach my $refH ( @ref_arrayH ) {

带有 for 循环

for (my $i = 0; $i <= $entriesH; $i++) {

?

4

2 回答 2

4 
#!/usr/bin/env perl

use strict;
use warnings;

my %seen;

open my $het, '<', 't.het' or die $!;
$seen{ $_ } = undef while <$het>;
close $het or die $!;

open my $cns, '<', 't.cns' or die $!;

while (my $line = <$cns>) {
    next if exists $seen{ $line };
    print $line;
}

close $cns or die $!;

如果您想匹配整行以外的内容，请提取字段并使用它（或它们的组合）键入%seen散列。

这将使用与HET文件中的行数成比例的内存。

为了减少内存使用，您可以绑定%seen到磁盘上的DBM_File。

于 2013-08-14T18:49:07.267 回答
0

如果您担心内存使用情况，您可以在进行比较时一次读取两个文件一行。以下是另一种方法。

注意：由于文件句柄的工作方式，我们每次在第二个嵌套循环中读取文件时都必须重置连接。

#!/usr/bin/env perl

use strict;
use warnings;

open my $cns, '<', 't.cns' or die $!;

CNS: 
while (my $line = <$cns>) { #read one line at a time from t.cns file.
    open my $het, '<', 't.het' or die $!;
    while (my $reference = <$het>){
            if ($line eq $reference) { #test if current t.cns line is equal to any line in t.hex file. 
                close $het; #reset conection to t.hex file. 
                next CNS; # skip current t.cns line.            
    }
}   
    print $line;  
    close $het; #reset conection to t.hex file.
}
close $cns or die $!;
于 2013-08-15T13:36:10.397 回答