Changeset 271

Timestamp:

09/30/09 14:05:39 (2 months ago)

Author:

cholt

Message:

make runlog softmask instensitive for old data. Fix fgenesh xdef errors

Files:

• lib/GI.pm (modified) (6 diffs)
• lib/Process/MpiChunk.pm (modified) (1 diff)
• lib/Widget/fgenesh.pm (modified) (6 diffs)
• lib/clean.pm (modified) (1 diff)
• lib/ds_utility.pm (modified) (3 diffs)
• lib/maker/auto_annotator.pm (modified) (26 diffs)
• lib/runlog.pm (modified) (7 diffs)

lib/GI.pm

@@ -352 +352 @@
-   }
-
-
-
+
+
+
+
+
+
+
-      return @holdovers, @keepers;
-   }
@@ -368 +373 @@
-         ($b, $e) = ($e, $b) if $b > $e;
-
-
-
+ && $p_cutoff > $chunk->offset +1
+ && $m_cutoff > $chunk->offset +1
-            ) {
+            push(@{$group_keepers}, $hit);
+         }
+         else {
+            push(@{$group_holdovers}, $hit);
+         }
+      }
+
+      push(@keepers, $group_keepers);
+      push(@holdovers, $group_holdovers);
+   }
+
+   #hit holdovers and keepers are returned in same order given by user
+   return @holdovers, @keepers;
+}
+#-----------------------------------------------------------------------------
+sub process_the_chunk_divide_temp{
+   my $chunk = shift @_;
+   my $split_hit = shift @_;
+   my $hit_groups = \@_; #processed and returned in order given by user
+
+   my $p_hits;
+
+   foreach my $group (@{$hit_groups}) {
+      push(@{$p_hits}, @{$group});
+   }
+
+   my $coors  = PhatHit_utils::to_begin_and_end_coors($p_hits, 'query');
+
+   foreach my $coor (@{$coors}) {
+      $coor->[0] -= $chunk->offset();
+      $coor->[1] -= $chunk->offset();
+      #fix coordinates for hits outside of chunk end   
+      $coor->[0] = $chunk->length if($coor->[0] > $chunk->length);
+      $coor->[1] = $chunk->length if($coor->[1] > $chunk->length);
+      #fix coordinates for hits outside of chunk begin
+      $coor->[0] = 0 if($coor->[0] < 0);
+      $coor->[1] = 0 if($coor->[1] < 0);
+   }
+
+   my $pieces = Shadower::getPieces(\($chunk->seq), $coors, 10);
+   $pieces = [sort {$b->{e} <=> $a->{e}} @{$pieces}];
+
+   my @keepers;
+   my @holdovers;
+
+   my $cutoff = $chunk->length + $chunk->offset - $split_hit;
+   my $p_cutoff = $chunk->length + $chunk->offset + 1;
+
+   foreach my $piece (@{$pieces}) {
+      if ($piece->{e} + $chunk->offset >= $cutoff) {
+         $p_cutoff = $piece->{b} + $chunk->offset;
+      }
+   }
+
+   #too small, all are heldover for next round
+   if ($p_cutoff <= 1 + $chunk->offset) {
+     foreach my $g (@{$hit_groups}){
+         push (@holdovers, $g);
+         push (@keepers, []);
+     }
+      return @holdovers, @keepers;
+   }
+
+   foreach my $group (@{$hit_groups}) {
+      my $group_keepers = [];
+      my $group_holdovers = [];
+
+      foreach my $hit (@{$group}) {
+         my $b = $hit->nB('query');
+         my $e = $hit->nE('query');
+         my $strand = $hit->strand;
+
+         ($b, $e) = ($e, $b) if $b > $e;
+
+         if (($e < $p_cutoff && $p_cutoff > $chunk->offset +1)) {
-            push(@{$group_keepers}, $hit);
-         }
@@ -598 +678 @@
-   #my $FA;
-   #open($FA, "> $t_full"); #full file
+
-   my $wflag = 1; #flag set so warnings gets printed only once 
-   while (my $fasta = $fasta_iterator->nextEntry()) {
@@ -842 +923 @@
-   $command .= " -p $pro";
-   $command .= " -o $out_file";
-
+#
-   $command .= " $in_file";
-
@@ -935 +1016 @@
-
-   my $command  = "$wrap $exe";
-
+#
-   $command .= " $org";
-   $command .= " $in_file";
@@ -2736 +2817 @@
-   ($CTL_OPT{out_name}) = $genome =~ /([^\/]+)$/;
-   $CTL_OPT{out_name} =~ s/\.[^\.]+$//;
+   $CTL_OPT{CWD} = Cwd::cwd();
-   if(! $CTL_OPT{out_base}){
-Cwd::cwd()
+$CTL_OPT{CWD}
-   }
-   mkdir($CTL_OPT{out_base}) if(! -d $CTL_OPT{out_base});

lib/Process/MpiChunk.pm

@@ -1221 +1221 @@
-                $pred_gff_keepers,
-                $model_gff_keepers
-
+_temp
-                                                $CTL_OPT{'split_hit'},
-                                                $blastn_keepers,

lib/Widget/fgenesh.pm

@@ -4 +4 @@
-package Widget::fgenesh;
-use strict;
-
-
-
- $TMP
+#
+#
+
+
-use PostData;
-use FileHandle;
@@ -144 +144 @@
-           }
-           waitpid $pid, 0;
-#
+
-        }
-        else {
@@ -293 +293 @@
-        $command = $wrap . " $command"; #prepend wrapper
-        $command .= " $file_name";
-
+#
-        $command .= ' -exon_table:'.$xdef_file if -e $xdef_file;
-                           
@@ -530 +530 @@
-        my @xdef;
-
-        my $num = (scalar @{$p_pieces} ) + (scalar @{$n_pieces});
-        push @xdef, "$num $s";
-        for (my $i = 0; $i < @{$p_pieces}; $i++){
-                my $p = $p_pieces->[$i];
@@ -537 +535 @@
-                my $c_e = $p->{e} - $offset;
-
-$c_b.' '.$c_e.' '.'1000'
+"$c_b $c_e 100"
-                
-                push(@xdef, $l);
-        }
-
-
-
-
+
-        for (my $i = @{$n_pieces} - 1; $i > 0;  $i--){
-                my $p_r = $n_pieces->[$i]; 
@@ -554 +549 @@
-                next if abs($i_b - $i_e) < 25;
-
-$i_b.' '.$i_e.' -1000'
+"$i_b $i_e -1000"
-
-                push(@xdef, $l);
-        }
+
+        my $num = @xdef;
+        unshift(@xdef, "$num $s");
-
-        # The next part also put ESTs as arbitrary exons.

lib/clean.pm

@@ -216 +216 @@
-
-        #if average basepair hits twice on both query and subject, skip
-2
+1.5
-    }
-    

lib/ds_utility.pm

@@ -14 +14 @@
-@ISA = qw(
-       );
+
+my $CWD = Cwd::getcwd();
+
-#------------------------------------------------------------------------
-#--------------------------- METHODS   ----------------------------------
@@ -32 +35 @@
-   my $self = shift @_;
-   my %CTL_OPTIONS = %{shift @_};
+
+   $CWD = $CTL_OPTIONS{CWD};
-
-   my $out_base = $CTL_OPTIONS{out_base};
@@ -127 +132 @@
-
-   #remove deep directory data so log is relative
-
+($CWD) ? $CWD : 
-   my $entry = join("\t", @F);
-
-   #while loop is used to solve weird incorrect cwd
-   #that happens randomly on the cluster
-   while($cwd ne Cwd::getcwd()){
-       $cwd = Cwd::getcwd();
-   }
-
-   if($entry =~ /\tFINISHED|\tSTARTED|\tDIED|\tSKIPPED|\tRETRY/){

lib/maker/auto_annotator.pm

@@ -605 +605 @@
-    my $the_void         = shift;
-    my $build            = shift;
-ION
+
-    $LOG                 = shift;
-    
@@ -615 +615 @@
-    
-    #reset gene names
-ION
+
-    $SEEN = {};#$GFF_DB->get_existing_gene_names($seq_id);
-    
@@ -625 +625 @@
-                                                  $models,
-                                                  $v_seq_ref,
-ION
-ION
-ION
+
+
+
-                                                  );
-    
@@ -642 +642 @@
-                                 'model_gff',
-                                 $predictions,
-ION
+
-                                 );
-        
@@ -652 +652 @@
-                                                      'model_gff',
-                                                      $the_void,
-ION
+
-                                                      );
-    }
-    
-    #---hint based gene prediction here (includes est2genome)
-ION
+
-        next if($prdr eq 'model_gff' || $prdr eq 'pred_gff');
-        print STDERR "Producing $prdr hint based annotations\n" unless($main::quiet);
@@ -668 +668 @@
-                                 $prdr,
-                                 $predictions,
-ION
+
-                                 );
-        
@@ -678 +678 @@
-                                      $prdr,
-                                      $the_void,
-ION
+
-                                      );
-        
@@ -694 +694 @@
-                            'abinit', #all abinits not just pref_gff
-                            $predictions,
-ION
+
-                            );
-    
@@ -704 +704 @@
-                                   'abinit', #all abinits not_just pred_gff
-                                   $the_void,
-ION
+
-                                   );
-
@@ -873 +873 @@
-    my $annotations = shift;
-    my $out_base = shift;
-ION
+
-
-    print STDERR "Choosing best annotations\n" unless($main::quiet);
@@ -881 +881 @@
-
-    #keep all gff3 passthrough if there's nothing else
-IONS->{_predictor}} == 1 && $CTL_OPTION
+S->{_predictor}} == 1 && $CTL_OPT
-        foreach my $g (@{$annotations->{'model_gff'}}){
-            if($g->{g_strand} == 1){
@@ -891 +891 @@
-        }
-    }
-ION
+
-        #set up lists for plus and minus strands as well as possible mergers
-        #predictor types are processed in the order given by control files
@@ -898 +898 @@
-        my @p_est2g;
-        my @m_est2g;
-ION
+
-            foreach my $g (@{$annotations->{$p}}){
-                if($p ne 'est2genome' && $p ne 'protein2genome' && $g->{g_strand} == 1){
-
+
+
+
-                }
-                elsif($p ne 'est2genome' && $p ne 'protein2genome' && $g->{g_strand} == -1) {
-
+
+
+
-                }
-                elsif($g->{g_strand} == 1){
-
+
+
-                }
-                elsif($g->{g_strand} == -1){
-
+
+
-                }
-                else{
@@ -1155 +1161 @@
-   my $g = shift;
-   
-($g->{AED} + $g->{abAED})/2
+$g->{AED} + ($g->{abAED} * 1/3)
-}
-#------------------------------------------------------------------------
@@ -1207 +1213 @@
-        my $model    = $set->{model};
-        my $gomiph   = $set->{gomiph};
+        my $blastx   = get_selected_types($gomiph,'blastx', 'protein_gff');
+        my $pol_p    = get_selected_types($gomiph,'protein2genome');
+        my $alt_ests = $set->{alt_ests};
+        my $preds    = $set->{preds};
-
-        #------gff passthrough
@@ -1224 +1234 @@
-        if ($predictor eq 'abinit') {
-            next if(! defined $model);
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-            #add UTR to ab-inits
-            my $transcript = pneu($ests, $model, $seq);
@@ -1296 +1330 @@
-        
-        #------default hint based behavior
-        my $blastx   = get_selected_types($gomiph,'blastx', 'protein_gff');
-        my $pol_p    = get_selected_types($gomiph,'protein2genome');
-        my $alt_ests = $set->{alt_ests};
-        my $preds    = $set->{preds};
-        
-        my ($pred_shots, $strand) = get_pred_shot($seq, 
-                                                  $def, 
@@ -1313 +1342 @@
-        my $on_right_strand = get_best_pred_shots($strand, $pred_shots);
-        
-blastx
-> 1 && ! defined $mia && ! @$pol_p && ! @$preds
+single exon blastx suport
+&& ! defined $mia && (!@$pol_p  || (@$pol_p == 1 && $pol_p->[0]->hsps == 1))
-            my $clean  = clean::purge_single_exon_hits($alt_ests);
+            next if(@$clean);
+
-            my $coors  = PhatHit_utils::get_hsp_coors($blastx, 'query');
-            my $pieces = Shadower::getPieces($seq, $coors, 0);
-            
-! @$clean && 
+
-                my $keep;
-                foreach my $h (@$on_right_strand){
-                    my $set = get_overlapping_hits($h, $predictions);
-                    $keep = 1 if(@$set);
+                    
+                    my $abAED = shadow_AED::get_abAED($set, $h);
+                    $keep = 0 if($abAED > 0.3);
+
+                    my $AED = shadow_AED::get_AED($blastx,$h);
+                    $keep = 0 if($AED > 0.3);
-                }
-
@@ -1355 +1392 @@
-   my $set         = shift;
-   my $predictor   = shift;
-ION
+
-   my $LOG         = shift;
-
-   if($predictor eq 'snap'){
-IONS->{snap}.' '.$CTL_OPTION
+S->{snap}.' '.$CTL_OPT
-      return Widget::snap::get_pred_shot($seq, 
-                                         $def, 
@@ -1365 +1402 @@
-                                         $the_void, 
-                                         $set, 
-ION
+
-                                         $pred_command,
-ION
+
-                                         $LOG
-                                        );
-   }
-   elsif($predictor eq 'augustus'){
-IONS->{augustus}.' --species='.$CTL_OPTION
+S->{augustus}.' --species='.$CTL_OPT
-      return Widget::augustus::get_pred_shot($seq, 
-                                             $def, 
@@ -1378 +1415 @@
-                                             $the_void, 
-                                             $set, 
-ION
+
-                                             $pred_command,
-ION
+
-                                             $LOG
-                                            );
-   }
-   elsif($predictor eq 'fgenesh'){
-IONS->{fgenesh}.' '.$CTL_OPTION
+S->{fgenesh}.' '.$CTL_OPT
-      return Widget::fgenesh::get_pred_shot($seq, 
-                                            $def, 
@@ -1391 +1428 @@
-                                            $the_void, 
-                                            $set, 
-ION
+
-                                            $pred_command,
-ION
+
-                                            $LOG
-                                           );
-   }
-   elsif($predictor eq 'twinscan'){
-ION
+
-      return Widget::twinscan::get_pred_shot($seq, 
-                                             $def, 
@@ -1404 +1441 @@
-                                             $the_void, 
-                                             $set, 
-ION
+
-                                             $pred_command,
-ION
+
-                                             $LOG
-                                            );
@@ -1427 +1464 @@
-        my $p_base       = shift;
-        my $the_void     = shift;
-ION
+
-
-        my $transcript_seq  = get_transcript_seq($f, $seq);
@@ -1565 +1602 @@
-   my $predictor    = shift;
-   my $the_void     = shift;
-ION
+
-
-   #fix weird sequence names
@@ -1688 +1725 @@
-         my $p_base = defined($f->{set_id}) ? $p_bases{$f->{set_id}} : undef;
-
-ION
-
-ION
-ION
+
+
+
+
-                                             );
-

lib/runlog.pm

@@ -77 +77 @@
-
-my %SEEN;
+my $CWD = Cwd::getcwd();
-
-#-------------------------------------------------------------------------------
@@ -107 +108 @@
-   print STDERR "\n\n\n--Next Contig--\n\n" unless($main::quiet);
-
+   $CWD = $self->{CTL_OPTIONS}->{CWD};
-   my $min_contig = $self->{CTL_OPTIONS}->{min_contig};
-   my $length = $self->{params}->{seq_length};
@@ -197 +199 @@
-        if($continue_flag >= 0 || $continue_flag == -1){
-            #CHECK CONTROL FILE OPTIONS FOR CHANGES
-
-
-
-
-
-
-
+
-
-            foreach my $key (@ctl_to_log) {
@@ -218 +214 @@
-
-                    #temp lamprey
-
-
-
-
-
-
-
+
+
-
-                    #only continue if change not already happening
@@ -269 +260 @@
-                if($key eq 'organism_type' && ! $log_val){
-                    $log_val = 'eukaryotic';
+                }
+
+                #softmask was always true before and not logged
+                if($key eq 'softmask' && $log_val eq ''){
+                    $log_val = 1;
-                }
-                
@@ -614 +610 @@
-
-   #log control file options
-
-
-
-
-
-
-
+
- 
-   foreach my $key (@ctl_to_log) {
@@ -650 +640 @@
-
-   my $log_file = $self->{file_name};
-
-
-
-
-
-
-
+
-
-   #this line hides unnecessarilly deep directory details