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

Обновлено 29.05.2009 Автор: Administrator

Задача: для набора сиквенсов запустить 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";

< Предыдущая		Следующая >