Perl/Linux, фильтрующий большой файл с содержимым другого файла

Question

Я фильтрую файл размером 580 МБ, используя содержимое другого меньшего файла. File1 (меньший размер файла)

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

File2 (большой файл)

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

Прочитайте небольшой файл в структуру данных и проверять каждую строку другого файла против него.

Здесь я прочитал его в массиве, каждый из которых является массивом с полями из строки. Затем каждая строка файла данных проверяется против arrayrefs в этом массиве, сравнивая поля для каждого требования.

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; 

Это выводит правильные строки из предоставленных образцов. Это позволяет использовать несколько строк. Если этого как-то не может быть, добавьте last в блок if, чтобы выйти из foreach, как только совпадение найдено.

Несколько комментариев к коду. Позвольте мне знать, может ли быть полезным больше.

При чтении ссылочного файла <$fh> используется в контексте списка, чтобы он возвращал все строки, а grep отфильтровывает пустые. map первый символ chomp с символом новой строки, а затем делает arrayref на [ ], причем элементы являются полями на линии, полученными split. Список выходных данных присваивается @ref.

При повторном использовании $fh он закрыт первым (если он был открыт), поэтому нет необходимости в close.

Я храню строки заголовков просто так, возможно, для печати или проверки. Нам действительно нужно их исключить.

Answer 2

Как говорилось ранее, на каждой итерации происходит очень медленное вызов awk. Полное awk решения было бы возможно, я только что видел решение Perl, вот мое решение Python как OP не будет возражать:

• создать словарь из маленького файла: список CHR => пара начала/конца
• итерации через большой файл и попытайтесь сопоставить 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 

Вы можете сделать это в Perl через DBI и драйвер DBD::SQLite.

Answer 5

awk power с одним проходом. Ваш код выполняет итерацию файла2 столько раз, сколько строк в файле1, поэтому время выполнения линейно увеличивается. Пожалуйста, дайте мне знать, если это однопроходное решение работает медленнее, чем другие решения.

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