Как найти числовой диапазон, в котором находится данная цифра?

Question

У меня есть два файла

file1

Данные SNP, состоящие из хромосомы и ее положения (приблизительно 400 000 записей)

chr pos
a1 456
a2 789
 . .
 . . 
so on

file2

Данные GTF, состоящие из хромосомы, position_start, position_end и подробностей (приблизительно 500 000 записей)

chr pos_start pos_end detail
a1 100 400 gene1
a1 401 700 gene2
a2 200 500 gene3
a2 501 900 gene4
 . .
 . . 
so on

Желаемый результат

chr pos chr pos_start pos_end detail
a1 456 a1 401 700 gene2
a2 789 a2 501 900 gene4

Я получаю этот результат с помощью сценария оболочки:

(grep "$chr" file2.gtf | awk '{if($2 <= '$pos' && $3 >= '$pos') print $0}') 

в цикле while, но обработка всех цифр в file1 занимает слишком много времени.

Кто-нибудь знает более эффективный способ в shell, Python или Perl для этого?

Answer 1

Вот версия Perl.Основная идея заключается в том, что он кэширует данные gtf в хеш-таблицу, а затем для каждой строки в файле snp он просматривает только записи gtf, соответствующие этой хромосоме.

#!/usr/bin/perl
use warnings;
use strict;
use feature qw/say/;
use autodie;

my $snp_file = "file1.txt";
my $gtf_file = "file2.txt";

# Read the gtf data into a hash of arrays
my %gtf;
open my $file, "<", $gtf_file;
my $hdr = <$file>; # Discard header line
while (<$file>) {
  chomp;
  my @cols = split /\s+/;
  push @{$gtf{$cols[0]}}, \@cols;
}
close $file;

open $file, "<", $snp_file;
$hdr = <$file>; # Discard header line
say "chr\tpos\tchr\tstart\tend\tdetail";
# Read the snp data
$" = "\t"; # Use tab for array element separator
while (<$file>) {
  chomp;
  my ($chr, $pos) = split /\s+/;
  # Look up all matches of this chromosome in the gtf hash and filter just
  # the ones where pos is in range.
  my @matches = grep { $pos >= $_->[1] && $pos <= $_->[2] } @{$gtf{$chr}};
  # And print them out.
  for my $match (@matches) {
    say "$chr\t$pos\t@$match";
  }
}
close $file;

Другой вариант, который я хотел быпродолжайте, если вы собираетесь много работать с этими данными, загрузите их все в sqlite или другую базу данных и ищите результаты с помощью SQL.Таким образом, вам не нужно продолжать читать файлы данных;вы просто просматриваете вещи в предварительно заполненной таблице (с соответствующими индексами, чтобы сделать вещи эффективными).

Answer 2

Я думаю, что это делает то, что вы хотите с awk:

awk '
   FNR==1  { next}
   FNR==NR { chr[FNR]=$1; start[FNR]=$2; end[FNR]=$3; det[FNR]=$4; N=FNR; next}
           { c=$1; p=$2;
             for(i=2;i<=N;i++){
                if((c==chr[i]) && (p>=start[i]) && (p<=end[i])){
                   print c, p, chr[i], start[i], end[i], det[i]
                   next
                }
             }
           }
   ' file2 file1

Итак, сначала обратите внимание на последнюю строку, что один вызов awk обрабатывает оба файла.

Внутри обработки первая строка каждого файла игнорируется, проверяя, равен ли номер строки в текущем файле 1, и пропуская, если так:

FNR==1  { next}

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

FNR==NR { chr[FNR]=$1; start[FNR]=$2; end[FNR]=$3; det[FNR]=$4; N=FNR; next}

В противном случае мы должны обрабатывать второй файл.В этом случае мы перебираем все массивы, которые мы сохранили из первого файла, чтобы найти подходящую запись.Если мы находим один в правильном диапазоне, мы печатаем нужные вам биты и кусочки и сразу переходим к следующей записи:

{ c=$1; p=$2;
  for(i=2;i<=N;i++){
     if((c==chr[i]) && (p>=start[i]) && (p<=end[i])){
         print c, p, chr[i], start[i], end[i], det[i]
         next
     }
   }
}