Использование remote BLAST в автоматическом режиме

Задача: для набора сиквенсов запустить remote BLAST из NCBI и получить окрестность для хита. Реализация на Perl + BioPerl. Вот попытка реализации:

1. чтение сиквенсов из текстового файла (каждая строка - сиквенс). Можно легко адаптировать на FASTA-формат.

2. Запуск remote blastn с сервера NCBI с заданием:
    - минимального e-value
    - база данных nr, ограничение txid10239 [ORGN] - все вирусы

3. Результат BLAST сохраняется в виде XML. Для каждого хита скачивается последовательность в которой содержится хит, вырезается окрестность заданной длинны ($DELTA), содержащая хит.

4. Сохранение в текстовый файл - один для каждого исходного сиквенса. В каждом файле может быть несколько хитов.

Код переделан из demo для BioPerl и мало оптимизирован. Например, последовательность в которой есть хит(т.е. геном) скачивается каждый раз заново, а хорошо бы это кэшировать.

Все параметры указаны в тексте. Используется дополнительный модуль misc.pm - только ради функции trim() - обрезать пробелы в строке. Это можно убрать.

Если этот код кому-то окажется полезным - постараюсь ответить на вопросы (в комментах)

       

#
# BLAST nucleotide sequences against Genbank
#
use List::Util qw[min max];
use Bio::SearchIO;
use Bio::Tools::Run::RemoteBlast;
use Bio::DB::GenBank;
use misc;

my $infile="seqs.txt"; # file with sequences
my $prog = 'blastn';
my $blastdb = 'nr';
my $e_val= '1';
$DELTA=100; #how many nucleotides before and after hit to preserve

$db = new Bio::DB::GenBank; #connect to NCBI
my @params = ( '-prog' => $prog,
'-data' => $blastdb,
'-expect' => $e_val,
'-readmethod' => 'xml' );
my $factory = Bio::Tools::Run::RemoteBlast->new(@params);

#SEE http://www.ncbi.nlm.nih.gov/BLAST/Doc/urlapi.html for more params
$factory->retrieve_parameter('FORMAT_TYPE', 'XML'); # tells NCBI to send XML back

#change a query paramter
$Bio::Tools::Run::RemoteBlast::HEADER{'ENTREZ_QUERY'} = 'txid10239 [ORGN]'; # all viruses
#$Bio::Tools::Run::RemoteBlast::HEADER{'ENTREZ_QUERY'} = 'txid629 [ORGN]'; # Yersinia

$Bio::Tools::Run::RemoteBlast::HEADER{'FILTER'} = 'L';
$Bio::Tools::Run::RemoteBlast::HEADER{'WORD_SIZE'} = '7';
$Bio::Tools::Run::RemoteBlast::HEADER{'FORMAT_TYPE'} = 'XML';
#see also OTHER_ADVANCED for specifying word size etc

#howto: change a retrieval parameter
#$Bio::Tools::Run::RemoteBlast::RETRIEVALHEADER{'DESCRIPTIONS'} = 1000;
#remove a parameter
#delete $Bio::Tools::Run::RemoteBlast::HEADER{'FILTER'};

open(INFILE,$infile) or die "error opening file";
while () {
my $str=$_; chomp($str);
my $pos=index $str, "#";
if ($pos==0) { next; } #skipping commentaries

$seq=trim($str);

$outfile=$seq.".txt";
# to allow re-run
if (-e $outfile) { print "file $outfile already exists - skipping\n"; next; }

open(OUTFILE,">$outfile");
$old_fh=select(OUTFILE);$|=1;select($old_fh); #disable buffering

my $seqobj = Bio::Seq->new(-seq =>$seq);
my $r = $factory->submit_blast($seqobj);

while ( my @rids = $factory->each_rid ) { #RID - request ID
# my @rids = $factory->each_rid ; #RID - request ID
$rid= $rids[-1]; #taking last one

my $rc = $factory->retrieve_blast($rid);
if( !ref($rc) ) { if($rc<0) { $factory->remove_rid($rid); } print "."; sleep 1;}
else {
my $result = $rc->next_result();
my $filename = $seq.".blres.xml";
$factory->save_output($filename); # save balst result to xml

# parsing xml file
my $searchio = new Bio::SearchIO( -format => 'blastxml', -file => $filename );
while ( my $res = $searchio->next_result() ) {
while( my $hit = $res->next_hit ) { # Bio::Search::Hit::HitI object
my $query_id=$hit->accession;
my $hitdescr=$hit->description;
while( my $hsp = $hit->next_hsp ) { # Bio::Search::HSP::HSPI object
$start=$hsp->{HIT_START};
$end=$hsp->{HIT_END};
$evalue=$hsp->{EVALUE};

#NB: do this before retrieving the sequence!
if ($start>$end) { print "Error: start $start > end $end Swapping them\n";
my $tmp=$start;
$start=$end; $end=$tmp;
}

$delta=$DELTA; #NB reassign it here each time

print "retrieving $query_id...";
$seq = $db->get_Seq_by_acc($query_id);
print "done\n";

if ((($start-$delta)<=0)||(($end+$delta)>length($seq->seq)))
{ $delta = min($start-1, length($seq->seq)-$end, 100); } #change delta

# for debug:
# print "debug start= $start end = $end delta=$delta seqlength=".length($seq->seq)." ";
# printing output
print OUTFILE "$query_id\t$hitdescr\n";
print OUTFILE "$start\t$end\t$hsp->{IDENTICAL}\t$evalue\n";

$subseq = $seq->primary_seq->subseq($start-$delta,$end+$delta);
print OUTFILE "$subseq\n";

for(my $i=0;$i<$delta;$i++) { print OUTFILE " ";};
print OUTFILE "$hsp->{QUERY_SEQ}\n";
for(my $i=0;$i<$delta;$i++) { print OUTFILE " ";};
print OUTFILE "$hsp->{HOMOLOGY_SEQ}\n";
for(my $i=0;$i<$delta;$i++) { print OUTFILE " ";};
print OUTFILE "$hsp->{HIT_SEQ}\n";
}
print OUTFILE "\n"; #separator between records
} # for each hsp
} #$result = $searchio->next_result()
last;
}
}# while ( my @rids = $factory->each_rid )
print "\nfinished for sequence $seq\n";
close(OUTFILE);
}#reading $infile
close(INFILE);
print "#work is finished";

Комментарии   

 
0 #1 emilio 10.10.2014 03:59
hello, I am very interested to your script, but I need to load the module named misc.pm Where can I find it? Best regards,
Emilio.

Привет! Я хочу использовать ваш скрипт, но я не знаю как скачать модуль misc.pm. Где я могу его найти? С уважением, Эмилио.
Цитировать
 

Добавить комментарий


Защитный код
Обновить