Changeset 151

Timestamp:

02/20/09 12:38:14 (9 months ago)

Author:

cholt

Message:

maker update CDS has stop codon and other opts

Files:

• Apollo/gff3.tiers (modified) (1 diff)
• MPI/mpi_maker (modified) (4 diffs)
• bin/gff3_merge (modified) (7 diffs)
• bin/maker (modified) (3 diffs)
• lib/Bio/Search/Hit/PhatHit/Base.pm (modified) (1 diff)
• lib/CGL/TranslationMachine.pm (modified) (1 diff)
• lib/Dumper/GFF/GFFV3.pm (modified) (4 diffs)
• lib/Fasta.pm (modified) (1 diff)
• lib/GFFDB.pm (modified) (10 diffs)
• lib/GI.pm (modified) (10 diffs)
• lib/Process/MpiChunk.pm (modified) (6 diffs)
• lib/Widget/augustus.pm (modified) (3 diffs)
• lib/clean.pm (modified) (2 diffs)
• lib/cluster.pm (modified) (1 diff)
• lib/evaluator/funs.pm (modified) (1 diff)
• lib/evaluator/quality_index.pm (modified) (1 diff)
• lib/maker/auto_annotator.pm (modified) (31 diffs)
• lib/maker/join.pm (modified) (1 diff)
• lib/shadow_AED.pm (modified) (2 diffs)

Apollo/gff3.tiers

@@ -515 +515 @@
-
-[Type]
+label : GENMARK
+tiername : Gene Prediction
+datatype : genemark:dummy
+datatype : GeneMark
+datatype : pred_gff:GeneMark
+glyph : DrawableResultFeatureSet
+color : 103,155,104
+column : SCORE
+column : GENOMIC_RANGE
+column : query_frame
+column : GENOMIC_LENGTH
+sortbycolumn : GENOMIC_RANGE
+
+
+[Type]
-label : SNAP
-tiername : Gene Prediction

MPI/mpi_maker

@@ -242 +242 @@
-my $root = 0; #define root node (only changed for debugging)
-
+#---Process options on the command line 
+try{
+    GetOptions("RM_off|R" => \$OPT{R},
+               "force|f" => \$OPT{force},
+               "genome|g=s" => \$OPT{genome},
+               "cpus|c=i" => \$OPT{cpus},
+               "predictor=s" =>\$OPT{predictor},
+               "retry=i" =>\$OPT{retry},
+               "quiet" =>\$main::quiet,
+               "CTL" => sub {GI::generate_control_files() if($rank == $root); MPI_Finalize(); exit(0);},
+               "help|?" => sub {print $usage if($rank == $root); MPI_Finalize(); exit(0)}
+               );
+}
+catch Error::Simple with{
+    my $E = shift;
+    
+    print STDERR $E->{-text};
+    die "\n\nMaker failed parsing command line options!!\n\n";
+};
+
-#--------------------------------------
-#---------PRIMARY MPI PROCCESS---------
@@ -248 +268 @@
-#--check if root node
-if ($rank == $root) {
-
-   try{
-      #--Process arguments and the command line 
-      GetOptions("RM_off|R" => \$OPT{R},
-                 "force|f" => \$OPT{force},
-                 "genome|g=s" => \$OPT{genome},
-                 "cpus|c=i" => \$OPT{cpus},
-                 "predictor=s" =>\$OPT{predictor},
-                 "retry=i" =>\$OPT{retry},
-                 "quiet" =>\$main::quiet,
-                 "CTL" => sub {GI::generate_control_files(); exit(0);},
-                 "help|?" => sub {print $usage; exit(0)}
-                );
-   }
-   catch Error::Simple with{
-      my $E = shift;
-      
-      print STDERR $E->{-text};
-      die "\n\nMaker failed parsing command line options!!\n\n";
-   };
-
-   #varibles that are persistent outside of try
-   my %CTL_OPT;
@@ -383 +382 @@
-      if($t_need_flag > 0){
-         my $tier;
-Entry
+Fasta
-            $tier = Process::MpiTiers->new({fasta =>$fasta,
-                                            CTL_OPT => \%CTL_OPT,
@@ -449 +448 @@
-            my $tier;
-
-Entry
+Fasta
-               $tier = Process::MpiTiers->new({fasta => $fasta,
-                                               CTL_OPT => \%CTL_OPT,

bin/gff3_merge

@@ -1 @@
-local
+bin
-
-use strict;
@@ -11 +11 @@
-my $outfile;
-
-
+=s
-            "i=s" => \@files,
-            "o=s" => \$outfile,
@@ -17 +17 @@
-            );
-
-! $datastore && ! @files
+(! $datastore && ! @files) || ! $outfile
-    print $usage;
-    exit();
@@ -26 +26 @@
-    open(IN, "< $datastore");
-
-
-
+
+
+
+
+
+
+
+
+
+
-
-    foreach my $file (@files){
@@ -35 +43 @@
-}
-
+@files = sort @files;
+
+open(my $GFF, "> $outfile");
+print $GFF "\#\#gff-version 3\n";
+
-my ($ANN, $ann_file) = tempfile();
-print $ANN "\#\#gff-version 3\n";
-my ($FAS, $fas_file) = tempfile();
-print $FAS "\#\#FASTA\n";
-
+my %uniq;
-foreach my $file (@files){
-    die "ERROR: The file \'$file\' does not exist\n" if (! -e $file);
@@ -48 +61 @@
-    while (defined(my $line = <IN>)){
-        next if ($line =~ /^\#\#gff-version 3/);
+        if($line =~ /^\#\#genome-build/){
+            next if exists $uniq{$line};
+            $uniq{$line}++;
+            print $GFF $line;
+            next;
+        }
+        if($line =~ /^\#\#sequence-region/){
+            die "ERROR: This contig has already been added\: $line\n" if exists $uniq{$line};
+            $uniq{$line}++;
+            print $GFF $line;
+            next;
+        }
-        if ($line =~ /^\#\#FASTA/){
-            $FH = $FAS;
@@ -56 +81 @@
-    }
-}
-
+
-close($ANN);
-close($FAS);
-
-mv $ann_file
+cat $ann_file >>
-system ("cat $fas_file >> $outfile");
+unlink("$ann_file");
-unlink("$fas_file");

bin/maker

@@ -7 +7 @@
-use lib "$FindBin::Bin/../lib";
-use vars qw($LOG $RANK);
-
+
-BEGIN{
-   $RANK = 'non_mpi';           #default value of rank
@@ -24 +24 @@
-   };
-
+   #supress warnings from storable module
+   $SIG{'__WARN__'} = sub {
+       warn $_[0] if ( $_[0] !~ /Not a CODE reference/ &&
+                       $_[0] !~ /Can\'t store item CODE/
+                       );
+   };
+   
-   #output to log file of seq that caused rank to die
-   $SIG{'__DIE__'} = sub {
-defined 
+
-         my $die_count = $LOG->get_die_count();
-         $die_count++;
@@ -514 +521 @@
-         $pred_gff_keepers = $GFF_DB->phathits_on_chunk($chunk,
-                                                        $q_seq_ref,
-repeat
+pred
-                                                       );
-      }

lib/Bio/Search/Hit/PhatHit/Base.pm

@@ -462 +462 @@
-         $hstr{ $h }++;
-      }
-
-
+
+
+
-      if ($seqType =~ /list|array/i) {
-         return ($qstr, $hstr);

lib/CGL/TranslationMachine.pm

@@ -179 +179 @@
-################################################ subroutine header begin ##
-
+=head2 longest_translation
+
+ Usage     :
+
+=for example
+  $seq = "atgaaaaaauaa";
+
+=for example begin
+
+  use CGL::TranslationMachine;
+  my $t = new CGL::TranslationMachine;
+  my ($longest_orf,$offset) = $t->longest_translation($seq);
+
+=for example end
+
+=for example_testing
+  is($longest_orf, "MKK", "Did it find the right orf?");
+  is($offset, 0, "Did it get the offset right?");
+
+ Purpose   : Find the longest open reading frame in any of three
+             frames in the input sequence.
+ Returns   : A list of the amino acid sequence of the longest
+             open reading frame and the offset at which it begins in
+             the input sequence.
+ Arguments : A translatable sequence.
+ Throws    :
+ Comments  :
+           :
+ See Also  : Bio::Tools::CodonTable
+
+=cut
+
+################################################## subroutine header end ##
+
+sub longest_translation_plus_stop {
+    my $self = shift;
+    my $seq  = shift;
+    my $seqlength = length($seq);
+
+    my ($longest,$longest_len,$offset,@frames) = ('',0,0);
+    for( my $i = 0; $i < FRAMES; $i++ ) {
+        my $counter = $i;
+        foreach my $orf ( split(/\*/,$self->translate(substr($seq,$i))) ) {
+            my $orflen = length($orf);
+            if( $orflen > $longest_len ) {
+                $longest = $orf;
+                $offset  = $counter;
+                $longest_len = $orflen;
+                my $c_start = ($orflen * 3) + $offset;
+                my $codon = ($c_start < $seqlength) ? substr($seq, $c_start, 3) : '';
+                if(length($codon) == 3 && $self->translate($codon) eq '*'){
+                    $longest .= '*';
+                }
+            }
+            # multi by 3 for aa->codon
+            # 1 for the stop codon (included
+            # in seq but not in length) 
+            $counter += ($orflen + 1)*3;
+        }
+    }
+
+    return ($longest,$offset);
+}
+
+################################################ subroutine header begin ##
+
-=head2 translate_from_offset
-

lib/Dumper/GFF/GFFV3.pm

@@ -10 +10 @@
-use PostData;
-use PhatHit_utils;
+use File::Copy;
-
-@ISA = qw(
@@ -113 +114 @@
-
-    system("cat $seq_f >> $ann_f");
-"$seq_f"
-system("mv $ann_f $gff_f"
+$seq_f
+move($ann_f, $gff_f
-
-    die "ERROR: GFF3 file not created\n" if(! -e $gff_f);
@@ -440 +441 @@
-        my $type;
-        if    ($class =~ /^blastx$/i){
-
-
-
-
-
+
+
+
-                $type = $k eq 'hit' ? 'protein_match' : 'match_part'; 
-        }
-
-
+
+
+
+
+
+
+
+
+
+
+
-        }
-        elsif ($class =~ /^blastn$/i){
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-        }
-        else {
-       
+    
-        }
-
@@ -686 +701 @@
-        push(@data, $seq_id, 'maker', 'mRNA', $t_b, $t_e, '.', $t_s, '.');
-        my $nine = 'ID='.$t_id.';Parent='.$g_id.';Name='.$t_name;
-;aed='.$AED.'eval_score='.$score
+'
-           $nine .= ';'.$t_hit->{-attrib} if($t_hit->{-attrib});
-        push(@data, $nine);

lib/Fasta.pm

@@ -149 +149 @@
-
-        my $def = Fasta::getDef($fasta_ref);
-D
+d
-
-        return $seq_id;

lib/GFFDB.pm

@@ -371 +371 @@
-
-    $data[1] = "$tag:$data[1]" if($tag && $data[1] !~ /^$tag\:/);
+    $data[6] = '+' if($data[6] =~ /^\.$|^0$/); #fixes some repeat entries
-
-    my %l = (seqid  => $data[0],
@@ -443 +444 @@
-    }
-    elsif($h_type eq 'pred'){
-
+
+
+
-    }
-    elsif($h_type eq 'other'){
@@ -600 +603 @@
-
-        my $tran_id   = $t->{id};
-
-
+
+
+
-
-        my $description = "g_name=$gene_name;g_id=$gene_id;t_name=$tran_name;t_id=$tran_id" unless(! $gene_name);
@@ -678 +682 @@
-    my $hit_strand = 1;
-    my $hit_name = $t->{name};
-
-
+
+
+
-    foreach my $e (@{$t->{cdss}}){
-        my @args;
@@ -786 +791 @@
-    my $hit_strand = 1;
-    my $hit_name = $t->{name};
-$hit_name =~ s/\s*QI\:([\d\|\.-]+)$/
-
+($hit_name) = $hit_name =~ /(^[^\s]+)
+
-
-    foreach my $e (@{$t->{exons}}){
@@ -907 +912 @@
-    my $cdss  = _grab(['CDS'],  $features);
-    my $mRNAs = _grab(['mRNA'], $features);
+    my $UTRs = _grab(['five_prime_UTR', 'three_prime_UTR'], $features);
-
-    foreach my $p_id (keys %{$mRNAs}){
@@ -926 +932 @@
-            my $id=_get_annotation($f,'ID');
-            my $name=_get_annotation($f,'Name');
+
+            #take care of wormbases incorrect parentage
+            if(exists $exons->{$id}){
+                foreach my $mRNA (@{$mRNAs->{$id}}){
+                    my $t_id = $mRNA->{id};
+                    $mRNA->{exons} = _fix_wormbase($exons->{$id}, $cdss->{$t_id}, $UTRs->{$t_id});
+                }
+            }
-            
-            push(@genes, {'f'       => $f,
@@ -949 +963 @@
-}
-#-------------------------------------------------------------------------------
+#try and discover the exon parantage for wormbase entries
+sub _fix_wormbase {
+    my $exons = shift || [];
+    my $cdss = shift || [];
+    my $UTRs = shift || [];
+
+    my @keepers;
+    foreach my $exon (@$exons){
+        my $e = $exon->{f};
+        my $eB = $e->start;
+        my $eE = $e->end;
+
+        my $ok = 0;
+        my $okB = 0;
+        my $okE = 0;
+
+        #check if exon goes with this CDS
+        foreach my $piece (@$cdss){
+            my $p = $piece->{f};
+            my $pB = $p->start;
+            my $pE = $p->end;
+
+            if($pB == $eB && $pE == $eE){
+                $ok = 1;
+            }
+            elsif($pB == $eB && $pE != $eE){
+                $okB = 1;
+            }
+            elsif($pB != $eB && $pE == $eE){
+                $okE = 1;
+            }
+
+            last if($ok);
+        }
+
+        if($ok){
+            push(@keepers, $exon);
+            next;
+        }
+
+        #check if exon goes with or is completed by this UTR
+        foreach my $piece (@$UTRs){
+            my $p = $piece->{f};
+            my $pB = $p->start;
+            my $pE = $p->end;
+
+            if($pB == $eB && $pE == $eE){
+                $ok = 1;
+            }
+            elsif($pB == $eB && $pE != $eE){
+                $okB = 1;
+                $ok = 1 if($okE);
+            }
+            elsif($pB != $eB && $pE == $eE){
+                $okE = 1;
+                $ok = 1 if($okB);
+            }
+
+            last if($ok);
+        }
+
+        if($ok){
+            push(@keepers, $exon);
+            next;
+        }
+    }
+
+    return \@keepers;
+}
+#-------------------------------------------------------------------------------
-#try too build a sructure similar to that producd by _get_genes
-#but for non gene hits in the gff3
@@ -959 +1043 @@
-    foreach my $f (@{$features}){
-        my $tag_t = $f->primary_tag();
+
+        next if($tag_t =~ /^gene$|^mRNA$|^exon$|^CDS$/);
+
-        my $id = _get_annotation($f,'ID');
-        my $name = _get_annotation($f, 'Name'); 
@@ -1266 +1353 @@
-
-    my @sorted;
-
+
+
-        @sorted = 
-            sort {$a->{f}->start <=> $b->{f}->start} @{$t->{exons}}; 

lib/GI.pm

@@ -13 +13 @@
-use URI::Escape;
-use File::Path;
+use File::Copy;
-use Data::Dumper;
-use Getopt::Long;
@@ -1046 +1047 @@
-sub build_fasta_index {
-   my $db = shift;
-   print "\n\n".$db."\n\n";
-   my $index = new Bio::DB::Fasta($db);
-   return $index;
@@ -1069 +1069 @@
-
-
-$
+
-      if (($type eq 'blastn' && ! -e $file.'.xnd') ||
-          ($type eq 'blastx' && ! -e $file.'.xpd') ||
@@ -1083 +1083 @@
-      }
-   }
-$
+
-      if (($type eq 'blastn' && ! -e $file.'.nsq') ||
-          ($type eq 'blastx' && ! -e $file.'.psq') ||
@@ -1098 +1098 @@
-   }
-   else {
-
+ not \'$command\'
-   }
-}
@@ -1142 +1142 @@
-   #copy db to local tmp dir and run xdformat or formatdb 
-   if (! @{[<$tmp_db.xn?*>]} && (! -e $blast_finished || $opt_f) ) {
-system("cp $db $tmp_db"
+copy($db, $tmp_db
-      dbformat($formater, $tmp_db, 'blastn');
-   }
@@ -1391 +1391 @@
-   #copy db to local tmp dir and run xdformat or format db 
-   if (! @{[<$tmp_db.xp?*>]} && (! -e $blast_finished || $opt_f) ) {
-system("cp $db $tmp_db"
-
+ copy($db, $tmp_db
+ 
-   }
-   elsif (-e $blast_finished && ! $opt_f) {
@@ -1652 +1652 @@
-   #copy db to local tmp dir and run xdformat or formatdb
-   if (! @{[<$tmp_db.xn?*>]} && (! -e $blast_finished || $opt_f) ) {
-system("cp $db $tmp_db"
+copy($db, $tmp_db
-      dbformat($formater, $tmp_db, 'tblastx');
-   }
@@ -2495 +2495 @@
-   print OUT "fgenesh:$O{fgenesh} #location of fgenesh executable\n";
-#   print OUT "twinscan:$O{twinscan} #location of twinscan executable\n";
-
+#
-   print OUT "\n";
-   print OUT "#-----Other Algorithms\n";
@@ -2514 +2514 @@
-   print OUT "eva_hspmax:$O{eva_hspmax}\n";
-   print OUT "eva_gspmax:$O{eva_gspmax}\n";
-
+#
-   close (OUT);
-}

lib/Process/MpiChunk.pm

@@ -5 +5 @@
-use strict;
-
-
+ :warndie
-use Error::Simple;
-use Storable;
@@ -623 +623 @@
-                        q_def
-                        masked_total_seq
+                        fasta
-                        LOG
-                        CTL_OPT}
@@ -633 +634 @@
-            my $q_def = $VARS->{q_def};
-            my $masked_total_seq = $VARS->{masked_total_seq};
+            my $fasta = $VARS->{fasta};
-            my $the_void = $VARS->{the_void};
-            my $safe_seq_id = $VARS->{safe_seq_id};
@@ -640 +642 @@
-            my $masked_fasta = Fasta::toFasta($q_def.' masked', \$masked_total_seq);
-            my $masked_file = $the_void."/query.masked.fasta";
- ,
+, 
-         
+            my $unmasked_file = $the_void."/query.fasta";
+            FastaFile::writeFile(\$fasta, $unmasked_file);
+
-            #==ab initio predictions here
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-            #==QRNA noncoding RNA prediction here
-            my $qra_preds = [];
@@ -1156 +1170 @@
-               $pred_gff_keepers = $GFF_DB->phathits_on_chunk($chunk,
-                                                              $q_seq_ref,
-repeat
+pred
-                                                             );
-            }
@@ -1669 +1683 @@
-            
-            #-- write to DS log the finished files
-
+#
-            
-            #--- clear the log variable

lib/Widget/augustus.pm

@@ -359 +359 @@
-                        $this_gene{score} = $score;
-                }
-
+els
-                        my $strand = $fields[6] eq '+' ? 1 : -1;
-
@@ -394 +394 @@
-                                                  };
-
-                }
-                elsif ($fields[1] eq 'AUGUSTUS' && $fields[2]  eq 'gene'){
-                }
-                elsif ($fields[1] eq 'AUGUSTUS' && $fields[2]  eq 'transcript'){
@@ -456 +454 @@
-                my @stuff = split(/\n/, $line);
-
+                my $t_count = grep {/^[^\#]/ && /^[^\s]+\s+AUGUSTUS\s+transcript\s+/} @stuff;
+
+                #augustus can call multiple transcripts for inconsistent hints
+                if($t_count > 1){
+                    my @multi_stuff;
+                    my $count;
+                    foreach my $l (@stuff){
+                        if($l =~ /^\#/){ #maybe important someday
+                            foreach my $s (@multi_stuff){
+                                push(@$s, $l);
+                            }
+                        }
+                        elsif($l =~ /^[^\s]+\s+AUGUSTUS\s+gene\s+/){
+                            foreach my $s (@multi_stuff){
+                                push(@$s, $l);
+                            }
+                        }
+                        elsif($l =~ /^[^\s]+\s+AUGUSTUS\s+transcript\s+/){
+                            $count++;
+                            push(@{$multi_stuff[$count-1]}, $l);
+                        }
+                        elsif($l =~ /^[^\s]+\s+AUGUSTUS\s+[^\s]+\s+/){
+                            push(@{$multi_stuff[$count-1]}, $l);
+                        }
+                    }
+
+                    foreach my $s (@multi_stuff){
+                        my $gene = parse_gene($s);
+                        push(@genes, $gene);
+                    }
+                    next;
+                }
+
-                my $gene = parse_gene(\@stuff);
-                push(@genes, $gene);

lib/clean.pm

@@ -19 +19 @@
-#------------------------------------------------------------------------
-#--------------------------- FUNCTIONS ----------------------------------
+#------------------------------------------------------------------------
+sub purge_short_single_exons {
+    my $phat_hits = shift;
+    my $min = shift;
+
+    my @keepers;
+    foreach my $hit (@{$phat_hits}){
+        if ($hit->num_hsps == 1){
+            push(@keepers, $hit) if($hit->{_hsps}->[0]->length('query') >= $min);
+        }
+        else{
+            push(@keepers, $hit);
+        }
+    }
+    
+    return \@keepers
+}
-#------------------------------------------------------------------------
-sub throw_out_bad_splicers {
@@ -58 +75 @@
-sub alt_sort {
-        $b->num_hsps <=> $a->num_hsps || $b->length  <=>  $a->length;
+}
+#------------------------------------------------------------------------
+sub is_identical_form {
+    my $hit1 = shift;
+    my $hit2 = shift;
+
+    my @hsps1 = sort {$a->start('query') <=> $b->start('query')} $hit1->hsps;
+    my @hsps2 = sort {$a->start('query') <=> $b->start('query')} $hit2->hsps;
+
+    return if (@hsps1 != @hsps2);
+
+    for(my $i = 0; $i < @hsps1; $i++){
+        return if($hsps1[$i]->start('query') != $hsps2[$i]->start('query'));
+        return if($hsps1[$i]->end('query') != $hsps2[$i]->end('query'));
+    }
+
+    return 1;
-}
-#------------------------------------------------------------------------

lib/cluster.pm

@@ -443 +443 @@
-    my $i = 0;
-    foreach my $c (keys %{$map}){
+        next if(! defined $map->{$c});
-        foreach my $m (@{$map->{$c}}){
-            my $hit = $lookup{$m};

lib/evaluator/funs.pm

@@ -493 +493 @@
-
-        my $seq = $box->{transcription_seq};
+        my $length = length($$seq);
+
-        my $first_codon = uc(substr($$seq, $box->{translation_offset}, 3));
-        my $last_codon  = uc(substr($$seq, $box->{translation_end}-4, 3));

lib/evaluator/quality_index.pm

@@ -258 +258 @@
-
-        my $l_trans = get_length_protein($t);
-
+
-        my ($length_5, $length_3) = get_length_of_utrs($t);
-

lib/maker/auto_annotator.pm

@@ -24 +24 @@
-use evaluator::AED;
-use shadow_AED;
+use Storable qw(dclone);
+
+$Storable::forgive_me = 1; 
-
-@ISA = qw(
@@ -44 +47 @@
-        #only ESTs from splice concious algorithms i.e. no blastn
-        my $clean_est = get_selected_types($est_hits,'est2genome', 'est_gff');
+        my $clean_altest = $alt_est_hits;
-        if ($single_exon == 1) {
-$clean_est = $clean_est;
+#do nothing
-        }else {
-            # don't use unpliced single exon ESTs-- usually genomic contamination
-            $clean_est = clean::purge_single_exon_hits($clean_est);
+            $clean_altest = clean::purge_single_exon_hits($clean_altest);
-        }
-
-        # throw out the exonerate est hits with weird splice sites
-        $clean_est = clean::throw_out_bad_splicers($clean_est, $seq);
-
-        
-        # combine puts type in order they are given
@@ -63 +67 @@
-                            $prot_hits,
-                            $clean_est, 
-alt_est_hits
+clean_altest
-                           );
-
@@ -70 +74 @@
-        my $p_clusters = cluster::shadow_cluster(0, $seq, $p, 10);
-        my $m_clusters = cluster::shadow_cluster(0, $seq, $m, 10);
+
+        #purge after clustering so as to still have the effect of evidence joining ESTs
+        $p_clusters = purge_short_ESTs_in_clusters($p_clusters, 250);
+        $m_clusters = purge_short_ESTs_in_clusters($m_clusters, 250);
+
-        my $careful_clusters = [];
-        push(@{$careful_clusters}, @{$p_clusters});
-        push(@{$careful_clusters}, @{$m_clusters});
-
+
-        # identify the ab-inits that fall within and between clusters
-        my ($c_index, $hit_one, $hit_none, $hit_mult) = segment_preds($predictions, 
@@ -83 +92 @@
-        }
-
+        #--new clusters built by merging clusters that overlap preds
+        #--used for ab-initio coring
+        #--must be ran before merge as merge alters cluster content
+        my $pred_clusters = join_clusters_on_pred($predictions, $careful_clusters, $c_index);
+
+        #--split preds across clusters that they overlap
-        merge_into_cluster($hit_one, $careful_clusters, $c_index);
-        merge_into_cluster($hit_mult, $careful_clusters, $c_index);
-
+        #--data prep
-        my $c_id = 0;
-        my @bx_data;
@@ -100 +116 @@
-           $c_id++;
-        }
-        #--- end c_bag
-
-        #--- processing for scoring ab-initio predictions
-        my @pr_data;
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
-
-        return (\@bx_data, \@gf_data, \@pr_data);
@@ -229 +228 @@
-      }
-   }
+}
+#------------------------------------------------------------------------
+sub join_clusters_on_pred {
+   my $hits = shift;
+   my $clusters = shift;
+   my $c_index = shift;
+
+   my @p_clusters;
+
+   foreach my $s (@{$hits}, ){
+      my @keys = keys %{$c_index->{$s->{temp_id}}};
+      my @new_c = ($s);
+      foreach my $i (@keys){
+         push(@new_c, @{$clusters->[$i]});
+      }
+      push(@p_clusters, \@new_c);
+   }
+
+   return \@p_clusters;
-}
-#------------------------------------------------------------------------
@@ -252 +270 @@
-# }
-#------------------------------------------------------------------------
+sub purge_short_ESTs_in_clusters{
+    my $clusters = shift;
+    my $min = shift;
+    
+    my @c_keepers;
+    foreach my $c (@$clusters){
+        my $ests_in_cluster  = get_selected_types($c,'est2genome', 'est_gff');
+        my $ps_in_cluster    = get_selected_types($c,'protein2genome');
+        my $bx_in_cluster    = get_selected_types($c,'blastx', 'protein_gff');
+        my $alt_ests_in_cluster = get_selected_types($c,'tblastx', 'altest_gff');
+        my $models_in_cluster = get_selected_types($c,'model_gff', 'maker');
+        my $preds_in_cluster = get_selected_types($c,'snap', 'augustus', 'fgenesh', 'twinscan', 'pred_gff');
+
+        $ests_in_cluster = clean::purge_short_single_exons($ests_in_cluster, $min);
+        $alt_ests_in_cluster = clean::purge_short_single_exons($alt_ests_in_cluster, $min);
+
+        my @new_c;
+        push(@new_c, @$ests_in_cluster);
+        push(@new_c, @$alt_ests_in_cluster);
+        push(@new_c, @$ps_in_cluster);
+        push(@new_c, @$bx_in_cluster);
+        push(@new_c, @$models_in_cluster);
+
+        if(@new_c){
+            push(@new_c, @$preds_in_cluster);
+            push(@c_keepers, \@new_c);
+        }
+    }
+
+    return \@c_keepers;
+}
+#------------------------------------------------------------------------
-#returns an array of hashes with the following atributes
-#ests => set of all ests
@@ -265 +315 @@
-        my $ps_in_cluster    = get_selected_types($c,'protein2genome');
-        my $bx_in_cluster    = get_selected_types($c,'blastx', 'protein_gff');
-
+, 'altest_gff'
-        my $models_in_cluster = get_selected_types($c,'model_gff', 'maker');
-        my $preds_in_cluster = get_selected_types($c,'snap', 'augustus', 'fgenesh', 'twinscan', 'pred_gff');
-
-s
+
-
-        # groups of most informative protein hits
@@ -327 +377 @@
-        my $ps_in_cluster    = get_selected_types($c,'protein2genome');
-        my $bx_in_cluster    = get_selected_types($c,'blastx', 'protein_gff');
-
+, 'altest_gff'
-        my $preds_in_cluster = get_selected_types($c,'snap', 'augustus', 'fgenesh', 'twinscan', 'pred_gff');
-
-s
+
-
-        # groups of most informative protein hits
@@ -368 +418 @@
-        my $ps_in_cluster    = get_selected_types($c,'protein2genome');
-        my $bx_in_cluster    = get_selected_types($c,'blastx', 'protein_gff');
-
+, 'altest_gff'
-        my $preds_in_cluster = get_selected_types($c,'snap', 'augustus', 'fgenesh', 'twinscan', 'pred_gff');
-
@@ -514 +564 @@
-        }
-
-
+hint based 
-        foreach my $prdr (@{$CTL_OPTIONS->{_predictor}}){
-
+ || $prdr eq 'abinit'
-           my $transcripts = run_it($bx_data,
-                                    $the_void,
@@ -538 +588 @@
-
-           $annotations{$prdr} =  $annot;
-           push(@bag, @{$annot});
-        }
-
@@ -562 +611 @@
-                                      );
-
-
-
-
-
-
+
-
-        return \%annotations;
+}
+#------------------------------------------------------------------------
+sub get_n_count{
+    my $g = shift;
+
+    my ($start, $end) = ($g->{g_start}, $g->{g_end});
+    my $offset = $start;
+    my $length = $end - $start + 1;
+    my @b_seq = map {0} (1..$length);
+
+    foreach my $s (@{$g->{t_structs}}){
+        my $hit = $s->{hit}; 
+        my @hsps = $hit->hsps() if defined $hit->hsps();
+        foreach my $hsp (@hsps){
+            my $s = $hsp->start('query') - $offset;
+            my $e = $hsp->end('query') - $offset;
+            
+            #array space coors
+            die "ERROR: Start value not permited!!\n" if($s >= $length || $s < 0);
+            die "ERROR: End value not permited!!\n" if($e < 0 || $e >= $length);
+            
+            @b_seq[$s..$e] = map {1} ($s..$e);
+        }
+    }
+
+    my @index = (0, 0);
+
+    foreach my $i (@b_seq){
+        $index[$i]++;
+    }
+
+    return $index[1];
-}
-#------------------------------------------------------------------------
@@ -579 +656 @@
-   my $CTL_OPTIONS = shift;
-
+   return $annotations->{gff} 
+      if(@{$CTL_OPTIONS->{_predictor}} == 1 && grep {/^gff$/} @{$CTL_OPTIONS->{_predictor}});
+
-   my %order; #this is important for sorting
-
-
+
+
+
+
-   my $i = @{$CTL_OPTIONS->{_predictor}};
-   foreach my $p (@{$CTL_OPTIONS->{_predictor}}){
@@ -588 +670 @@
-      foreach my $g (@{$annotations->{$p}}){
-         if($g->{g_strand} == 1){
-
+
+
+
-         }
-         elsif($g->{g_strand} == -1){
-
+
+
+
-         }
-         else{
@@ -600 +686 @@
-   %SCORE = %order; #must be set for sort to work correctly
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
-       my $g_B = $g->{g_start};
-       my $g_E = $g->{g_end};
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-       }
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-   #write evaluator reports
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-}
-#------------------------------------------------------------------------
@@ -709 +850 @@
-         next if(! defined $model);
-         my $transcript = $model;
+
+         my $utr_trans; 
+         if($predictor eq 'abinit' && defined($ests)){
+             my $copy = dclone($transcript);
+             $utr_trans = pneu($ests, $copy, $seq);
+             push(@transcripts, [$utr_trans, $set, $transcript]);
+             if (! clean::is_identical_form($utr_trans, $transcript)){
+                 push(@transcripts, [$utr_trans, $set, undef]);
+             }
+         }
-         push(@transcripts, [$transcript, $set, undef]);
+
-         $i++;
-         next;
@@ -717 +869 @@
-         next if (! defined $mia);
-         my $transcript = pneu($ests, $mia, $seq);
-
+ if $transcript
-         $i++;
-         next;
-          
+
-
-      my ($pred_shots, $strand) = get_pred_shot($seq, 
@@ -850 +1002 @@
-        
-        my $t_name = "$g_name-mRNA-$i"; #affects GFFV3.pm
-        $f->name($t_name);
-
-        #my $qi = maker::quality_index::get_transcript_qi($f,$evi,$offset,$len_3_utr,$l_trans);
-
-        #----evaluator here
@@ -858 +1007 @@
-        my $pol_e_hits  = get_selected_types($evi->{ests}, 'est2genome', 'est_gff');
-        my $blastx_hits = get_selected_types($evi->{gomiph},'blastx', 'protein_gff');
+        my $tblastx_hits = get_selected_types($evi->{alt_ests},'tblastx', 'altest_gff');
-        my $abinits = $evi->{all_preds};
-
-        my @bag = (@$pol_p_hits,
-                   @$pol_e_hits,
-
+
+
-                  );
-
-        #holds evalutor struct
-
-
-
-
-
-
-
-
-
-
-
-
-
+#
+#
+#
+#
+#
+#
+#
+#
+#
+#
+#
+#
+#
-        
-
-
-
+
+
+
+
+
-        #----
-        $t_name .= " AED:";
-        $t_name .= sprintf '%.2f', $AED; # two decimal places
-        $t_name .= " $qi";
+        $f->name($t_name); #give name to hit
-
-        my $t_struct = {'hit'      => $f,
@@ -897 +1051 @@
-                        'AED'      => $AED,
-                        'score'    => $score,
-eva->{report}
+report
-                    };
-
@@ -1069 +1223 @@
-      }
-      
-
-
-
+'';#
+'';#
+'';#
-      #----
-      
-      my $score = 0;
+      my $AED = 1;
-      my $i = 1;
-      foreach my $f (@{$c}) {
@@ -1082 +1237 @@
-         load_transcript_struct($f, $g_name, $i, $seq, $evidence, $so_code,
-                                $geneAED, $alt_spli_sup, $the_void, $CTL_OPTIONS);
+
+         if(!$f || ! $t_struct){
+             sleep 1;
+         }
+
-         push(@t_structs, $t_struct);
-         $score = $t_struct->{score} if($t_struct->{score} > $score);
+         $AED = $t_struct->{AED} if($t_struct->{AED} < $AED);
-         $i++;
-      }
@@ -1095 +1256 @@
-                         'g_strand'  => $g_strand,
-                         'score'     => $score,
-
+AED,#$
-                         'predictor' => $predictor,
-                         'so_code'   => $so_code
@@ -1151 +1312 @@
-                }
-        }
+        if(! @exons){
+            return;
+        }
-        my @sorted_b = sort {$a->start('query') <=> $b->start('query')} @exons;
-        my @sorted_e = sort {$b->end('query')   <=> $a->end('query')}   @exons;
+
+        my $ref = ref($sorted_b[0]);
-
-        my $g_start = $sorted_b[0]->start('query');
@@ -1227 +1393 @@
-        my $best = shift(@sorted);
-
-
-
+
-}
-#------------------------------------------------------------------------
@@ -1259 +1424 @@
-        my $p_seq_2 = $tM->translate_from_offset($seq, $offset);
-
-)\*?
+\*?)
-
-        die "logic error in auto_annotate::get_off_and_str!\n"
@@ -1272 +1437 @@
-        my $tM = new CGL::TranslationMachine();
-
-
+
+
+
-
-        $p_seq =~ s/\*$//; # added 11/21/06
-
-        my $end = length($p_seq)*3 + $offset + 1 + 3;
-
-        my ($trim, $leader);
-        if ($p_seq =~ /^M/){
-                # easy  it begins with an M....
@@ -1287 +1451 @@
-                my ($off_new, $p_seq_new) = get_longest_m_seq($seq);
-
- + 3
+
-                            if defined($p_seq_new);
-
+                $p_seq_new =~ s/\*$// if($p_seq_new); # added 11/21/06
-
-                if (!defined($p_seq_new)){
@@ -1295 +1460 @@
-                }
-                elsif ($offset < 3 && length($p_seq) - length($p_seq_new) > 30){
-               
+
-                }
-                else {

lib/maker/join.pm

@@ -348 +348 @@
-        my $length = 0;
-        foreach my $hsp (@anno_hsps){
-
-
+
+
+
+
-        }
-

lib/shadow_AED.pm

@@ -9 +9 @@
-   my $tran = shift;
-
-0.5
+1
-
-   my ($start, $end) = ($tran->start('query'), $tran->end('query'));
@@ -68 +68 @@
-   my $spec = $index{3}/($index{2} + $index{3}); #specificity
-   my $sens = $index{3}/($index{1} + $index{3}); #sensitivity
-
+1 - 
-
-   return $AED;