如何优化独特行的搜索？

Question

我有大的制表符分隔的文件，如下面的例子：

scaffold1443 182629 182998 chr1.1.1.241051.241420 367  99.80 
scaffold1443 131948 132412 chr1.1.2.291778.292242 462  99.80 
scaffold1443 96142 96474 chr1.1.3.327471.327803 330  99.70 
scaffold1443 53153 53479 chr1.1.4.370342.370668 322  99.40 
scaffold526  2870014 2870523 chr1.1.5.488372.488881 507  99.90 
scaffold526  2865956 2866314 chr1.1.6.490869.491234 357  98.10 
scaffold526  2867666 2868024 chr1.1.6.490869.491234 357  98.10 
scaffold526  2485557 2485867 chr1.1.7.610677.610987 310  100.00 

我想在一个新的文件只有行的第4列是唯一的打印。 在前面的示例中，除了第4列中包含“chr1.1.6.490869.491234”的两行外，应打印所有行。

我编写的以下脚本（它是较大管道的一部分）完成这项工作，但速度非常慢，尤其是当输入文件非常大时。

#!/usr/bin/perl 
use strict; 
use warnings; 
#This script takes the best hits output and finds the unique elements that up to only one scaffold. 

my $target = $ARGV[0]; 
my $chromosome = $ARGV[1]; 
my @mykeys = `cat OUTPUT_$target/psl_score_byname_$target/$chromosome.table| awk '{print \$4}'| sort -u`; 

foreach (@mykeys) 
{ 
    my $key = $_; 
    chomp($key); 
    my $command = "cat OUTPUT_$target/psl_score_byname_$target/$chromosome.table|grep -w $key"; 
    my @belongs= `$command`; 
    chomp(@belongs); 
    my $count = scalar(@belongs); 
    if ($count == 1) 
    { 
      open FILE, ">>OUTPUT_$target/unique_hces_$target/$chromosome.txt" or die $!; 
      print FILE "@belongs\n"; 

      @belongs =(); 
    } 
    else { 

      @belongs =(); 
    } 
} 

有没有更智能，更快捷的方法来做到这一点？ 非常感谢您提前。

Answer 1

鉴于您不想打印具有重复项目的行，您需要在打印之前查看整个文件，以首先找到具有重复项的行。然后返回并打印其他人。

这可以通过将整个文件与辅助数据结构一起保存在内存中，或者通过两次传递来完成。由于该文件是“很大”这里是更少的内存，使劲方式

use warnings; 
use strict; 

my $file = 'skip.txt'; 
open my $fh, '<', $file or die "Can't open $file: $!"; 

my (%seen, %dupe); 
while (<$fh>) 
{ 
    my $patt = (split)[3]; 

    # Record line numbers if the 4th field has been seen 
    if (exists $seen{$patt}) { 
     $dupe{ $seen{$patt} }++; # num of line with it seen first, with count 
     $dupe{$.} = 1;    # this line's number as well 
    } 
    else { $seen{$patt} = $. }  # first time this 4th field is seen 
} 

# Now we know all lines which carry duplicate fourth field 
my $outfile = 'filtered_' . $file; 
open my $fh_out, '>', $outfile or die "Can't open $outfile: $!"; 

seek $fh, 0, 0; # rewind to the beginning 
$. = 0;   # seek doesn't reset $. 
while (<$fh>) { 
    print $fh_out $_ if not exists $dupe{$.} 
} 
close $fh_out; 

重复发现其原线也需要被记录，$dupe{$seen{$patt}}++，在该分支的第一次。这个需求只需要完成一次，而我们可以检查（是否已经被记录），我们可能会选择一个可能有用的重复计数。

我已经添加了一些更多的副本（一些超过两倍）到您的发布示例，这会产生正确的输出。

---

评论上张贴的代码

张贴的代码检查关于对整个文件中的每一行中的第四字段中，从而处理该文件作为有线多次。这是很多工作，需要花时间，特别是对于大文件。

此外，没有理由为该作业使用外部程序。

Answer 2

由于oneliner：

perl -F"\t" -lanE 'push @l,[@F];$s{$F[3]}++}{say join"\t",@$_ for grep{$s{$_->[3]}==1}@l' <<EOF 
scaffold1443 182629 182998 chr1.1.1.241051.241420 367 99.80 
scaffold1443 131948 132412 chr1.1.2.291778.292242 462 99.80 
scaffold1443 96142 96474 chr1.1.3.327471.327803 330 99.70 
scaffold1443 53153 53479 chr1.1.4.370342.370668 322 99.40 
scaffold526 2870014 2870523 chr1.1.5.488372.488881 507 99.90 
scaffold526 2865956 2866314 chr1.1.6.490869.491234 357 98.10 
scaffold526 2867666 2868024 chr1.1.6.490869.491234 357 98.10 
scaffold526 2485557 2485867 chr1.1.7.610677.610987 310 100.00 
EOF 

输出

scaffold1443 182629 182998 chr1.1.1.241051.241420 367 99.80 
scaffold1443 131948 132412 chr1.1.2.291778.292242 462 99.80 
scaffold1443 96142 96474 chr1.1.3.327471.327803 330 99.70 
scaffold1443 53153 53479 chr1.1.4.370342.370668 322 99.40 
scaffold526 2870014 2870523 chr1.1.5.488372.488881 507 99.90 
scaffold526 2485557 2485867 chr1.1.7.610677.610987 310 100.00 

更具可读性：

perl -F"\t" -lanE ' 
    push @lines, [ @F ]; $seen{ $F[3] }++; 
    END { 
     say join("\t",@$_) for grep { $seen{ $_->[3] } == 1 } @lines 
    } 
' 

你可以，如果愿意，我创造了这个作为oneliner因为你说其翻译为完整的脚本：它是更大的管道的一部分。

另外需要注意的是，上面的内容首先将整个文件读入内存 - 所以非常大的文件可能会导致问题。

Answer 3

简单的方法包括使用关联数组来识别重复项。

perl -F'\t' -lane' 
    push @{ $h{ $F[3] } }, $_; 
    END { 
     for (values(%h)) { 
      print(@$_) if @$_ == 1; 
     } 
    } 
' file.tsv 

上述方法需要尽可能多的内存，因为文件很大。如果你的文件真的很大，那么这是一个不行。

---

如果你有真正的大文件，简单的办法是使用sort命令行实用程序对文件进行排序（这是相当快，能够处理任意大的文件）。通过首先重新排列文件以使重复文件彼此相邻，我们可以轻松过滤出重复文件，而不用担心内存问题。

sort -t$'\t' -k 4,4 file.tsv | perl -F'\t' -lane' 
    if ($key ne $F[3]) { 
     print(@buf) if @buf == 1; 
     @buf =(); 
    } 

    $key = $F[3]; 
    push @buf, $_; 

    END { print(@buf) if @buf == 1; } 
' 

---

如果你有真正的大文件，另一种比较简单的办法是加载在数据库（例如一个sqlite3的数据库）中的数据。您可以使用这种方法轻松维护原始订单。