Загрузка больших файлов в хэши в Perl (таблицы BLAST)

Я новичок в Perl, пожалуйста, помогите мне с моим запросом... Я пытаюсь извлечь информацию из таблицы blast (фрагмент того, как это выглядит ниже): Это стандартный ввод таблицы blast... I в основном хочу извлечь любую информацию из списка прочтений (посмотрите на мой второй скрипт ниже, чтобы получить представление о том, что я хочу сделать).... Во всяком случае, это именно то, что я сделал во втором скрипте:

ВХОДЫ:

1) взрывная таблица:

38.1    0.53    59544   GH8NFLV01A02ED  GH8NFLV01A02ED rank=0113471 x=305.0 y=211.5 length=345  1   YP_003242370    Dynamin family protein [Paenibacillus sp. Y412MC10] -1  0   48.936170212766 40.4255319148936    47  345 1213    13.6231884057971    3.87469084913438    31  171 544 590
34.3    7.5 123828  GH8NFLV01A03QJ  GH8NFLV01A03QJ rank=0239249 x=305.0 y=1945.5 length=452 1   XP_002639994    Hypothetical protein CBG10824 [Caenorhabditis briggsae] 3   0   52.1739130434783    32.6086956521739    46  452 367 10.1769911504425    12.5340599455041    111 248 79  124
37.7    0.70    62716   GH8NFLV01A09B8  GH8NFLV01A09B8 rank=0119267 x=307.0 y=1014.0 length=512 1   XP_002756773    PREDICTED: probable G-protein coupled receptor 123-like, partial [Callithrix jacchus]   1   0   73.5294117647059    52.9411764705882    34  512 703 6.640625    4.83641536273115    43  144 273 306
37.7    0.98    33114   GH8NFLV01A0H5C  GH8NFLV01A0H5C rank=0066011 x=298.0 y=2638.5 length=573 1   XP_002756773    PREDICTED: probable G-protein coupled receptor 123-like, partial [Callithrix jacchus]   -3  0   73.5294117647059    52.9411764705882    34  573 703 5.93368237347295    4.83641536273115    131 232 273 306
103 1e-020  65742   GH8NFLV01A0MXI  GH8NFLV01A0MXI rank=0124865 x=300.5 y=644.0 length=475  1   ABZ08973    hypothetical protein ALOHA_HF4000APKG6B14ctg1g18 [uncultured marine crenarchaeote HF4000_APKG6B14]  2   0   77.9411764705882    77.9411764705882    68  475 151 14.3157894736842    45.0331125827815    2   205 1   68
41.6    0.053   36083   GH8NFLV01A0QKX  GH8NFLV01A0QKX rank=0071366 x=301.0 y=1279.0 length=526 1   XP_766153   hypothetical protein [Theileria parva strain Muguga]    -1  0   66.6666666666667    56.6666666666667    30  526 304 5.70342205323194    9.86842105263158    392 481 31  60
45.4    0.003   78246   GH8NFLV01A0Z29  GH8NFLV01A0Z29 rank=0148293 x=304.0 y=1315.0 length=432 1   ZP_04111769 hypothetical protein bthur0007_56280 [Bacillus thuringiensis serovar monterrey BGSC 4AJ1]   3   0   51.8518518518518    38.8888888888889    54  432 193 12.5    27.979274611399 48  209 97  150
71.6    4e-011  97250   GH8NFLV01A14MR  GH8NFLV01A14MR rank=0184885 x=317.5 y=609.5 length=314  1   ZP_03823721 DNA replication protein [Acinetobacter sp. ATCC 27244]  1   0   92.5    92.5    40  314 311 12.7388535031847    12.8617363344051    193 312 13  52
58.2    5e-007  154555  GH8NFLV01A1KCH  GH8NFLV01A1KCH rank=0309994 x=310.0 y=2991.0 length=267 1   ZP_03823721 DNA replication protein [Acinetobacter sp. ATCC 27244]  1   0   82.051282051282 82.051282051282 39  267 311 14.6067415730337    12.540192926045 142 258 1   39

2) Список прочтений:

GH8NFLV01A09B8
GH8NFLV01A02ED
etc
etc

3) результат, который я хочу:

37.7    0.70    62716   GH8NFLV01A09B8  GH8NFLV01A09B8 rank=0119267 x=307.0 y=1014.0 length=512 1   XP_002756773    PREDICTED: probable G-protein coupled receptor 123-like, partial [Callithrix jacchus]   1   0   73.5294117647059    52.9411764705882    34  512 703 6.640625    4.83641536273115    43  144 273 306
38.1    0.53    59544   GH8NFLV01A02ED  GH8NFLV01A02ED rank=0113471 x=305.0 y=211.5 length=345  1   YP_003242370    Dynamin family protein [Paenibacillus sp. Y412MC10] -1  0   48.936170212766 40.4255319148936    47  345 1213    13.6231884057971    3.87469084913438    31  171 544 590

Мне нужно подмножество информации в первом списке, учитывая список прочитанных имен, которые я хочу извлечь (который находится в 4-м столбце). Вместо хэширования списка прочтений (только?) я хочу хэшировать саму таблицу blast, и используйте информацию в столбце 4 (таблицы blast) в качестве ключей для извлечения значений каждого ключа, даже если этот ключ может иметь более одного значения (т. е. каждое прочитанное имя может фактически иметь более одного попадания или связанного blast результат в таблице), имея в виду, что значение включает ВСЮ строку с этим ключом (имя чтения) в нем.

Мой скрипт greplist.pl делает это, но я думаю, что он очень-очень медленный (и поправьте меня, если я ошибаюсь), что, загружая всю таблицу в хэш, это должно значительно ускорить работу...

Спасибо за помощь.

Мои скрипты: Сломанный (mambo5.pl)

#!/usr/bin/perl -w
# purpose:  extract blastX data from a list of readnames
use strict;
open (DATA,$ARGV[0]) or die ("Usage: ./mambo5.pl BlastXTable readslist");
open (LIST,$ARGV[1]) or die ("Usage: ./mambo5.pl BlastXTable readslist");
my %hash = <DATA>;
close (DATA);
my $filename=$ARGV[0];
open(OUT, "> $filename.bololom");

my $readName;

while ( <LIST> )
{
    #########;
    if(/^(.*?)$/)#
    {
        $readName=$1;#
        chomp $readName;
        if (exists $hash{$readName})
        {
            print "bingo!";
            my $output =$hash{$readName};
            print OUT "$output\n";
        }
        else 
        {
            print "it aint workin\n";
            #print %hash;
        }           
    }
}
close (LIST);

Медленный и быстрый чит (который работает) и очень медленный (мои бласт-таблицы могут иметь размер от 400 МБ до 2 ГБ, я уверен, вы понимаете, почему он такой медленный)

#!/usr/bin/perl -w
## 
# This script finds a list of names in a blast table and outputs the result in a new file
# name must exist and list must be correctly formatted
# will not output anything using a "normal" blast file, must be a table blast
# if you have the standard blast output use blast2table script

use strict;
my $filein=$ARGV[0] or die ("usage: ./listgrep.pl readslist blast_table\n");
my $db=$ARGV[1] or die ("usage: ./listgrep.pl readslist blast_table\n");
#open the reads you want to grep
my $read;
my $line;
open(READSLIST,$filein);
while($line=<READSLIST>)
{
    if ($line=~/^(.*)$/) 
    {
        $read = $1;
        print "$read\n";
        system("grep \"$read\" $db >$read\_.out\n");
    }


    #system("grep $read $db >$read\_.out\n");
}
system("cat *\_.out >$filein\_greps.txt\n");
system("rm *.out\n");

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

Спасибо !


split hash perl blast
person Gimly_Gloin    schedule 09.03.2012    source источник
comment
Парсинг выглядит как работа для Text::CSV.   -  person matthias krull    schedule 09.03.2012
comment
Ваша строка my %hash = <DATA> не делает ничего полезного, хотя было бы неплохо, если бы все работало именно так! Вы явно понятия не имеете, что такое хэши Perl, и предлагаете взглянуть на Ускоренный курс Hash на сайте www.perl.com.   -  person Borodin    schedule 09.03.2012

Ответы (3)


arrow_upward
4
arrow_downward

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

#!/usr/bin/perl 
use strict;
use warnings;
use 5.010;

# Read the readslist file into a hash
open my $fh, '<', 'readslist' or die "Can't open 'readslist' for reading:$!";
my %readslist = map { chomp; $_ => 1 }<$fh>;
close $fh;

open my $fh_blast, '<', 'blastfile' or die "Can't open 'blastfile' for reading:$!";
# loop on all the blastfile lines
while (<$fh_blast>) {
    chomp;
    # retrieve the key (4th column)
    my ($key) = (split/\s+/)[3];
    # print the line if the key exists in the hash
    say $_ if exists $readslist{$key};
}
close $fh_blast;
person Toto    schedule 09.03.2012

arrow_upward
0
arrow_downward

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

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

use strict;
use warnings;

use Fcntl qw/SEEK_SET/;

my %index;

open my $blast, '<', 'blast.txt' or die $!;

until (eof $blast) {
  my $place = tell $blast;
  my $line = <$blast>;
  my $key = (split ' ', $line, 5)[3];
  push @{$index{$key}}, $place;
}

open my $reads, '<', 'reads.txt' or die $!;

while (<$reads>) {

  next unless my ($key) = /(\S+)/;
  next unless my $places = $index{$key};

  foreach my $place (@$places) {
    seek $blast, $place, SEEK_SET;
    my $line = <$blast>;
    print $line;
  }
}
person Borodin    schedule 09.03.2012
comment
да, я думаю, это работает, и да, я понятия не имею, как работают хэши ... спасибо за ваш вклад. - person Gimly_Gloin; 10.03.2012
comment
Было бы интересно узнать, насколько это быстрее, чем ваше существующее решение, и сколько времени требуется для индексации файла размером 2 ГБ (только от первого цикла до второго открытия). Спасибо за интересный вопрос. - person Borodin; 11.03.2012
comment
они оба довольно быстрые, я не могу сказать разницу... еще раз спасибо. - person Gimly_Gloin; 18.03.2012

arrow_upward
0
arrow_downward

Вуаля, 2 способа сделать это, один не имеет ничего общего с perl:

awk 'BEGIN {while ( i = getline < "reads_list") ar[$i] = $1;} {if ($4 in ar) print $0;}' blast_table > new_blast_table

Мамбо6.pl

#!/usr/bin/perl -w
# purpose:  extract blastX data from a list of readnames. HINT: Make sure your list file only has unique names , that way you save time. 
use strict;
open (DATA,$ARGV[0]) or die ("Usage: ./mambo5.pl BlastXTable readslist");
open (LIST,$ARGV[1]) or die ("Usage: ./mambo5.pl BlastXTable readslist");
my %hash;
my $val;
my $key;
while (<DATA>)
{
    #chomp;
    if(/((.*?)\t(.*?)\t(.*?)\t(.*?)\t(.*?)\t(.*?)\t(.*?)\t(.*?)\t(.*?)\t(.*?)\t(.*?)\t(.*?)\t(.*?)\t(.*?)\t(.*?)\t(.*?)\t(.*?)\t(.*?)\t(.*?)\t(.*?)\t(.*?))$/)
    {
        #print "$1\n";
        $key= $5;#read
        $val= $1;#whole row; notice the brackets around the whole match.
        $hash{$key} .= exists $hash{$key} ? "$val\n" : $val;
    }
    else {
        print "something wrong with format";
    }
}
close (DATA);
open(OUT, "> $ARGV[1]\_out\.txt");

my $readName;

while ( <LIST> )
{
    #########;
    if(/^(.*?)$/)#
    {
        $readName=$1;#
        chomp $readName;
        if (exists $hash{$readName})
        {
            print "$readName\n";
            my $output =$hash{$readName};
            print OUT "$output";
        }
        else 
        {
            #print "it aint workin\n";
        }           
    }
}
close (LIST);
close (OUT);

Oneliner быстрее и, вероятно, лучше, чем мой сценарий, я уверен, что некоторые люди могут найти более простые способы сделать это... Я просто подумал, что выложу это, так как он делает то, что я хочу.

person Gimly_Gloin    schedule 10.03.2012