Changeset 160

Timestamp:

03/16/09 14:18:16 (8 months ago)

Author:

cholt

Message:

Maker Updates for OOmycete genome, clone hsps, hits in join.pm

Files:

• MPI/mpi_maker (modified) (1 diff)
• lib/GI.pm (modified) (12 diffs)
• lib/Process/MpiChunk.pm (modified) (1 diff)
• lib/maker/auto_annotator.pm (modified) (20 diffs)
• lib/maker/join.pm (modified) (8 diffs)
• lib/runlog.pm (modified) (4 diffs)
• lib/shadow_AED.pm (modified) (1 diff)

MPI/mpi_maker

@@ -250 +250 @@
-               "predictor=s" =>\$OPT{predictor},
-               "retry=i" =>\$OPT{retry},
-
+
+
-               "quiet" =>\$main::quiet,
-               "CTL" => sub {GI::generate_control_files() if($rank == $root); MPI_Finalize(); exit(0);},

lib/GI.pm

@@ -694 +694 @@
-   my $snaphmm = $CTL_OPT->{snaphmm};
-   
-not $exe
-
+! grep {/snap/} @{$CTL_OPT->{_run}}
+
-   my %params;
-   my $out_file = "$the_void/$seq_id\.all\.snap";
-
+
-   $LOG->add_entry("STARTED", $out_file, "");   
-
@@ -740 +740 @@
-   my $org = $CTL_OPT->{augustus_species};
-
-not $exe
+! grep {/augustus/} @{$CTL_OPT->{_run}}
-
-   my %params;
@@ -786 +786 @@
-   my $org = $CTL_OPT->{fgenesh_par_file};
-
-not $exe
+! grep {/fgenesh/} @{$CTL_OPT->{_run}}
-
-   my %params;
@@ -832 +832 @@
-   my $org = $CTL_OPT->{twinscan_species};
-
-not $exe
+! grep {/twinscan/} @{$CTL_OPT->{_run}}
-
-   my %params;
@@ -2017 +2017 @@
-      $CTL_OPT{'rm_gff'} = '';
-      $CTL_OPT{'predictor'} = 'est2genome';
-
+model_
-      $CTL_OPT{'snaphmm'} = 'fly';
-      $CTL_OPT{'augustus_species'} = 'fly';
@@ -2038 +2038 @@
-      $CTL_OPT{'clean_try'} = 0;
-      $CTL_OPT{'TMP'} = '';
-eval_pred
-eval_pred'} = 'snap' if($main::eva)
+run
+unmask'} = 1
-      $CTL_OPT{'clean_up'} = 0;
-   }
@@ -2184 +2184 @@
-
-   #required for evaluator to work
-
+model_
-   $CTL_OPT{model_pass} = 1 if($main::eva);
-
@@ -2190 +2190 @@
-   $CTL_OPT{predictor} =~ s/\s+//g;
-   my @predictors = split(',', $CTL_OPT{predictor});
-
-
+
+
-   foreach my $p (@predictors) {
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-      }
-   }
@@ -2283 +2291 @@
-   push (@infiles, 'RepeatMasker') if($CTL_OPT{model_org});
-   push (@infiles, 'rmlib') if ($CTL_OPT{rmlib});
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
-   push (@infiles, 'fathom') if ($CTL_OPT{enable_fathom});
-   push (@infiles, 'rm_gff') if($CTL_OPT{rm_gff});
@@ -2306 +2309 @@
-
-   #uniq the array
-my %uniq
+%uniq = ()
-   foreach my $in (@infiles) {
-      $uniq{$in}++;
@@ -2493 +2496 @@
-   print OUT "#-----Gene Prediction Options\n" if(!$ev);
-   print OUT "#-----Evaluator Ab-Initio Comparison Options\n" if($ev);
-
-
+
+
+
-   print OUT "snaphmm:$O{snaphmm} #SNAP HMM model\n";
-   print OUT "augustus_species:$O{augustus_species} #Augustus gene prediction model\n";
@@ -2509 +2513 @@
-   print OUT "\n";
-   print OUT "#-----Maker Specific Options\n";
-
+, 1 = yes, 0 = no
-   print OUT "max_dna_len:$O{max_dna_len} #length for dividing up contigs into chunks (larger values increase memory usage)\n";
-   print OUT "min_contig:$O{min_contig} #all contigs from the input genome file below this size will be skipped\n";

lib/Process/MpiChunk.pm

@@ -648 +648 @@
-
-            #==ab initio predictions here
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-            #==QRNA noncoding RNA prediction here

lib/maker/auto_annotator.pm

@@ -24 +24 @@
-#use evaluator::AED;
-use shadow_AED;
-use Storable qw(dclone);
-
-$Storable::forgive_me = 1; 
@@ -582 +581 @@
-                                    $v_seq_ref,
-                                    $def,
-
-
+
-                                    $the_void,
-                                    $CTL_OPTIONS
-                                   );
-           
-
-  
-  
-  
-  
-   '
-  
-  
- 
+model_
+       
+       
+       
+       
+        'model_
+       
+       
+       
-        }
-
-        #---hint based gene prediction here
-        foreach my $prdr (@{$CTL_OPTIONS->{_predictor}}){
-
+model_
-           print STDERR "Producing $prdr hint based annotations\n" unless($main::quiet);
-
@@ -609 +607 @@
-                                    $v_seq_ref,
-                                    $def,
-                                    'bx',
-                                    $prdr,
-                                    $CTL_OPTIONS
@@ -635 +632 @@
-                                $v_seq_ref,
-                                $def,
-                                'pr',
-                                'abinit',
-                                $CTL_OPTIONS
@@ -687 +683 @@
-}
-#------------------------------------------------------------------------
+sub gene2all_phathits {
+    my $g = shift;
+    
+    my @hits;
+    foreach my $t (@{$g->{t_structs}}){
+        push (@hits, $t->{hit});
+    }
+
+    return \@hits;
+}
+#------------------------------------------------------------------------
-#this subrutine returns finished MAKER annotations
-sub best_annotations {
@@ -702 +709 @@
-        #set up lists for plus and minus strands as well as possible mergers
-        #predictor types are processed in the order given by control files
+        my @p_bag;
+        my @m_bag;
-        foreach my $p (@{$CTL_OPTIONS->{_predictor}}){
-            foreach my $g (@{$annotations->{$p}}){
-                next unless($g->{AED} < 1 || $p eq 'gff');
-
-                if($g->{g_strand} == 1){
-
+
+
+
-                }
-
-
+
+
+
+
-                }
-                else{
@@ -717 +728 @@
-            }
-        }
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-        #remove low scoring overlaping genes
-        #@$p_list  = sort {crit1($a) <=> crit1($b) || crit2($a) <=> crit2($b) || crit3($a) <=> crit3($b)} @$p_list;
@@ -733 +766 @@
-        my $key = $CTL_OPTIONS->{_predictor}->[0];
-        foreach my $g (@{$annotations->{$key}}){
-
+model_
-        }
-    }
@@ -803 +836 @@
-    my $v_seq        = shift;
-    my $def          = shift;
-    my $id           = shift;
-    my $predictor    = shift;
-    my $CTL_OPTIONS  = shift;
@@ -816 +848 @@
-
-        #------gff passthrough and ab-init passthrough
-
+model_
-            next if(! defined $model);
-            my $transcript = $model;
@@ -822 +854 @@
-            #add UTR to ab-inits
-            if($predictor eq 'abinit' && defined($ests)){
-
-
+
-                if (! clean::is_identical_form($utr_trans, $transcript)){
-                    push(@transcripts, [$utr_trans, $set, $transcript]);
@@ -834 +865 @@
-            next;
-        }
+
+        #------ab-init passthrough
+        if ($predictor eq 'model_gff' || $predictor eq 'abinit') {
+            next if(! defined $model);
+            my $transcript = $model;
+            
+            #add UTR to ab-inits
+            if($predictor eq 'abinit' && defined($ests)){
+                my $utr_trans = pneu($ests, $transcript, $seq);
+                if (! clean::is_identical_form($utr_trans, $transcript)){
+                    push(@transcripts, [$utr_trans, $set, $transcript]);
+                }
+            }
+            
+            push(@transcripts, [$transcript, $set, undef]);
+            
+            $i++;
+            next;
+        }
-        
-        #------est2genome
-        if ($predictor eq 'est2genome') {
-            next if (! defined $mia);
-
-
+
-            push(@transcripts, [$transcript, $set, $mia]) if $transcript;
-            $i++;
@@ -854 +903 @@
-        my ($pred_shots, $strand) = get_pred_shot($seq, 
-                                                  $def, 
-d.'.'.$i
+
-                                                  $the_void, 
-                                                  $set,
@@ -881 +930 @@
-                my $transcript = pneu($ests, $copy, $seq);
-                if (defined($transcript)) {
-
-    
+;
+
-                else {
-
-    
+;
+
-            }
-            elsif (defined($pred_shot) && !defined($mia)) {
@@ -994 +1043 @@
-                  );
-        
-        
-        #evidence AED
-        my $AED = shadow_AED::get_AED(\@bag, $f);
-
-        #average the AED from ab-inits
-        my $abinit_AED = 0;
-        foreach my $ab (@$abinits){
-            $abinit_AED += shadow_AED::get_AED([$ab], $f);
-        }
-        $abinit_AED = (@$abinits) ? $abinit_AED/@$abinits : 1;
-
-        #combined abinit and evidence AED
-        my $X_AED = ($abinit_AED + $AED) / 2;
-
-        my $qi    = maker::quality_index::get_transcript_qi($f,$evi,$offset,$len_3_utr,$l_trans);
@@ -1023 +1061 @@
-                        't_qi'     => $qi,
-                        'AED'      => $AED,
-                        'ab_AED'   => $abinit_AED,
-                        'X_AED'    => $X_AED,
-                        'evi'      => $evi
-                    };
@@ -1031 +1067 @@
-}
-#------------------------------------------------------------------------                                                           
-overlapping_hits
-eat
+hits_overlapping_gene
+g
-   my $hits = shift;
-
+
+
+
+
-   my @keepers;
-   foreach my $hit (@{$hits}){
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
-   return \@keepers;
-}
@@ -1090 +1137 @@
-   
-   my $careful_clusters = [];
-
+model_
-      my %index;
-      my $i = 0;
@@ -1108 +1155 @@
-   elsif ($predictor eq 'abinit' ) {
-      foreach my $t (@transcripts) {
-push(@{$careful_clusters}, [$t])
+ push(@{$careful_clusters}, [$t]);
-      }
-   }
@@ -1136 +1183 @@
-      
-      my $g_name;
-
+model_
-         $g_name = $c->[0]->{gene_name}; #affects GFFV3.pm
-         if ($g_name =~ /^maker-$seq_id|$seq_id-abinit/) {

lib/maker/join.pm

@@ -270 +270 @@
-        #print "PPPPPPPPPPPPP\n";
-
-
+#
-
-        my $sorted_g = PhatHit_utils::sort_hits($g, 'query');
@@ -289 +289 @@
-        foreach my $hsp (@{$sorted_5}){
-                last if $i ==  $b_5->{five_j};
-$hsp
+clone_hsp($hsp)
-                $i++;
-        }
-        #-- merge the 5-prime overlaping features
-        if ( !defined($join_offset_5) || $join_offset_5 == -1){
-$sorted_g->[0]
+clone_hsp($sorted_g->[0])
-        }
-        else {
@@ -308 +308 @@
-        #-- push on the gene pred's interior hsps
-        for (my $i = 1; $i < @{$sorted_g} -1; $i++){
-$sorted_g->[$i]
+clone_hsp($sorted_g->[$i])
-        }
-
@@ -314 +314 @@
-        my $omega = $#{$sorted_g};
-        if (!defined($join_offset_3) || $join_offset_3 == -1){
-$sorted_g->[$omega]
+clone_hsp($sorted_g->[$omega])
-        }
-        else {
@@ -337 +337 @@
-        if (defined($b_3->{three_j})){
-                for (my $i = $b_3->{three_j} + 1; $i < @{$sorted_3}; $i++){
-$sorted_3->[$i]
+clone_hsp($sorted_3->[$i])
-                } 
-        }
@@ -378 +378 @@
-}
-#------------------------------------------------------------------------
+sub clone_hsp{
+    my $hsp = shift;
+
+    my @args;
+
+    push(@args, '-query_start');
+    push(@args, $hsp->{QUERY_START});
+
+    push(@args, '-query_seq');
+    push(@args, $hsp->{QUERY_SEQ});
+
+    push(@args, '-score');
+    push(@args, $hsp->{SCORE});
+
+    push(@args, '-homology_seq');
+    push(@args, $hsp->{HOMOLOGY_SEQ});
+
+    push(@args, '-hit_start');
+    push(@args, $hsp->{HIT_START});
+
+    push(@args, '-hit_seq');
+    push(@args, $hsp->{HIT_SEQ});
+
+    push(@args, '-hsp_length');
+    push(@args, $hsp->{HSP_LENGTH});
+
+    push(@args, '-identical');
+    push(@args, $hsp->{IDENTICAL});
+
+    push(@args, '-hit_length');
+    push(@args, $hsp->{HIT_LENGTH});
+
+    push(@args, '-query_name');
+    push(@args, $hsp->{QUERY_NAME});
+
+    push(@args, '-algorithm');
+    push(@args, $hsp->{ALGORITHM});
+
+    push(@args, '-bits');
+    push(@args, $hsp->{BITS});
+
+    push(@args, '-evalue');
+    push(@args, $hsp->{EVALUE});
+
+    push(@args, '-pvalue');
+    push(@args, $hsp->{PVALUE});
+
+    push(@args, '-query_length');
+    push(@args, $hsp->{QUERY_LENGTH});
+
+    push(@args, '-query_end');
+    push(@args, $hsp->{QUERY_END});
+
+    push(@args, '-conserved');
+    push(@args, $hsp->{CONSERVED});
+
+    push(@args, '-hit_name');
+    push(@args, $hsp->{HIT_NAME});
+
+    push(@args, '-hit_end');
+    push(@args, $hsp->{HIT_END});
+
+    push(@args, '-query_gaps');
+    push(@args, $hsp->{QUERY_GAPS});
+
+    push(@args, '-hit_gaps');
+    push(@args, $hsp->{HIT_GAPS});
+
+    my $REF = ref($hsp);
+    my $clone = new $REF(@args);
+
+    $clone->{queryName} = $hsp->{queryName};
+    #-------------------------------------------------
+    # setting strand because bioperl is all messed up!
+    #------------------------------------------------
+
+    $clone->{_strand_hack}->{query} = $hsp->{_strand_hack}->{query};
+    $clone->{_strand_hack}->{hit}   = $hsp->{_strand_hack}->{query};
+
+    return $clone;
+}
+#------------------------------------------------------------------------
-sub merge_hsp {
-        my $ext   = shift;
@@ -384 +466 @@
-        my $type  = shift;
-
-
-        my $offset;
-        my $class;
@@ -398 +479 @@
-        my $sorted_ext = PhatHit_utils::sort_hits($ext->{f}, 'query');
-        my $ext_hsp = $sorted_ext->[$offset];
+
+        $g_hsp = clone_hsp($g_hsp);
-
-        if     ($class eq '1'){

lib/runlog.pm

@@ -24 +24 @@
-                  'rm_gff',
-                  'predictor',
-eval_pred
+run
-                  'sort_base',
-                  'snaphmm',
@@ -58 +58 @@
-                  'ep_score_limit',
-                  'en_score_limit',
-
+
+
-                 );
-
@@ -181 +182 @@
-                   $log_val =~ s/.*\/(maker\/data\/te_proteins.fasta)$/$1/;
-                }
+                elsif($key eq 'run'){
+                   #don't care about order
+                   my @set = sort @{$logged_vals{CTL_OPTIONS}{_run}};
+                   $log_val = join(",", @set);
+                }
-            }
-            
@@ -190 +196 @@
-                  #don't care about absolute location
-                  $ctl_val =~ s/.*\/(maker\/data\/te_proteins.fasta)$/$1/;
+               }
+               elsif($key eq 'run'){
+                  #don't care about order
+                  my @set = sort @{$CTL_OPTIONS{_run}};
+                  $ctl_val = join(",", @set);
-               }
-            }

lib/shadow_AED.pm

@@ -4 +4 @@
-package shadow_AED;
-use strict;
+
+sub get_abAED{
+    my $hits = shift;
+    my $tran = shift;
+
+    my $sum;
+    foreach my $h (@$hits){
+        $sum += get_AED([$h], $tran);
+    }
+    return 1 if(! defined $sum);
+    return $sum/@$hits;
+}
-
-sub get_AED {