【问题标题】:Removing lines in one file that are present in another file删除一个文件中存在于另一个文件中的行
【发布时间】:2013-08-14 18:13:33
【问题描述】:

我有 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;

 }

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++) {

?

【问题讨论】:

  • 我看到您要过滤掉的所有行在第 5 列都有一个 T。为什么不使用该模式来跳过这些行?
  • 这些是包含数千行的巨大文件,我只显示第一行的快照。除了 T 之外,还有其他字母要过滤掉。
  • 数千行并不是我认为的庞大。如果您有 1,000,000 行,每行 100 个字符,那仍然只有 100,000,000 个字节。即使有 100% 的开销,我们也只是说不到 200 MB。

标签: perl multidimensional-array nested-loops vcf-vcard


【解决方案1】:
#!/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

【讨论】:

  • 太好了,这段代码工作得很好。因为这些文件可能非常大 使用更少内存的最佳方法是什么?
  • 为什么这个嵌套循环不能正常工作?如果我想保留这种数组数组结构,为什么只在第一次进行迭代?将 foreach my $refH ( @ref_arrayH ) { 更改为 for (my $i = 0; $i
【解决方案2】:

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

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

#!/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 $!;

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 2011-06-14
    • 1970-01-01
    • 1970-01-01
    • 2021-09-18
    • 2012-10-07
    • 1970-01-01
    • 2015-02-04
    • 1970-01-01
    相关资源
    最近更新 更多