Как быстро найти общие предметы из двух массивов?

Question

Я пытаюсь найти общие строки среди двух разделенных табуляцией файлов на одном поле. Одна строка первого файла:

1  52854 s64199.1  A  .  .  .  PR  GT  0/0 

Одна строка второго файла:

chr1 52854  .  C  T  215.302 .  AB=0.692308;ABP=7.18621;AC=1;AF=0.5;AN=2;AO=9;CIGAR=1X;DP=13;DPB=13;DPRA=0;EPP=3.25157;EPPR=3.0103;GTI=0;LEN=1;MEANALT=1;MQM=60;MQMR=60;NS=1;NUMALT=1;ODDS=17.5429;PAIRED=0;PAIREDR=0.25;PAO=0;PQA=0;PQR=0;PRO=0;QA=318;QR=138;RO=4;RPP=3.25157;RPPR=5.18177;RUN=1;SAF=0;SAP=22.5536;SAR=9;SRF=1;SRP=5.18177;SRR=3;TYPE=snp;technology.illumina=1;BVAR GT:DP:RO:QR:AO:QA:GL 0/1:13:4:138:9:318:-5,0,-5 

Основываясь на втором поле (52854) в этом примере у меня есть много. Вот мой код, который находит общие, но мои файлы довольно большие и занимают много времени. Есть ли способ ускорить процесс? Спасибо вам большое заблаговременно.

#!/app/languages/perl/5.14.2/bin/perl 
use strict; 
use warnings; 
my $map_file = $ARGV[0]; 
my $vcf_file = $ARGV[1]; 
open my $map_info, $map_file or die "Could not open $map_file: $!"; 

my @map_array =(); 
my @vcf_array =(); 
while(my $mline = <$map_info>) { 
    chomp $mline; 
    my @data1 = split('\t', $mline); 
    my $pos1 = $data1[1]; 
    push (@map_array, $pos1); 
} 
open my $vcf_info, $vcf_file or die "Could not open $vcf_file: $!"; 
while(my $line = <$vcf_info>) { 
    if ($line !~ m/^#/) { 
      push (@vcf_array, $line); 
    } 
} 
foreach my $a (@map_array) { 
    chomp $a; 
foreach my $b (@vcf_array) { 
      chomp $b; 
      my @data = split('\t', $b); 
      my $pos2 = $data[1]; 
      my $ref2 = $data[3]; 
      my $allele = $data[4]; 
      my $genotype = $data[9]; 
      if ($a == $pos2) { 
       print $pos2 . "\t" . $ref2. "\t".$allele."\t".$genotype. "\n";  
      #print "$b\n"; 
      } 

    } 
} 

Answer 1

Ниже приведены минимальные изменения вашего сценария для хэш на основе поиска

use strict; 
use warnings; 
my $map_file = $ARGV[0]; 
my $vcf_file = $ARGV[1]; 

my %vcf_hash; 
open(my $vcf_info, $vcf_file) or die "Could not open $vcf_file: $!"; 
while(my $line = <$vcf_info>) { 
    next if $line =~ m/^#/; # Skip comment lines 
    chomp $line; 
    my (@data) = split(/\t/, $line); 
    die unless @data >= 10; # Check number of fields in the input line 
    my ($pos) = $data[1]; 
    # $. - line number in the file 
    $vcf_hash{$pos}{$.} = \@data; 
} 

open(my $map_info, $map_file) or die "Could not open $map_file: $!"; 
while(my $mline = <$map_info>) { 
    chomp $mline; 
    my (@data) = split(/\t/, $mline); 
    die unless @data >= 2; # Check number of fields in the input line 
    my ($pos) = $data[1]; 
    if(exists $vcf_hash{$pos}) { 
     my $hash_ref = $vcf_hash{$pos}; 
     for my $n (sort{$a<=>$b} keys %$hash_ref) { 
     my $array_ref = $hash_ref->{$n}; 
     my $pos2  = $array_ref->[1]; 
     my $ref2  = $array_ref->[3]; 
     my $allele = $array_ref->[4]; 
     my $genotype = $array_ref->[9]; 
     print $pos2 . "\t" . $ref2. "\t".$allele."\t".$genotype. "\n"; 
     } 
    } 
} 

Сценарий может быть улучшена дополнительно уменьшить использование памяти при использовании больших файлов данных.

Answer 2

Вот версия, которая должна работать намного быстрее, чем ваш собственный

Она считывает файл карты и сохраняет каждое pos поле в хэш %wanted. Затем он считывает второй файл и проверяет, находится ли запись в списке желаемых значений. Если да, то он разделяет записи и выводит поля, которые требуют

Обратите внимание, что я не был в состоянии проверить это за убедившись, что он компилирует

use strict; 
use warnings; 
use 5.010; 
use autodie; 

my ($map_file, $vcf_file) = @ARGV; 

my %wanted; 

{ 
    open my $map_fh, '<', $map_file; 

    while (<$map_fh>) { 
     chomp; 
     my $pos = (split /\t/, $_, 3)[1]; 
     ++$wanted{$pos}; 
    } 
} 

{ 
    open my $vcf_fh, '<', $vcf_file; 

    while (<$vcf_fh>) { 

     next if /^#/; 

     chomp; 
     my $pos = (split /\t/, $_, 3)[1]; 
     next unless $wanted{$pos}; 

     my ($ref, $allele, $genotype) = (split /\t/)[3, 4, 9]; 
     print join("\t", $pos, $ref, $allele, $genotype), "\n"; 

    } 
} 

Answer 3

Существует не нужно держать map_file в памяти, но только клавиши. Хорошо сделать их ключами в хеше, которые вы используете для проверки существования. Вам также не нужно хранить ваш vcf_file в памяти, но вы можете просто принять решение о выходе или нет.

#!/app/languages/perl/5.14.2/bin/perl 
use strict; 
use warnings; 
use autodie; 

use constant KEY => 1; 
use constant FIELDS => (1, 3, 4, 9); 

my ($map_file, $vcf_file) = @ARGV; 

my %map; 
{ 
    my $fh; 
    open $fh, '<', $map_file; 

    while (<$fh>) { 
     $map{ (split /\t/, $_, KEY + 2)[KEY] } = undef; 
    } 
} 

{ 
    my $fh; 
    open $fh, '<', $vcf_file; 
    while (<$fh>) { 
     next if /^#/; 
     chomp; 
     my @data = split /\t/; 
     print join "\t", @data[FIELDS] if exists $map{ $data[KEY] }; 
    } 
}