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