附件：Easyseq源码（[3]maplus）

#!/usr/bin/perl
 use threads; 
 use Thread::Queue;
 use Bio::SeqIO;
 use Bio::Seq;
 use Cwd;
 use Tie::File;
     $q = Thread::Queue->new();
     $q -> limit = 0;
     $location_dir=getcwd;
     $location_dir=~s/\//\\/;
     @file = glob("$location_dir\\*.filter");
    foreach $file(@file){
    open(DATA,"<$file") or die "Can't open $file";
     @list = <DATA>;
    foreach $list (@list){
     @list_elm=split(" ",$list);
     $pairs=1;
     if($list_elm[0]=~/primer_pairs_count/){
     $primer_pairs_count=$list_elm[1];
    }
    while($pairs<$primer_pairs_count+1){
      if($list_elm[0]=~/F_$pairs/){
          $ORFF[$pairs]=$list_elm[1];
          print "F_$pairs : $ORFF[$pairs]\n";
      }
      if($list_elm[0]=~/P_$pairs/){
          $ORFP[$pairs]=$list_elm[1];
          print "P_$pairs : $ORFP[$pairs]\n";
      }
      if($list_elm[0]=~/R_$pairs/){
          $rev_seq= reverse ($list_elm[1]);
          $rev_seq=~tr/[ATCGatcg]/[TAGCtagc]/;
          $ORFR[$pairs]=$rev_seq;
          print "R_$pairs : $ORFR[$pairs]\n";
      }
      $pairs=$pairs+1;
    }
     if($list_elm[0]=~/mapping_max_nohit/){
     $mapping_length=$list_elm[1];
    } 
    if($list_elm[0]=~/max_threads/){
     $max_threads=$list_elm[1];
    }
    if($list_elm[0]=~/max_target_length/){
     $max_target_length=$list_elm[1];
    }
    }
    close DATA;
    } 
 sub produce {       
    my $location_dir=getcwd;
    $location_dir=~s/\//\\/;
    my @file = glob("$location_dir\\classified\\*.fa");
      my $dirpath="$location_dir\\mapped";
    mkdir ($dirpath); 
      foreach my $file(@file){
      print " successful import $file\n";
    my $fa=Bio::SeqIO->new(-file=>$file,-format=>'fasta');
    while($seq_obj=$fa->next_seq){
        my $id=$seq_obj->id;
        my $seq=$seq_obj->seq;
        my $desc=$seq_obj->desc;
        my @r=($id, $seq, $desc, $file);
        $r=\@r;
        $q->enqueue($r);        
    }
    }
}
 sub consume {     
    my($id, $seq, $desc, $file) = @_;
    print("Sequence number $uid $file $id\n");        
    $consume = threads -> create(\&score, $id, $seq, $desc, $file);
    $uid = $consume->tid();      
}
 sub score {     
    my($id, $seq, $desc, $file)=@_;
    $pai=1; 
    $mm=0;
    $n=0;
  
    LOOP:do{
      if($n==0){
        $name="Target_Forward_Primer[$pai]";
        $IN=$ORFF[$pai];
      }
      if($n==1){
        $name="Target_Reverse_Primer[$pai]";
        $IN=$ORFR[$pai];
      }
      if($n==2){
        $name="Target_Probe_Primer[$pai]";
        $IN=$ORFP[$pai];
      }
      @IN=split("",$IN);
      $c=0;
      $c1=0;
      $c2=0;
      $c3=0;
      $g=0;
      %score=();
      $score=0;
      %seq=();
      #print "loading $id successiful,matching $name\n";
      @seq=split('', $seq);
      $seq_length=@seq;
      $primer_length=@IN;
      while($c3<=$seq_length-$primer_length){
        $end=$c3+$primer_length-1;
      while($c2<=$primer_length-1){
        if($seq[$c1] eq $IN[$c2]){
          $score=$score+2;
          $c2=$c2+1;
          $c1=$c1+1;
        }else{
          $c2=$c2+1;
          $c1=$c1+1;
        }
      }
          $score{$c3}=$score;
          $c2=0;
          $c3=$c3+1;
          $c1=$c3;
          $score=0;
      }
          @goat= values %score;
      foreach (@goat){
        if($_>$g) {
          $g=$_;
        }
      }
      until( $g eq $score{$c} ){
        $c = $c + 1;
        }
      $L=2*($primer_length-$mapping_length);
      $end=$c+$primer_length-1;
      @now=@seq[$c..$end];
      $now=join('', @now);  
    if($score{$c}>=$L){ 
      $mmm=1;
      $mm=$mm+$mmm;
    }else{
      $mmm=0;
      $mm=$mm+$mmm;
    }  
    $seq2[$n]=$now;
    if($n<2){       
          #print "Sequences $id,matching $name,best match region is $now \n";
          #print "$mm\n";
          
          $n=$n+1;
          goto LOOP;        
    }else{
      if($mm==3){
          $n=0;
          $mm=0;
          $seq =~ /$seq2[0]/;
          $m2=$&;
          $m3=$';
          $m3 =~ /$seq2[1]/;
          $m5=$&;
          $m6=$`;
          $m6 =~ /$seq2[2]/;
          $m7=$&;
          $file_txt=$file;    
              $file_txt=~s/\.fa//;
              $file_txt=~s/classified/mapped/;
          #print "$file_txt\n";
          $Target_amplicon=$m2.$m6.$m5;
          #print  "$m2$m7$m5\n";
          #print  "$m2$m6$m5\n";
        if(length($Target_amplicon)<=$max_target_length-1){                    
          open WH, ">> $file_txt\_Primer$pai.hit";
          print WH ">$id $desc\n";
          print WH "$m2$m7$m5\n";
          #print  "$m2$m7$m5\n";
          close WH;                
          open WH, ">> $file_txt\_Primer$pai.fa";
          print WH ">$id $desc\n";
          print WH "$m2$m6$m5\n";
          close WH;
          #print  "$m2$m6$m5\n";       
        }
      }else{
          $n=0;
          $mm=0;
      }
          if($pai<$primer_pairs_count){
            $pai=$pai+1;         
            goto LOOP;
          }else{
            $pai=1;
          }
    
    }#提取最佳匹配序      
    }#LOOP      
  } 

 $produce = threads->create({'exit'=>'thread_only'},\&produce, "Main");
 $produce -> join();
 $queue_length=$q->pending();
 while($queue_length>1){  
  $queue_length=$q->pending();
 if($queue_length>=$max_threads){ 
    if(scalar(threads->list())<$max_threads) {
        my $r=$q->dequeue();
        my($id, $seq, $desc, $file)=@$r; 
        consume($id, $seq, $desc, $file); 
      }else{
        foreach $thread(threads->list(threads::all)) {
        $thread->join();
        }
      }
 }else{
    if($queue_length>1) {
        my $r=$q->dequeue();
        my($id, $seq, $desc, $file)=@$r; 
        consume($id, $seq, $desc, $file);
      }else{
        foreach $thread(threads->list(threads::all)) {
        $thread->join();
        }
      }       
 }     
 }
        $r=$q->extract(); 
        $q->end();
        my($id, $seq, $desc, $file)=@$r; 
        score($id, $seq, $desc, $file);
        $uid = $consume->tid();
        print("consume $uid $file $id\n");