Perl/Linux过滤大文件与其他文件的内容

Question

我正在使用另一个较小文件的内容过滤580 MB文件。 File1中（较小的文件）

chr start End 
1 123 150 
2 245 320 
2 450 600 

文件2（大文件）

chr pos RS ID A B C D E F 
1 124 r2 3 s 4 s 2 s 2 
1 165 r6 4 t 2 k 1 r 2 
2 455 t2 4 2 4 t 3 w 3 
3 234 r4 2 5 w 4 t 2 4 

我想如果下列条件满足，以捕获来自文件2行。 File2.Chr == File1.Chr && File2.Pos > File1.Start && File2.Pos < File1.End 我试过使用awk，但它运行速度非常慢，我也想知道是否有更好的方法来实现相同？

谢谢。

这里是我正在使用的代码：

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

my $bed_file = "/data/1000G/Hotspots.bed";#File1 smaller file 
my $SNP_file = "/data/1000G/SNP_file.txt";#File2 larger file 
my $final_file = "/data/1000G/final_file.txt"; #final output file 

open my $in_fh, '<', $bed_file 
     or die qq{Unable to open "$bed_file" for input: $!}; 

    while (<$in_fh>) { 

    my $line_str = $_; 

    my @data = split(/\t/, $line_str); 

    next if /\b(?:track)\b/;# skip header line 
    my $chr = $data[0]; $chr =~ s/chr//g; print "chr is $chr\n"; 
    my $start = $data[1]-1; print "start is $start\n"; 
    my $end = $data[2]+1; print "end is $end\n"; 

    my $cmd1 = "awk '{if(\$1==chr && \$2>$start && \$2</$end) print (\"chr\"\$1\"_\"\$2\"_\"\$3\"_\"\$4\"_\"\$5\"_\"\$6\"_\"\$7\"_\"\$8)}' $SNP_file >> $final_file"; print "cmd1\n"; 
    my $cmd2 = `awk '{if(\$1==chr && \$2>$start && \$2</$end) print (\"chr\"\$1\"_\"\$2\"_\"\$3\"_\"\$4\"_\"\$5\"_\"\$6\"_\"\$7\"_\"\$8)}' $SNP_file >> $final_file`; print "cmd2\n"; 

} 

Answer 1

阅读小文件到数据结构和对证其他文件的每一行。

这里我将它读入一个数组中，每个元素都是一个arrayref字段中的字段。然后根据数组中的每个行检查数组文件，比较每个需求的字段。

use warnings 'all'; 
use strict; 

my $ref_file = 'reference.txt'; 
open my $fh, '<', $ref_file or die "Can't open $ref_file: $!"; 
my @ref = map { chomp; [ split ] } grep { /\S/ } <$fh>; 

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

# Drop header lines 
my $ref_header = shift @ref;  
my $data_header = <$fh>; 

while (<$fh>) 
{ 
    next if not /\S/; # skip empty lines 
    my @line = split; 

    foreach my $refline (@ref) 
    { 
     next if $line[0] != $refline->[0]; 
     if ($line[1] > $refline->[1] and $line[1] < $refline->[2]) { 
      print "@line\n"; 
     } 
    } 
} 
close $fh; 

这将从所提供的样本中打印出正确的线条。它允许多行匹配。如果不知何故，请在if区块中添加last，以便在找到匹配项后退出foreach。

关于代码的几点意见。让我知道更多是否有用。

读取参考文件时，在列表上下文中使用<$fh>，因此它将返回所有行，并且grep将过滤掉空白行。 map第一个chomp是换行符，然后通过[ ]生成一个arrayref，元素是由split获得的行上的字段。输出列表分配给@ref。

当我们重用$fh它首先关闭（如果它是开放的），所以不需要close。

我只是这样存储标题行，可能是为了打印或检查。我们真的只需要排除它们。

Answer 2

如前所述，在每次迭代中调用awk是非常缓慢的。全awk解决方案是可行的，我刚才看到一个Perl的解决方案，这是我的Python的溶液，作为OP不会介意：

• 创建小文件的字典：夫妻CHR =>列表开始/结束
• 遍历大文件并尝试匹配其中一个开始/结束元组之间的位置。

代码：

with open("smallfile.txt") as f: 
    next(f) # skip title 
    # build a dictionary with chr as key, and list of start,end as values 
    d = collections.defaultdict(list) 
    for line in f: 
     toks = line.split() 
     if len(toks)==3: 
      d[int(toks[0])].append((int(toks[1]),int(toks[2]))) 


with open("largefile.txt") as f: 
    next(f) # skip title 
    for line in f: 
     toks = line.split() 
     chr_tok = int(toks[0]) 
     if chr_tok in d: 
      # key is in dictionary 
      pos = int(toks[1]) 
      if any(lambda x : t[0]<pos<t[1] for t in d[chr_tok]): 
       print(line.strip()) 

我们可以通过排序元组的列表，并appyling bisect避免线性搜索速度稍快。只有在“小”文件中元组列表很大时才需要这样做。

Answer 3

的另一种方式，这一次存储阵列（HOA）的哈希基于“CHR”字段中较小的文件：

use strict; 
use warnings; 

my $small_file = 'small.txt'; 
my $large_file = 'large.txt'; 

open my $small_fh, '<', $small_file or die $!; 

my %small; 

while (<$small_fh>){ 
    next if $. == 1; 
    my ($chr, $start, $end) = split /\s+/, $_; 
    push @{ $small{$chr} }, [$start, $end]; 
} 

close $small_fh; 

open my $large_fh, '<', $large_file or die $!; 

while (my $line = <$large_fh>){ 
    my ($chr, $pos) = (split /\s+/, $line)[0, 1]; 

    if (defined $small{$chr}){ 
     for (@{ $small{$chr} }){ 
      if ($pos > $_->[0] && $pos < $_->[1]){ 
       print $line; 
      } 
     } 
    } 
} 

Answer 4

把它们放入一个SQLite数据库，做加盟。这会比尝试自己写一些东西快得多，少用多少内存，使用的内存也少。而且它更加灵活，现在您只需对数据执行SQL查询即可，无需继续编写新脚本并重新分析文件。

您可以通过解析和插入自己来导入它们，也可以将它们转换为CSV并使用SQLite's CSV import ability。使用这些简单的数据转换为CSV可以像s{ +}{,}g一样简单，或者您可以使用全面且非常快速的Text::CSV_XS。

你的表看起来像这样（你会希望使用更好的表名和字段名称）。

create table file1 (
    chr integer not null, 
    start integer not null, 
    end integer not null 
); 

create table file2 (
    chr integer not null, 
    pos integer not null, 
    rs integer not null, 
    id integer not null, 
    a char not null, 
    b char not null, 
    c char not null, 
    d char not null, 
    e char not null, 
    f char not null 
); 

在您要搜索的列上创建一些索引。索引会减慢导入速度，所以请确保在导入后执行。

create index chr_file1 on file1 (chr); 
create index chr_file2 on file2 (chr); 
create index pos_file2 on file2 (pos); 
create index start_file1 on file1 (start); 
create index end_file1 on file1 (end); 

然后进行连接。

select * 
from file2 
join file1 on file1.chr == file2.chr 
where file2.pos between file1.start and file1.end; 

1,124,r2,3,s,4,s,2,s,2,1,123,150 
2,455,t2,4,2,4,t,3,w,3,2,450,600 

您可以通过DBI在Perl做到这一点和DBD::SQLite驱动程序。

Answer 5

awk power with single pass。你的代码迭代file1的次数与file1中的行数一样多，所以执行时间线性增加。请让我知道这个单通解决方案是否比其他解决方案慢。

awk 'NR==FNR { 
    i = b[$1];  # get the next index for the chr 
    a[$1][i][0] = $2; # store start 
    a[$1][i][1] = $3; # store end 
    b[$1]++;   # increment the next index 
    next; 
} 

{ 
    p = 0; 
    if ($1 in a) { 
     for (i in a[$1]) { 
      if ($2 > a[$1][i][0] && \ 
       $2 < a[$1][i][1]) 
       p = 1     # set p if $2 in range 
     } 
    } 
} 

p {print}' 

一个班轮

awk 'NR==FNR {i = b[$1];a[$1][i][0] = $2; a[$1][i][1] = $3; b[$1]++;next; }{p = 0;if ($1 in a){for(i in a[$1]){if($2>a[$1][i][0] && $2<a[$1][i][1])p=1}}}p' file1 file2