<commit_before>#!/usr/bin/perl -w
use strict;
use warnings;
use CGI;
use CGI::Carp 'fatalsToBrowser';
use CGI::Cookie;
use CGI::Ajax;
use HTML::Template;
use Data::Dumper;
use File::Basename;
use File::Temp;
use CoGe::Accessory::GenBank;
use CoGe::Accessory::LogUser;
use CoGe::Accessory::Web;
use CoGe::Accessory::Restricted_orgs;
use CoGe::Accessory::bl2seq_report;
use CoGe::Accessory::blastz_report;
use CoGe::Accessory::lagan_report;
use CoGe::Accessory::chaos_report;
use CoGe::Accessory::dialign_report;
use CoGe::Accessory::dialign_report::anchors;
use CoGe::Graphics;
use CoGe::Graphics::Chromosome;
use CoGe::Graphics::Feature;
use CoGe::Graphics::Feature::Gene;
use CoGe::Graphics::Feature::NucTide;
use CoGe::Graphics::Feature::GAGA;
use CoGe::Graphics::Feature::Exon_motifs;
use CoGe::Graphics::Feature::AminoAcid;
use CoGe::Graphics::Feature::Domain;
use CoGe::Graphics::Feature::HSP;
use CoGeX;
use CoGeX::Feature;
use DBIxProfiler;
use DBI;
use LWP::Simple;
use Parallel::ForkManager;
use Statistics::Basic;
use Benchmark qw(:all);

# for security purposes

$ENV{PATH} = "/opt/apache/CoGe/";
delete @ENV{ 'IFS', 'CDPATH', 'ENV', 'BASH_ENV' };
#for chaos
$ENV{'LAGAN_DIR'} = '/opt/apache/CoGe/bin/lagan/';
#for dialign
$ENV{'DIALIGN2_DIR'} = '/opt/apache/CoGe/bin/dialign2_dir/';
use vars qw( $PAGE_NAME $DATE $DEBUG $BL2SEQ $BLASTZ $LAGAN $CHAOS $DIALIGN $TEMPDIR $TEMPURL $USER $FORM $cogeweb $BENCHMARK $coge $NUM_SEQS $MAX_SEQS $REPEATMASKER $MAX_PROC);
$PAGE_NAME = "GEvo.pl";
$BL2SEQ = "/usr/bin/bl2seq ";
$BLASTZ = "/usr/bin/blastz ";
$LAGAN = "/opt/apache/CoGe/bin/lagan/lagan.pl";
$CHAOS = "/opt/apache/CoGe/bin/lagan/chaos_coge";
$DIALIGN = "/opt/apache/CoGe/bin/dialign2_dir/dialign2-2_coge";
$REPEATMASKER = "/opt/apache/CoGe/bin/RepeatMasker/RepeatMasker";
$TEMPDIR = "/opt/apache/CoGe/tmp/GEvo";
$TEMPURL = "/CoGe/tmp/GEvo";
$MAX_PROC=8;
# set this to 1 to print verbose messages to logs
$DEBUG = 0;
$BENCHMARK = 0;
$NUM_SEQS = 3;
$MAX_SEQS = 21;
$| = 1; # turn off buffering 
$DATE = sprintf( "%04d-%02d-%02d %02d:%02d:%02d",
		 sub { ($_[5]+1900, $_[4]+1, $_[3]),$_[2],$_[1],$_[0] }->(localtime));

$FORM = new CGI;
$CGI::POST_MAX= 60 * 1024 * 1024; # 24MB
$CGI::DISABLE_UPLOADS = 0; 
($USER) = CoGe::Accessory::LogUser->get_user();

$coge = CoGeX->dbconnect();
#$coge->storage->debugobj(new DBIxProfiler());
#$coge->storage->debug(1);

my %ajax = CoGe::Accessory::Web::ajax_func();
$ajax{dataset_search} = \&dataset_search; #override this method from Accessory::Web for restricted organisms
$ajax{feat_search} = \&feat_search; 
my $pj = new CGI::Ajax(
		       run=>\&run,
		       loading=>\&loading,
		       add_seq=>\&add_seq,
		       get_file=>\&get_file,
		       gen_go_run=>\&gen_go_run,
		       gen_hsp_colors =>\&gen_hsp_colors,
		       save_settings_gevo=>\&save_settings_gevo,
		       reset_settings_gevo=>\&reset_settings_gevo,
		       %ajax,
		      );
$pj->JSDEBUG(0);
$pj->DEBUG(0);
#$pj->js_encode_function('escape');
print $pj->build_html($FORM, \&gen_html);

#$USER="elyons";print gen_html();

sub loading
  {
    my $message = shift || "Generating results. . .";
    return qq{<div class="dna"><div id="loading">$message</div></div>}; 
  }

sub gen_html
  {
    my $html;# =  "Content-Type: text/html\n\n";
    unless ($USER)
      {
	$html = login();
      }
    else
      {
	my $template = HTML::Template->new(filename=>'/opt/apache/CoGe/tmpl/generic_page.tmpl');
	$template->param(LOGO_PNG=>"GEvo-logo.png");
	$template->param(TITLE=>'Genome Evolution Analysis');
	$template->param(HELP=>'GEvo');
	$template->param(USER=>$USER);
	$template->param(DATE=>$DATE);
	$template->param(NO_BOX=>1);
	$template->param(BODY=>gen_body());
	my $prebox = HTML::Template->new(filename=>'/opt/apache/CoGe/tmpl/GEvo.tmpl');
	$prebox->param(RESULTS_DIV=>1);
	$template->param(PREBOX=>$prebox->output);
	
	$html .= $template->output;
      }
    return $html;
  }

sub gen_body
  {
    my $form = $FORM;
    my $prefs = load_settings(user=>$USER, page=>$PAGE_NAME);
    my $num_seqs = get_opt(params=>$prefs, form=>$form, param=>'num_seqs');
    $num_seqs = $NUM_SEQS unless defined $num_seqs;
    $MAX_SEQS = 10 if $form->param('override');
    my $message;
    if (! ($num_seqs =~ /^\d+$/) )
      {
	$message .= "Problem with requested number of sequences: '$num_seqs'.  Defaulting to $NUM_SEQS input sequences.";
	$num_seqs = $NUM_SEQS;
      }
    elsif ($num_seqs < 2)
      {
	$message .= "Minimum number of sequences to compare is two.";
	$num_seqs = 2;
      }
    elsif ($num_seqs > $MAX_SEQS)
      {
	$message .= "Maximum number of sequence is $MAX_SEQS.";
	$num_seqs = $MAX_SEQS;
      }
    my @seq_nums;
    my @seq_sub;
    my $autosearch_string;
    for (my $i = 1; $i <= $num_seqs; $i++)
      {
	my $draccn;
	$draccn= $form->param("accn".$i) if $form->param("accn".$i);
	my $pos = $form->param("x".$i) if $form->param("x".$i);
	my $chr = $form->param("chr".$i) if $form->param("chr".$i);
	my $drfid = $form->param("fid".$i) if $form->param("fid".$i);
	if ($drfid && !$draccn)
	  {
	    my $feat = $coge->resultset('Feature')->find($drfid);
	    ($draccn) = $feat->primary_name if $feat;
	  }
	my $drup = $form->param('dr'.$i.'up') if defined $form->param('dr'.$i.'up');
	my $drdown = $form->param('dr'.$i.'down') if defined $form->param('dr'.$i.'down');
	$drup = $form->param('drup'.$i) if defined $form->param('drup'.$i);
	$drdown = $form->param('drdown'.$i) if defined $form->param('drdown'.$i);
	$drup = 10000 unless defined $drup;
	$drdown = 10000 unless defined $drdown;
#	$drup += $pad_gs if $pad_gs;
#	$drdown += $pad_gs if $pad_gs;
	my $dsid = $form->param('dsid'.$i) if $form->param('dsid'.$i);
	my $gbaccn = $form->param("gbaccn".$i) if $form->param("gbaccn".$i);
	my $gbstart = $form->param("gbstart".$i) if $form->param("gbstart".$i);
	$gbstart = 1 unless defined $gbstart;
	my $gblength = $form->param("gblength".$i) if $form->param("gblength".$i);
	my $revy = "checked" if $form->param('rev'.$i);
	my $revn = "checked" unless $revy;
	my $maskexonon = "checked" if $form->param('maskexon'.$i);
	my $maskexonoff = "checked" unless $maskexonon;
	my $refn = "checked" if $form->param('nref'.$i);
	my $refy = "checked" unless $refn;
	$autosearch_string .= 'if ($'.qq!('#accn$i').val()) {dataset_search(['args__accn','accn$i','args__num', 'args__$i'!;
	$autosearch_string .= ", 'args__dsid', 'args__$dsid'" if $dsid;
	$autosearch_string .= ", 'args__featid', 'args__$drfid'" if $drfid;
	$autosearch_string .=qq!],[feat_search_chain]);}!;
	my $dsinfo = $dsid ? qq{<input type="hidden" id="posdsid$i" value="$dsid">} : qq{<input type="hidden" id="posdsid$i">};
	$dsinfo .= $chr ? qq{<input type="hidden" id="chr$i" value="$chr">} : qq{<input type="hidden" id="chr$i">};
	$dsinfo .=get_dataset_info(dsid=>$dsid, chr=>$chr) if $dsid;
	push @seq_nums, {
			  SEQ_NUM=>$i,
			 };
	my %opts = (
		    SEQ_NUM=>$i,
		    REV_YES=>$revy,
		    REV_NO=>$revn,
		    REF_YES=>$refy,
		    REF_NO=>$refn,
		    DRUP=>$drup,
		    DRDOWN=>$drdown,
		    DRACCN=>$draccn,
		    DSID=>$dsid,
		    GBACCN=>$gbaccn,
		    GBSTART=>$gbstart,
		    GBLENGTH=>$gblength,
		    POS=>$pos,
		    DSINFO=>$dsinfo,
		    EXON_MASK_ON=>$maskexonon,
		    EXON_MASK_OFF=>$maskexonoff,

		   );
#	print STDERR "pos $i: $pos\n";
	$opts{COGEPOS} = qq{<option value="cogepos$i" selected="selected">CoGe Database Position</option>} if $pos;
	push @seq_sub, {%opts}
	  
      }
    #page preferences
    my $pad_gs = $form->param("pad_gs") if $form->param("pad_gs");
    $pad_gs = 0 unless $pad_gs;
    my $prog = get_opt(params=>$prefs, form=>$form, param=>'prog');
    $prog = "blastz" unless $prog;
    my $image_width = get_opt(params=>$prefs, form=>$form, param=>'iw');
    $image_width = 1000 unless $image_width;
    my $image_height = get_opt(params=>$prefs, form=>$form, param=>'ih');
    $image_height = 100 unless $image_height;
    my $feature_height = get_opt(params=>$prefs, form=>$form, param=>'fh');
    $feature_height = 10 unless $feature_height;
    my $padding  = get_opt(params=>$prefs, form=>$form, param=>'padding');
    $padding = 2 unless defined $padding;
    my $gc_color = get_opt(params=>$prefs, form=>$form, param=>'gc');
    $gc_color = 0 unless $gc_color;
    my $nt_color = get_opt(params=>$prefs, form=>$form, param=>'nt');
    $nt_color = 1 unless $nt_color;
    my $auto_adjust_feats = get_opt(params=>$prefs, form=>$form, param=>'overlap');
    $auto_adjust_feats = 0 unless defined $auto_adjust_feats;
    my $hiqual = get_opt(params=>$prefs, form=>$form, param=>'hiqual');
    $hiqual = 0 unless $hiqual;
    my $color_hsp = get_opt(params=>$prefs, form=>$form, param=>'colorhsp');
    $color_hsp = 0 unless $color_hsp;
    my $hsp_label = get_opt(params=>$prefs, form=>$form, param=>'hsplabel');
    $hsp_label = undef unless defined $hsp_label;
    my $hsp_limit = get_opt(params=>$prefs, form=>$form, param=>'hsplim');
    $hsp_limit = 0 unless $hsp_limit;
    my $hsp_limit_num = get_opt(params=>$prefs, form=>$form, param=>'hsplimnum');
    $hsp_limit_num = 20 unless defined $hsp_limit_num;
    my $draw_model = get_opt(params=>$prefs, form=>$form, param=>'draw_model');
    $draw_model = "full" unless $draw_model;
    my $hsp_overlap_limit = get_opt(params=>$prefs, form=>$form, param=>'hsp_overlap_limit');
    $hsp_overlap_limit = 0 unless $hsp_overlap_limit;

    my $template = HTML::Template->new(filename=>'/opt/apache/CoGe/tmpl/GEvo.tmpl');
    $template->param(PAD_GS=>$pad_gs);
    $template->param(IMAGE_WIDTH=>$image_width);
    $template->param(IMAGE_HEIGHT=>$image_height);
    $template->param(FEAT_HEIGHT=>$feature_height);
    $template->param(PADDING=>$padding);
    $template->param(HSP_OVERLAP_LIMIT=>$hsp_overlap_limit);
    if ($gc_color) {$template->param(GC_COLOR_YES=>"checked");}
    else {$template->param(GC_COLOR_NO=>"checked");}
    if ($nt_color) {$template->param(NT_COLOR_YES=>"checked");}
    else {$template->param(NT_COLOR_NO=>"checked");}
    if ($auto_adjust_feats) {$template->param(OVERLAP_YES=>"checked");}
    else {$template->param(OVERLAP_NO=>"checked");}
    if ($hiqual) {$template->param(HIQUAL_YES=>"checked");}
    else {$template->param(HIQUAL_NO=>"checked");}
    if ($color_hsp) {$template->param(COLOR_HSP_YES=>"checked");}
    else {$template->param(COLOR_HSP_NO=>"checked");}
    if ($hsp_label eq "staggered") {$template->param(HSP_LABELS_STAG=>"selected");}
    elsif ($hsp_label eq "linear") {$template->param(HSP_LABELS_LIN=>"selected");}
    else {$template->param(HSP_LABELS_NO=>"selected");}
    if ($hsp_limit) {$template->param(HSP_LIMIT_YES=>"checked");}
    else {$template->param(HSP_LIMIT_NO=>"checked");}
    $template->param(HSP_LIMIT_NUM=>$hsp_limit_num);
    if ($draw_model eq "full") {$template->param(DRAW_MODEL_FULL=>"selected");}
    elsif ($draw_model eq "gene") {$template->param(DRAW_MODEL_GENE=>"selected");}
    elsif ($draw_model eq "mRNA") {$template->param(DRAW_MODEL_mRNA=>"selected");}
    elsif ($draw_model eq "CDS") {$template->param(DRAW_MODEL_CDS=>"selected");}
    elsif ($draw_model eq "RNA") {$template->param(DRAW_MODEL_RNA=>"selected");}
    else {$template->param(DRAW_MODEL_NO=>"selected");}

    my $box = HTML::Template->new(filename=>'/opt/apache/CoGe/tmpl/box.tmpl');
    #generate sequence submission selector
    $template->param(SEQ_SELECT=>1);
    $template->param(SEQ_SELECT_LOOP=>\@seq_sub);
    my $seq_submission = $template->output;
    $template->param(SEQ_SELECT=>0);

    #generate the hsp color option
    my $hsp_colors = gen_hsp_colors(num_seqs=>$num_seqs, prefs=>$prefs);
    my $html;
    my $spike_len = spike_filter_select();

    $template->param(SPIKE_LEN=>$spike_len);
    $template->param(SEQ_RETRIEVAL=>1);
    $template->param(NUM_SEQS=>$num_seqs);
    $message .= "<BR/>" if $message;
    $template->param(MESSAGE=>$message);
    $template->param(SEQ_SUB=>$seq_submission);
    $template->param(HSP_COLOR=>$hsp_colors);
    $template->param(GO_RUN=>gen_go_run($num_seqs));
    $template->param(AUTOSEARCH=>$autosearch_string);
    my $gobe_version = `svnversion /opt/apache/CoGe/gobe/flash`;
    $gobe_version =~ s/\n//g;;
    $template->param(GOBE_VERSION=>$gobe_version);
    $box->param(BOX_NAME=>"Options:");
    $template->param(OPTIONS=>1);
    $template->param(ALIGNMENT_PROGRAMS=>algorithm_list($prog));
    $template->param(SAVE_SETTINGS=>gen_save_settings($num_seqs)) unless !$USER || $USER =~ /public/i;
    $box->param(BODY=>$template->output);
    $html .= $box->output;
    return $html;
  }

sub run
  {
    my %opts = @_;
    my $num_seqs = $opts{num_seqs} || $NUM_SEQS;
    my $spike_len = $opts{spike};
    my $iw = $opts{iw};
    my $ih = $opts{ih};
    my $feat_h = $opts{fh};
    my $show_gc = $opts{gc};
    my $show_nt = $opts{nt};
    my $color_hsp = $opts{colorhsp};
    my $hsp_label = $opts{hsplabel};
    my $draw_model = $opts{draw_model};
    my $hsp_overlap_limit = $opts{hsp_overlap_limit};
    my $overlap_adjustment = $opts{overlap};
    my $hiqual = $opts{hiqual};
    my $hsp_limit = $opts{hsplim};
    my $hsp_limit_num = $opts{hsplimnum};
    my $show_hsps_with_stop_codon = $opts{showallhsps};
    my $padding = $opts{padding};
    my ($analysis_program, $param_string, $parser_opts) = get_algorithm_options(%opts);
    my $basefilename = $opts{basefile};
    my $show_spike = $opts{show_spike} || 0;
    my $pad_gs = $opts{pad_gs} || 0;
    my $color_overlapped_features = $opts{color_overlapped_features};
    my $hsp_overlap_length = $opts{hsp_overlap_length};
    my $message;
    $cogeweb = initialize_basefile(basename=>$basefilename, prog=>"GEvo");
    my @hsp_colors;
    for (my $i = 1; $i <= num_colors($num_seqs); $i++)
      {
	my $r = $opts{"r$i"};
	my $g = $opts{"g$i"};
	my $b = $opts{"b$i"};
	my @tmp;
	foreach my $color ($r, $g, $b)
	  {
	    $color = 0 unless $color && $color =~ /^\d+$/;
	    $color = 0 if $color < 0;
	    $color = 255 if $color > 255;
	    push @tmp, $color;
	  }
	push @hsp_colors,\@tmp;
      }



    my $stagger_label = $hsp_label && $hsp_label =~ /staggered/i ? 1 : 0;
    my $feature_labels = !$hsp_label ? 0 : 1;
    my $form = $FORM;
    my $gevo_link = $form->url."?prog=$analysis_program";
    $gevo_link .= ";spike_len=$spike_len";
    
    my @sets;
    my $html;
    my $t1 = new Benchmark;
    my $spike_seq;
    my $seqcount = 1;
    for (my $i = 1; $i <= $num_seqs; $i++)
      {
	my $skip_seq =$opts{"skip_seq$i"};
	next if $skip_seq;
	my $accn = $opts{"draccn$i"};
	my $featid = $opts{"featid$i"};
	my $feat = $coge->resultset('Feature')->find($featid) if $featid;
	my $dsid;
	$dsid = $feat->dataset_id if $feat;
	$dsid = $opts{"dsid$i"} unless $dsid;
	$dsid = $opts{"posdsid$i"} unless $dsid;
	my $chr;
	$chr = $feat->chromosome if $feat;
	$chr = $opts{"chr$i"} unless $chr;

	my $drup = $opts{"drup$i"};
	my $drdown = $opts{"drdown$i"};
	$drup += $pad_gs;
	$drdown += $pad_gs;
	my $pos = $opts{"pos$i"};
	my $gbaccn = $opts{"gbaccn$i"};
	my $gbstart = $opts{"gbstart$i"};
	$gbstart = 1 unless defined $gbstart;
	my $gblength = $opts{"gblength$i"};

	my $dirseq = $opts{"dirseq$i"};
	my $dirstart = $opts{"dirstart$i"};
	my $dirlength = $opts{"dirlength$i"};

	my $rev = $opts{"rev$i"};
	my $mask_cds_flag = $opts{"maskcds$i"};
	my $mask_ncs_flag = $opts{"maskncs$i"};
	my ($up, $down, $seq);
	my ($file, $obj);
	my $reference_seq =$opts{"ref_seq$i"};
	my $repeat_mask =$opts{"repmask$i"};
	next unless $accn || $featid || $gbaccn || $dirseq|| $pos;
	$gevo_link .= ";accn$seqcount=".CGI::escape($accn) if $accn;
	$gevo_link .= ";x$seqcount=".CGI::escape($pos) if $pos;
	$gevo_link .= ";fid$seqcount=".CGI::escape($featid) if $featid;
	$gevo_link .= ";dsid$seqcount=".CGI::escape($dsid) if $dsid;
	$gevo_link .= ";chr$seqcount=".CGI::escape($chr) if $chr;
	$gevo_link .= ";dr$seqcount"."up=$drup" if defined $drup;
	$gevo_link .= ";dr$seqcount"."down=$drdown" if defined $drdown;
	$gevo_link .= ";gbaccn$seqcount=".CGI::escape($gbaccn) if $gbaccn;
	$gevo_link .= ";gbstart$seqcount=$gbstart" if $gbstart;
	$gevo_link .= ";gblength$seqcount=$gblength" if $gblength;
	$gevo_link .= ";rev$seqcount=1" if $rev;
	$gevo_link .= ";nref$seqcount=1" unless $reference_seq;
	$seqcount++;
#	print STDERR "pos: $pos, dsid: $dsid, chr: $chr\n";
	if ($featid || $pos)
	  {
	    $obj = get_obj_from_genome_db( accn=>$accn, featid=>$featid, pos=>$pos, dsid=>$dsid, rev=>$rev, up=>$drup, down=>$drdown, chr=>$chr );
	    if ($obj)
	      {
		($file, $spike_seq, $seq) = 
		  generate_seq_file(obj=>$obj,
				    mask_cds=>$mask_cds_flag,
				    mask_ncs=>$mask_ncs_flag,
				    spike_len=>$spike_len,
				    seq_num=>$i,
				    repeat_mask=>$repeat_mask,
				   );
		$up = $drup;
		$down = $drdown;
	      }
	    else
	      {
		$message .=  "Unable to generate sequence for sequence $i.  (Skipped.)";
	      }
	  }
 	elsif ($dirseq )
 	  {
 	    ($obj) = generate_obj_from_seq($dirseq, $i, $rev);
	    $dirlength = length($obj->sequence)-$dirstart+1 unless $dirlength;
 	    if ($obj)
	      {
		#add an anchor
		$obj->add_feature(
			      type=>"direct sequence submission",
			      location=>(1-$dirstart+1)."..".(length($obj->sequence)-$dirstart+1),
			      strand=>1,
			      qualifiers=>{
					   type=>"anchor",
					   names=>[$obj->accn],
					  }
			     );

		($file, $spike_seq, $seq) = 
		  generate_seq_file (
				     obj=>$obj,
				     mask_cds=>$mask_cds_flag,
				     mask_ncs=>$mask_ncs_flag,
				     startpos=>$dirstart,
				     length=>$dirlength,
				     spike_len=>$spike_len, 
				     seq_num=>$i,
				     repeat_mask=>$repeat_mask,
				    );
		$obj->start($dirstart);
		$obj->stop($dirstart+length($seq)-1);
		$obj->chromosome(1);
#		$obj->dataset("NA");
		$up = $dirstart;
		$down = $dirlength;
	      }
	    else
	      {
		$message .= "Problem with direct sequence submission\n";
	      }
	  }
	elsif ($gbaccn )
	  {
	    $obj = new CoGe::Accessory::GenBank;
	    $obj->add_gene_models(1); #may want to make a user selectable option
 	    my ($res, $error) = $obj->get_genbank_from_ncbi($gbaccn, $rev);
	    $message .= $error."\n" unless $res;
	    
	    if ($obj->accn)
	      {
		#add an anchor
		$obj->add_feature(
				  type=>"genbank entry",
				  location=>(1-$gbstart+1)."..".(length($obj->sequence)-$dirstart+1),
				  strand=>1,
				  qualifiers=>{
					       type=>"anchor",
					       names=>[$gbaccn],
					      }
				 );

		($file, $spike_seq, $seq) = 
		  generate_seq_file (
				     obj=>$obj,
				     mask_cds=>$mask_cds_flag,
				     mask_ncs=>$mask_ncs_flag,
				     startpos=>$gbstart,
				     length=>$gblength,
				     spike_len=>$spike_len,
				     seq_num=>$i,
				     repeat_mask=>$repeat_mask,
				    );
		$obj->start($gbstart);
		$obj->stop($gbstart + length($seq)-1);
		$obj->chromosome("?") unless $obj->chromosome;
		$up = $gbstart;
		$down = $gblength;
	      }
	    else
	      {
		$message .= "No GenBank entry found for $gbaccn\n";
	      }
 	  }
	next unless $obj;
	unless ($show_spike)
	  {
	    $seq =~ s/N*$spike_seq$//;
	    $seq =~ s/N*$spike_seq$//;
	  }
	$obj->sequence($seq);

	if ($obj && $obj->sequence)
	  {
	    #need to check for duplicate accession names -- sometimes happens and major pain in the ass for other parts of the code
	    my $accn = $obj->accn;
	    my $count = 0;
	    foreach my $accn2 (map {$_->{obj}->accn()} @sets)
	      {
		$accn2 =~ s/\*\*\d+\*\*$//;
		$count++ if $accn eq $accn2;
	      }
	    $accn .= "**".($count+1)."**" if $count;
	    $obj->accn($accn);
	    push @sets, {
			 obj=>$obj,
			 file=>$file,
#			 file_begin=>$file_begin,
#			 file_end=>$file_end,
			 accn=>$accn,
			 rev=>$rev,
			 up=>$up,
			 down=>$down,
			 spike_seq=>$spike_seq,
			 reference_seq=>$reference_seq,
			 seq_num=>$i,
			};
	  }
	else
	  {
#	    push @sets, {seq_num=>$i};
	  }
      }
    $seqcount--;

    $gevo_link .= ";num_seqs=".$seqcount;
    $gevo_link .= ";hsp_overlap_limit=".$hsp_overlap_limit if defined $hsp_overlap_limit;
    unless (@sets >1)
      {
	$message .= "Problem retrieving information.  Please check submissions.\n";
	return '', '', '', '',0,$message;
      }
    my $t2 = new Benchmark;
    # set up output page
    
    # run bl2seq
    
    my $analysis_reports;

    if ($analysis_program eq "blastz")
      {
	$analysis_reports = run_blastz(sets=>\@sets, params=>$param_string, parser_opts=>$parser_opts);
      }
    elsif ($analysis_program eq "LAGAN")
      {
	$analysis_reports = run_lagan (sets=>\@sets, params=>$param_string, parser_opts=>$parser_opts);
      }
    elsif ($analysis_program eq "CHAOS")
      {
	$analysis_reports = run_chaos (sets=>\@sets, params=>$param_string, parser_opts=>$parser_opts);
      }
    elsif ($analysis_program eq "DiAlign_2")
      {
	($analysis_reports, $analysis_program,$message) = run_dialign (sets=>\@sets, params=>$param_string, parser_opts=>$parser_opts);
      }
    else
      {
	$analysis_reports = run_bl2seq(sets=>\@sets, params=>$param_string, parser_opts=>$parser_opts, blast_program=>$analysis_program, spike_seq=>$spike_seq);
      }
    $analysis_reports = [] unless ref($analysis_reports) =~ /ARRAY/i;
    write_log($message, $cogeweb->logfile) if $message;
   #sets => array or data for blast
   #blast_reports => array of arrays (report, accn1, accn2, parsed data)
    my $t3 = new Benchmark;
    my $count = 1;
    my $frame_height;
    initialize_sqlite();
    my @gfxs;
    foreach my $item (@sets)
      {
	my $obj = $item->{obj};
	next unless $obj->sequence;
	my $rev = $item->{rev};
	my $spike_seq = $item->{spike_seq};
	if ($obj)
	  {
	    write_log("generating image ($count/".scalar @sets.")for ".$obj->accn, $cogeweb->logfile);
	    $count++;
#	    my ($image, $map, $mapname, $gfx, $eval_cutoff, $feat_count, $overlap_count) = 
	    my ($gfx, $stats) = 
	      generate_image(
			     set=>$item,
			     start=>1,
			     stop=>length($obj->sequence),
			     data=>$analysis_reports,
			     iw=>$iw,
			     ih=>$ih,
			     fh=>$feat_h,
			     show_gc=>$show_gc,
			     show_nt=>$show_nt,
			     stagger_label=>$stagger_label,
			     overlap_adjustment=>$overlap_adjustment,
			     feature_labels=>$feature_labels,
			     hsp_limit=>$hsp_limit,
			     hsp_limit_num=>$hsp_limit_num,
			     color_hsp=>$color_hsp,
			     hsp_colors=>\@hsp_colors,
			     spike_sequence=>$spike_seq,
			     show_hsps_with_stop_codon => $show_hsps_with_stop_codon,
			     hiqual=>$hiqual,
			     padding=>$padding,
#			     seq_num=>$item->{seq_num},
			     reverse_image=>$rev,
			     draw_model=>$draw_model,
			     hsp_overlap_limit=>$hsp_overlap_limit,
			     color_overlapped_features=>$color_overlapped_features,
			     hsp_overlap_length=>$hsp_overlap_length,
			    );
#	    $item->{feat_count} = $feat_count;
#	    $item->{overlap_count} = $overlap_count;
	    my $filename = $cogeweb->basefile."_".$item->{seq_num}.".png";
	    $filename = check_filename_taint($filename);
	    my $image = basename($filename);
	    $item->{image}=$image;
	    $item->{stats}=$stats;
	    $gfx->generate_png(file=>$filename);
	    generate_image_db(set=>$item, gfx=>$gfx);
	    $frame_height += $gfx->ih + $gfx->ih*.1;
	  }
      }
    my $t3_5 = new Benchmark;
    image_db_create_hsp_pairs();

    my $t4 = new Benchmark;
    $html .= qq{<DIV id=flash_viewer></DIV>};
    $html .= qq{<table>};
    $html .= qq{<tr valign=top><td class = small>Alignment reports};
    my $stats_file = $cogeweb->basefile."_stats.txt";
    $stats_file = check_filename_taint($stats_file);

    if ($analysis_reports && @$analysis_reports)
      {
	foreach my $item (@$analysis_reports)
	  {
	    my $report = $item->[0];
	    my $accn1 = $item->[1];
	    my $accn2 = $item->[2];
	    my $basereportname = basename( $report );
	    $basereportname = $TEMPURL . "/$basereportname\n";
	    $html .= "<div><font class=small><A HREF=\"$basereportname\" target=_new>View alignment output for $accn1 versus $accn2</A></font></DIV>\n";
	  }
      }
    else
      {
	$html .= "<div class=small>No alignment reports were generated</DIV>\n";
      }
    $html .= qq{<td class = small>Fasta files};
    foreach my $item (@sets)
      {
	next unless $item->{file};
	my $basename = $TEMPURL."/".basename ($item->{file});
	print STDERR "basename is undefined: $basename\n" if $basename =~ /defined/i;
	my $accn = $item->{accn};
	$html .= "<div><font class=small><A HREF=\"$basename\" target=_new>Fasta file for $accn</A></font></DIV>\n";
      }
    $html .= qq{<td class = small><a href = "http://baboon.math.berkeley.edu/mavid/gaf.html">GAF</a> annotation files};
    foreach my $item (@sets)
      {
	my $anno_file = generate_annotation(%$item);
	next unless $anno_file;
	my $basename = $TEMPURL."/".basename ($anno_file);
	my $accn = $item->{accn};
	$html .= "<div><font class=small><A HREF=\"$basename\" target=_new>Annotation file for $accn</A></font></DIV>\n";
     }
    $html .= qq{<td class = small>SQLite db};
    my $dbname = $TEMPURL."/".basename($cogeweb->sqlitefile);
    
    $html .= "<div class=small><A HREF=\"$dbname\" target=_new>SQLite DB file</A></DIV>\n";
    $html .= qq{<td class = small>Log File};
    my $logfile = $TEMPURL."/".basename($cogeweb->logfile);
    $html .= "<div class=small><A HREF=\"$logfile\" target=_new>Log</A></DIV>\n";
    my $tiny = get("http://tinyurl.com/create.php?url=$gevo_link");
    ($tiny) = $tiny =~ /<b>(http:\/\/tinyurl.com\/\w+)<\/b>/;
    $html .= qq{<td class = small>GEvo Link<div class=small><a href=$tiny target=_new>$tiny<br>(See log file for full link)</a></div>};
    $html .=qq{<a href="GEvo_direct.pl?name=$basefilename" target=_new>Results only</a>};
    $html .= qq{<td class = small>Gene Homology Stats:<table>};


    write_log("Summary of overlap data for genes only:", $stats_file);
    write_log(join ("\t", qw(REGION1 REGION2 LENGTH OVERLAP/TOTAL PERCENT_OVERLAP)),$stats_file);
    foreach my $item (@sets)
      {
	$html .= "<tr class=small>";
	my @data = ("<td>".$item->{obj}->accn,commify(length $item->{obj}->sequence)." bp",$item->{stats}{feat_counts}{gene}{overlap}."/".$item->{stats}{feat_counts}{gene}{count});
	push @data, sprintf("%.2f", $item->{stats}{feat_counts}{gene}{overlap}/$item->{stats}{feat_counts}{gene}{count}*100)."%" if $item->{stats}{feat_counts}{gene}{count};
	$html .= join ("<td nowrap>",@data);
	shift @data;
	write_log(join ("\t",$item->{obj}->accn,"all",@data), $stats_file);
      }

    write_log("\n",$stats_file);
    write_log("Features overlapped by HSPs:", $stats_file);
    write_log(join ("\t", qw(REGION1 REGION2 TYPE OVERLAP_COUNT TOTAL PERCENT_W_OVERLAP)),$stats_file);
    foreach my $item (@sets)
      {
	my $stats = $item->{stats};
	foreach my $type (sort keys %{$stats->{feat_counts}})
	  {
	    foreach my $accn1 (sort keys %{$stats->{data}})
	      {
		foreach my $accn2 (sort keys %{$stats->{data}{$accn1}})
		  {
		    next if $type eq "hsp";
		    my $overlap_count = $stats->{data}{$accn1}{$accn2}{$type}{overlapped_hsps} || 0;
		    my $percent = sprintf("%.2f",$overlap_count/$stats->{feat_counts}{$type}{count}*100) if $stats->{feat_counts}{$type}{count};
		    write_log(join ("\t", $accn1, $accn2, $type, $overlap_count, $stats->{feat_counts}{$type}{count}, $percent), $stats_file);
		  }
	      }
	  }
      }


    write_log("\n",$stats_file);
    write_log("Average percent identity for all HSPs (normalized to HSP length):", $stats_file);
    write_log(join ("\t", qw(REGION1 REGION2 AVG_PERCENT_ID)),$stats_file);
    foreach my $item (@sets)
      {
	my $stats = $item->{stats};
	foreach my $accn1 (sort keys %{$stats->{data}})
	  {
	    foreach my $accn2 (sort keys %{$stats->{data}{$accn1}})
	      {
		my $len= 0;
		map {$len+=$_->[1]} @{$stats->{data}{$accn1}{$accn2}{hsp}};
		my $pid = 0;
		map {$pid+=$_->[1]*$_->[0]} @{$stats->{data}{$accn1}{$accn2}{hsp}};
		write_log("$accn1\t$accn2\t".sprintf("%.2f",$pid/$len), $stats_file);
	      }
	  }
      }
    my $stats_url = $TEMPURL."/".basename($stats_file);
    $html .= qq{<tr class=small><td><a href=$stats_url target=_new>Stats file</a>};
    $html .= qq{</table>};
    $html .= qq{</table>};

    my $template = HTML::Template->new(filename=>'/opt/apache/CoGe/tmpl/box.tmpl');
    my $results_name = "Results: $analysis_program";
    $results_name .= qq{ <span class="small">(spike sequence filter length: $spike_len)</span>} if $spike_len;
    $template->param(BOX_NAME=>$results_name);
    $template->param(BODY=>$html);
    my $outhtml = $template->output;
    my $t5 = new Benchmark;
    my $db_time = timestr(timediff($t2,$t1));
    my $blast_time = timestr(timediff($t3,$t2));
    my $image_time = timestr(timediff($t3_5,$t3));
    my $db_hsp_time = timestr(timediff($t4,$t3_5));
    my $html_time = timestr(timediff($t5,$t4));
    my $bench =  qq{
GEvo Benchmark: $DATE
Time to get sequence                              : $db_time
Time to run $analysis_program                     : $blast_time
Time to generate images, maps, and sqlite database: $image_time
Time to find and update sqlite database for HSPs  : $db_hsp_time
Time to process html                              : $html_time
};
    print STDERR $bench if $BENCHMARK;
    write_log($bench, $cogeweb->logfile);
    write_log("Finished!", $cogeweb->logfile);
    write_log("\nGEvo link: $gevo_link\n", $cogeweb->logfile);
    write_log("Tiny url: $tiny", $cogeweb->logfile);
    $count--;
    return $outhtml, $iw+400, $frame_height, $cogeweb->basefilename,$count,$message;

}


sub generate_image
  {
    my %opts = @_;
    my $set = $opts{set};
    my $gbobj = $set->{obj};
    my $start = $opts{start};
    my $stop = $opts{stop};
    my $data = $opts{data};
    my $masked_exons = $opts{mask};
    my $masked_ncs = $opts{mask_ncs};
    my $spike_seq = $opts{spike_sequence};
    my $iw = $opts{iw} || 1600;
    my $ih = $opts{ih} || 200;
    my $fh = $opts{fh} || 25;
    my $show_gc = $opts{show_gc};
    my $show_nt = $opts{show_nt};
    my $reverse_image = $opts{reverse_image};
    my $stagger_label = $opts{stagger_label};
    my $overlap_adjustment = $opts{overlap_adjustment};
    my $feature_labels = $opts{feature_labels};
    my $hsp_limit = $opts{hsp_limit};
    my $color_hsp = $opts{color_hsp};
    my $hsp_limit_num = $opts{hsp_limit_num};
    my $eval_cutoff = $opts{eval_cutoff};
    my $hsp_colors = $opts{hsp_colors};
    my $show_hsps_with_stop_codon = $opts{show_hsps_with_stop_codon};
    my $hiqual = $opts{hiqual};
    my $padding = $opts{padding} || 5;
#    my $seq_num = $opts{seq_num};
    my $draw_model = $opts{draw_model};
    my $hsp_overlap_limit = $opts{hsp_overlap_limit};
    my $color_overlapped_features = $opts{color_overlapped_features};
    my $hsp_overlap_length = $opts{hsp_overlap_length};
    my $graphic = new CoGe::Graphics;
    my $gfx = new CoGe::Graphics::Chromosome;
#    print STDERR "$start"."::"."$stop"." -- ".length($gbobj->sequence)."\n";
    $graphic->initialize_c (
			    c=>$gfx,
			    iw=>$iw,
			    start=> $start,
			    stop => $stop,
			    draw_chr=>1,
			    draw_ruler=>1,
			    draw_chr_end=>0,
			    chr_start_height=>$ih,
			    chr_mag_height=>5,
			    feature_start_height=>$fh,
			    mag=>0,
			    mag_off=>1,
			    chr_length => length($gbobj->sequence),
			    fill_labels=>1,
			    forcefit=>1,
			    minor_tick_labels=>1,
			    feature_labels=>$feature_labels,
			    draw_hi_qual=>$hiqual,
			    padding=>$padding,
			   );
 #   print STDERR $gfx->_region_start ."--".$gfx->_region_stop."\n";
    $gfx->overlap_adjustment(1);
    $gfx->top_padding(15);
    $gfx->skip_duplicate_features(1);
    $gfx->DEBUG(0);
    $gfx->major_tick_labels(0);
    $graphic->process_nucleotides(c=>$gfx, seq=>$gbobj->sequence, layers=>{gc=>$show_gc, nt=>$show_nt});
    my $stats = process_hsps(
			  c=>$gfx, 
			  data=>$data, 
			  accn=>$gbobj->accn, 
			  rev=>$reverse_image, 
			  seq_length=> length($gbobj->sequence), 
			  stagger_label=>$stagger_label, 
			  hsp_limit=>$hsp_limit,
			  hsp_limit_num=>$hsp_limit_num,
			  gbobj=>$gbobj,
			  spike_seq=>$spike_seq,
			  color_hsp=>$color_hsp,
			  colors=>$hsp_colors,
			  show_hsps_with_stop_codon=>$show_hsps_with_stop_codon,
			  hsp_overlap_limit=>$hsp_overlap_limit,
			  hsp_overlap_length=>$hsp_overlap_length,
			 );
    my ($feat_counts) = process_features(c=>$gfx, obj=>$gbobj, start=>$start, stop=>$stop, overlap_adjustment=>$overlap_adjustment, draw_model=>$draw_model, color_overlapped_features=>$color_overlapped_features);
    $stats->{feat_counts} = $feat_counts;
    return ($gfx, $stats);
  }

sub initialize_sqlite
  {
    my %opts = @_;
    my $dbfile = $cogeweb->sqlitefile;
    return if -r $dbfile;
    my $dbh = DBI->connect("dbi:SQLite:dbname=$dbfile","","");
    my $create = qq{
CREATE TABLE image_data
(
id INTEGER PRIMARY KEY,
name varchar,
type varchar,
xmin integer,
xmax integer,
ymin integer,
ymax integer,
bpmin integer,
bpmax integer,
image_track integer,
image_id integer,
pair_id integer,
link varchar,
annotation text,
color varchar
)
};
    $dbh->do($create);
    $dbh->do('CREATE INDEX type ON image_data (type)');
    $dbh->do('CREATE INDEX image_id ON image_data (image_id)');
    $dbh->do('CREATE INDEX x_idx ON image_data (xmin,xmax)');
    $dbh->do('CREATE INDEX y_idx ON image_data (ymin,ymax)');
    $dbh->do('CREATE INDEX image_track ON image_data (image_track)');
    $dbh->do('CREATE INDEX pair_id ON image_data (pair_id)');

#id INTEGER PRIMARY KEY AUTOINCREMENT,
    $create = qq{
CREATE TABLE image_info
(
id INTEGER,
iname varchar,
title varchar,
px_width integer,
bpmin integer,
bpmax integer,
dsid integer,
chromosome varchar
)
};
#TODO: make sure to populate the bpmin and pbmax and image width!!!!!!!!!!!!!!!!!
    $dbh->do($create);
    $dbh->do('CREATE INDEX id ON image_info (id)');
    system "chmod +rw $dbfile";
  }

sub generate_image_db
  {
    my %args = @_;
    my $gfx = $args{gfx};
    my $set = $args{set};
    next unless $gfx;
    my $dbh = DBI->connect("dbi:SQLite:dbname=".$cogeweb->sqlitefile,"","");
    my $image = $set->{image};
    my $accn = $set->{accn};
    my $title;
    $title = $accn;
    $title .= " ".$set->{obj}->organism() if $set->{obj}->organism();
    $title .= "(chr: ".$set->{obj}->chromosome." ". $set->{obj}->start."-".$set->{obj}->stop.")" if defined $set->{up};
    $title .= qq! Reverse Complement! if $set->{rev};
    my $width = $gfx->image_width;
    my $dsid = $set->{obj}->dataset;
    my ($chr) =  $set->{obj}->chromosome =~ /(\d+)/;
    $chr = "NULL" unless $chr;
    $dsid = "NULL" unless $dsid;
    my $image_start = $set->{obj}->start;
    my $image_stop = $set->{obj}->stop;
    my $image_id = $set->{seq_num};
    my $statement = qq{
INSERT INTO image_info (id, iname, title, px_width,dsid, chromosome, bpmin, bpmax) values ($image_id, "$image", "$title", $width, "$dsid", $chr, $image_start, $image_stop)
};
    print STDERR $statement unless $dbh->do($statement);
    foreach my $feat ($gfx->get_feats)
      {
	if ($feat->fill)
	  {
	    next unless $feat->type eq "anchor";
	  }
	next unless $feat->image_coordinates;
	my $type = $feat->type;
	my $pair_id = "-99";
	my $coords = $feat->image_coordinates;
	$coords =~ s/\s//g;
	next if $feat->type eq "unknown";
	my $name = $feat->type =~ /HSP/i ? $feat->alt : $feat->label;
	$name .= "_".$feat->type;# unless $name;

	my $color = "NULL";
	if ($feat->type =~ /HSP/)
	  {
	    $color ="#";
	    foreach my $c (@{$feat->color})
	      {
		$color .= sprintf("%X",$c);
	      }
	  }
	#generate link
	my $link = $feat->link;
	$link = " " unless $link;
	$link =~ s/'//g if $link;
	#generate image track
	my $image_track = $feat->track;
	$image_track = "-".$image_track if $feat->strand =~ /-/;
	
	my ($xmin, $ymin, $xmax, $ymax) = split /,/, $coords;
	$xmin++;
	$xmax++;
	my $anno = $feat->description;
	#	    $anno =~ s/'|"//g;
	$anno =~ s/'//g if $anno;
	$anno =~ s/<br\/?>/&#10;/ig if $anno;
	$anno =~ s/\n/&#10;/g if $anno;
	$anno =~ s/[\[\]\(\)]/ /g if $anno;
	$anno = " " unless $anno;
	my $start = $feat->start;
	my $stop = $feat->stop;
	my $length_nt = $stop-$start+1;
	my $length_pix = $xmax-$xmin;
	next if !$feat->{anchor} && ($length_nt == 0 );
	$type = "anchor" if $feat->{anchor};
	my $bpmin = $set->{obj}->start+$start-1;
	my $bpmax = $set->{obj}->start+$stop-1;
	$statement = qq{
INSERT INTO image_data (name, type, xmin, xmax, ymin, ymax, bpmin,bpmax,image_id, image_track,pair_id, link, annotation, color) values ("$name", "$type", $xmin, $xmax, $ymin, $ymax, $bpmin, $bpmax,$image_id, "$image_track",$pair_id, '$link', '$anno', '$color')
};
	print STDERR $statement unless $dbh->do($statement);
      }
  }

sub image_db_create_hsp_pairs
  {
    my $dbh = DBI->connect("dbi:SQLite:dbname=".$cogeweb->sqlitefile,"","");
    my $query = qq{select id,name from image_data where type = "HSP"};
    my %ids;
    foreach my $item (@{$dbh->selectall_arrayref($query)})
      {
	push @{$ids{$item->[1]}},$item->[0];
      }
    foreach my $name (keys %ids)
      {
	if (@{$ids{$name}} != 2)
	  {
	    print STDERR "Problem with $name.  Retrieved", Dumper $ids{$name} if $DEBUG;
	    next;
	  }
	my ($id1, $id2) = @{$ids{$name}};
	my $statement = "update image_data set pair_id = $id1 where id = $id2";
	$dbh->do($statement);
	$statement = "update image_data set pair_id = $id2 where id = $id1";
	$dbh->do($statement);
      }
  }


sub process_features
  {
    #process features
    my %opts = @_;
    my $c = $opts{c};
    my $obj=$opts{obj};
    my $start=$opts{start};
    my $stop = $opts{stop};
    my $overlap = $opts{overlap_adjustment};
    my $draw_model = $opts{draw_model};
    my $color_overlapped_features = $opts{color_overlapped_features};
#    return unless $draw_model;
    my $accn = $obj->accn;
    my $track = 1;
    my %feat_counts;
    my @opts = ($start, $stop) if $start && $stop;
    unless (ref $obj)
      {
	warn "Possible problem with the object in process_features.  Returning";
	return 0;
      }
    my %prior_feat;
    foreach my $feat (sort {$a->start <=> $b->start} $obj->get_features())
      {
        my $f;
	my $type = $feat->type;
	my ($name) = sort { length ($b) <=> length ($a) || $a cmp $b} @{$feat->qualifiers->{names}} if ref ($feat->qualifiers) =~ /hash/i;
	my $anchor = $feat if ref ($feat->qualifiers) =~ /hash/i && $feat->qualifiers->{type} && $feat->qualifiers->{type} eq "anchor";
        if ($type =~ /pseudogene/i)
          {
	    next unless $draw_model eq "full";
	    $f = CoGe::Graphics::Feature::Gene->new();
	    $f->color([255,33,0,50]);
	    $f->order($track);
	    $f->overlay(1);
	    $f->mag(0.5);
          }
        elsif ($type =~ /Gene/i)
          {
	    next unless $draw_model eq "full" || $draw_model eq "gene" || $anchor;
	    $f = CoGe::Graphics::Feature::Gene->new();
	    #$f->color([255,0,0,50]);
        $f->color([219,219,219,50]);
	    $f->order($track);
	    $f->overlay(1);
	    $f->mag(0.5);
          }
        elsif ($type =~ /CDS/i)
          {
	    next unless $draw_model eq "full" || $draw_model eq "CDS";
	    $f = CoGe::Graphics::Feature::Gene->new();
	    $f->color([0,255,0, 50]);
	    $f->order($track);
	    $f->overlay(3);
	    if ($accn)
	      {
		foreach my $name (@{$feat->qualifiers->{names}})
		  {
            my $cleaned_name = $name;
            $cleaned_name =~ s/[\(\)]//g;
		    if ($accn =~ /^$cleaned_name\(?\d*\)?$/i)
		      {
			$f->color([255,255,0]) ;
			$f->label($name);
		      }
		  }
	      }
	    
		
          }
        elsif ($type =~ /mrna/i)
          {
	    next unless $draw_model eq "full" || $draw_model eq "mRNA";
        	$f = CoGe::Graphics::Feature::Gene->new();
        	$f->color([0,0,255, 50]);
        	$f->order($track);
		$f->overlay(2);
		$f->mag(0.75);
          }
        elsif ($type =~ /rna/i)
          {
	    next unless $draw_model eq "full" || $draw_model eq "RNA";
        	$f = CoGe::Graphics::Feature::Gene->new();
        	$f->color([200,200,200, 50]);
        	$f->order($track);
		$f->overlay(2);
		if ($accn)
		  {
		    foreach my $name (@{$feat->qualifiers->{names}})
		      {
			if ($accn =~ /^$name\(?\d*\)?$/i)
			  {
			    $f->color([255,255,0]) ;
			    $f->label($name);
			  }
		      }
		  }

          }
	elsif ($type =~ /anchor/)
	  {
	    $f = CoGe::Graphics::Feature::NucTide->new({type=>'anchor',start=>$feat->blocks->[0][0], strand=>1});
	    $f->color([0,0,255]);
	    $f->type($type);
	    $f->description($feat->annotation);
	    $c->add_feature($f);
	    $f = CoGe::Graphics::Feature::NucTide->new({type=>'anchor',start=>$feat->blocks->[0][0], strand=>-1});
	    $f->color([0,0,255]);
	    $f->type($type);
	    $f->description($feat->annotation);
	    $c->add_feature($f);
	    next;
	  }
	elsif ($type =~ /cns/i)
	  {
	    $f = CoGe::Graphics::Feature::HSP->new({start=>$feat->blocks->[0][0], stop=>$feat->blocks->[0][1]});
	    $f->color([255,51,0]);
	    $f->order($track);
	    $f->overlay(4);
	    $f->type($type);
	    $f->force_draw(1);
	    $c->add_feature($f);
	    next;
	  }
	#need to create an anchor if this feature is an anchor, but not to be drawn
	if ($anchor && !$f)
	  {
	    $f = CoGe::Graphics::Feature->new({type=>'anchor',start=>$feat->blocks->[0][0], stop=>$feat->blocks->[0][1], strand=>1});
	    $f->{anchor}=1;
	    $f->order(0);
	    $f->skip_overlap_search(1);
	    $f->description("auto generated anchor for GEvo");
	    $c->add_feature($f);
	    next;
	  }
        next unless $f;
	$f->color([255,0,255]) if $color_overlapped_features && $feat->qualifiers->{overlapped_hsp};
	my $strand = 1;
 	$strand = -1 if $feat->location =~ /complement/;
	foreach my $block (@{$feat->blocks})
	  {
	    $block->[0] =1 unless $block->[0]; #in case $block is set to 0
	    $f->add_segment(start=>$block->[0], stop=>$block->[1]);
	    $f->strand($strand);
	    print STDERR "\t", join ("-", @$block),"\n" if $DEBUG;
	  }

	print STDERR $name,"\n\n" if $DEBUG;
        $f->type($type);
	$f->description($feat->annotation);
	$f->link("FeatList.pl?fid=".$feat->qualifiers->{id}) if $feat->qualifiers->{id};
	my $foverlap = $overlap ? 0 : 1; #I think I need this to get this to work as expected
	$f->skip_overlap_search($foverlap);
	$f->{anchor}=1 if $anchor;
        $c->add_feature($f);
	print STDERR $feat->type, " ", $feat->start,"-", $feat->stop," ",$strand," :: ",$start,"-",$stop,"\n";

	if ($prior_feat{$feat->type}{$strand})
	  {
	    if ($feat->start < $prior_feat{$feat->type}{$strand}->stop)
	      {
		print STDERR "IN skip\n";
		$prior_feat{$feat->type}{$strand} = $feat;
		next;
	      }
	  }
	next if $feat->start > $stop;
	next if $feat->stop < $start;
	print STDERR "!!\n";
	$feat_counts{$feat->type}{count}++;
	$feat_counts{$feat->type}{overlap}++ if $feat->qualifiers->{overlapped_hsp};
	$prior_feat{$feat->type}{$strand} = $feat;
      }
    return (\%feat_counts);
  }

sub process_hsps
  {
    my %opts = @_;
    my $c = $opts{c};
    my $data = $opts{data};
    my $accn = $opts{accn};
    my $reverse = $opts{rev};
    my $seq_len = $opts{seq_length};
    my $stagger_label = $opts{stagger_label};<commit_msg>removed debuggin messages<commit_after>#!/usr/bin/perl -w
use strict;
use warnings;
use CGI;
use CGI::Carp 'fatalsToBrowser';
use CGI::Cookie;
use CGI::Ajax;
use HTML::Template;
use Data::Dumper;
use File::Basename;
use File::Temp;
use CoGe::Accessory::GenBank;
use CoGe::Accessory::LogUser;
use CoGe::Accessory::Web;
use CoGe::Accessory::Restricted_orgs;
use CoGe::Accessory::bl2seq_report;
use CoGe::Accessory::blastz_report;
use CoGe::Accessory::lagan_report;
use CoGe::Accessory::chaos_report;
use CoGe::Accessory::dialign_report;
use CoGe::Accessory::dialign_report::anchors;
use CoGe::Graphics;
use CoGe::Graphics::Chromosome;
use CoGe::Graphics::Feature;
use CoGe::Graphics::Feature::Gene;
use CoGe::Graphics::Feature::NucTide;
use CoGe::Graphics::Feature::GAGA;
use CoGe::Graphics::Feature::Exon_motifs;
use CoGe::Graphics::Feature::AminoAcid;
use CoGe::Graphics::Feature::Domain;
use CoGe::Graphics::Feature::HSP;
use CoGeX;
use CoGeX::Feature;
use DBIxProfiler;
use DBI;
use LWP::Simple;
use Parallel::ForkManager;
use Statistics::Basic;
use Benchmark qw(:all);

# for security purposes

$ENV{PATH} = "/opt/apache/CoGe/";
delete @ENV{ 'IFS', 'CDPATH', 'ENV', 'BASH_ENV' };
#for chaos
$ENV{'LAGAN_DIR'} = '/opt/apache/CoGe/bin/lagan/';
#for dialign
$ENV{'DIALIGN2_DIR'} = '/opt/apache/CoGe/bin/dialign2_dir/';
use vars qw( $PAGE_NAME $DATE $DEBUG $BL2SEQ $BLASTZ $LAGAN $CHAOS $DIALIGN $TEMPDIR $TEMPURL $USER $FORM $cogeweb $BENCHMARK $coge $NUM_SEQS $MAX_SEQS $REPEATMASKER $MAX_PROC);
$PAGE_NAME = "GEvo.pl";
$BL2SEQ = "/usr/bin/bl2seq ";
$BLASTZ = "/usr/bin/blastz ";
$LAGAN = "/opt/apache/CoGe/bin/lagan/lagan.pl";
$CHAOS = "/opt/apache/CoGe/bin/lagan/chaos_coge";
$DIALIGN = "/opt/apache/CoGe/bin/dialign2_dir/dialign2-2_coge";
$REPEATMASKER = "/opt/apache/CoGe/bin/RepeatMasker/RepeatMasker";
$TEMPDIR = "/opt/apache/CoGe/tmp/GEvo";
$TEMPURL = "/CoGe/tmp/GEvo";
$MAX_PROC=8;
# set this to 1 to print verbose messages to logs
$DEBUG = 0;
$BENCHMARK = 0;
$NUM_SEQS = 3;
$MAX_SEQS = 21;
$| = 1; # turn off buffering 
$DATE = sprintf( "%04d-%02d-%02d %02d:%02d:%02d",
		 sub { ($_[5]+1900, $_[4]+1, $_[3]),$_[2],$_[1],$_[0] }->(localtime));

$FORM = new CGI;
$CGI::POST_MAX= 60 * 1024 * 1024; # 24MB
$CGI::DISABLE_UPLOADS = 0; 
($USER) = CoGe::Accessory::LogUser->get_user();

$coge = CoGeX->dbconnect();
#$coge->storage->debugobj(new DBIxProfiler());
#$coge->storage->debug(1);

my %ajax = CoGe::Accessory::Web::ajax_func();
$ajax{dataset_search} = \&dataset_search; #override this method from Accessory::Web for restricted organisms
$ajax{feat_search} = \&feat_search; 
my $pj = new CGI::Ajax(
		       run=>\&run,
		       loading=>\&loading,
		       add_seq=>\&add_seq,
		       get_file=>\&get_file,
		       gen_go_run=>\&gen_go_run,
		       gen_hsp_colors =>\&gen_hsp_colors,
		       save_settings_gevo=>\&save_settings_gevo,
		       reset_settings_gevo=>\&reset_settings_gevo,
		       %ajax,
		      );
$pj->JSDEBUG(0);
$pj->DEBUG(0);
#$pj->js_encode_function('escape');
print $pj->build_html($FORM, \&gen_html);

#$USER="elyons";print gen_html();

sub loading
  {
    my $message = shift || "Generating results. . .";
    return qq{<div class="dna"><div id="loading">$message</div></div>}; 
  }

sub gen_html
  {
    my $html;# =  "Content-Type: text/html\n\n";
    unless ($USER)
      {
	$html = login();
      }
    else
      {
	my $template = HTML::Template->new(filename=>'/opt/apache/CoGe/tmpl/generic_page.tmpl');
	$template->param(LOGO_PNG=>"GEvo-logo.png");
	$template->param(TITLE=>'Genome Evolution Analysis');
	$template->param(HELP=>'GEvo');
	$template->param(USER=>$USER);
	$template->param(DATE=>$DATE);
	$template->param(NO_BOX=>1);
	$template->param(BODY=>gen_body());
	my $prebox = HTML::Template->new(filename=>'/opt/apache/CoGe/tmpl/GEvo.tmpl');
	$prebox->param(RESULTS_DIV=>1);
	$template->param(PREBOX=>$prebox->output);
	
	$html .= $template->output;
      }
    return $html;
  }

sub gen_body
  {
    my $form = $FORM;
    my $prefs = load_settings(user=>$USER, page=>$PAGE_NAME);
    my $num_seqs = get_opt(params=>$prefs, form=>$form, param=>'num_seqs');
    $num_seqs = $NUM_SEQS unless defined $num_seqs;
    $MAX_SEQS = 10 if $form->param('override');
    my $message;
    if (! ($num_seqs =~ /^\d+$/) )
      {
	$message .= "Problem with requested number of sequences: '$num_seqs'.  Defaulting to $NUM_SEQS input sequences.";
	$num_seqs = $NUM_SEQS;
      }
    elsif ($num_seqs < 2)
      {
	$message .= "Minimum number of sequences to compare is two.";
	$num_seqs = 2;
      }
    elsif ($num_seqs > $MAX_SEQS)
      {
	$message .= "Maximum number of sequence is $MAX_SEQS.";
	$num_seqs = $MAX_SEQS;
      }
    my @seq_nums;
    my @seq_sub;
    my $autosearch_string;
    for (my $i = 1; $i <= $num_seqs; $i++)
      {
	my $draccn;
	$draccn= $form->param("accn".$i) if $form->param("accn".$i);
	my $pos = $form->param("x".$i) if $form->param("x".$i);
	my $chr = $form->param("chr".$i) if $form->param("chr".$i);
	my $drfid = $form->param("fid".$i) if $form->param("fid".$i);
	if ($drfid && !$draccn)
	  {
	    my $feat = $coge->resultset('Feature')->find($drfid);
	    ($draccn) = $feat->primary_name if $feat;
	  }
	my $drup = $form->param('dr'.$i.'up') if defined $form->param('dr'.$i.'up');
	my $drdown = $form->param('dr'.$i.'down') if defined $form->param('dr'.$i.'down');
	$drup = $form->param('drup'.$i) if defined $form->param('drup'.$i);
	$drdown = $form->param('drdown'.$i) if defined $form->param('drdown'.$i);
	$drup = 10000 unless defined $drup;
	$drdown = 10000 unless defined $drdown;
#	$drup += $pad_gs if $pad_gs;
#	$drdown += $pad_gs if $pad_gs;
	my $dsid = $form->param('dsid'.$i) if $form->param('dsid'.$i);
	my $gbaccn = $form->param("gbaccn".$i) if $form->param("gbaccn".$i);
	my $gbstart = $form->param("gbstart".$i) if $form->param("gbstart".$i);
	$gbstart = 1 unless defined $gbstart;
	my $gblength = $form->param("gblength".$i) if $form->param("gblength".$i);
	my $revy = "checked" if $form->param('rev'.$i);
	my $revn = "checked" unless $revy;
	my $maskexonon = "checked" if $form->param('maskexon'.$i);
	my $maskexonoff = "checked" unless $maskexonon;
	my $refn = "checked" if $form->param('nref'.$i);
	my $refy = "checked" unless $refn;
	$autosearch_string .= 'if ($'.qq!('#accn$i').val()) {dataset_search(['args__accn','accn$i','args__num', 'args__$i'!;
	$autosearch_string .= ", 'args__dsid', 'args__$dsid'" if $dsid;
	$autosearch_string .= ", 'args__featid', 'args__$drfid'" if $drfid;
	$autosearch_string .=qq!],[feat_search_chain]);}!;
	my $dsinfo = $dsid ? qq{<input type="hidden" id="posdsid$i" value="$dsid">} : qq{<input type="hidden" id="posdsid$i">};
	$dsinfo .= $chr ? qq{<input type="hidden" id="chr$i" value="$chr">} : qq{<input type="hidden" id="chr$i">};
	$dsinfo .=get_dataset_info(dsid=>$dsid, chr=>$chr) if $dsid;
	push @seq_nums, {
			  SEQ_NUM=>$i,
			 };
	my %opts = (
		    SEQ_NUM=>$i,
		    REV_YES=>$revy,
		    REV_NO=>$revn,
		    REF_YES=>$refy,
		    REF_NO=>$refn,
		    DRUP=>$drup,
		    DRDOWN=>$drdown,
		    DRACCN=>$draccn,
		    DSID=>$dsid,
		    GBACCN=>$gbaccn,
		    GBSTART=>$gbstart,
		    GBLENGTH=>$gblength,
		    POS=>$pos,
		    DSINFO=>$dsinfo,
		    EXON_MASK_ON=>$maskexonon,
		    EXON_MASK_OFF=>$maskexonoff,

		   );
#	print STDERR "pos $i: $pos\n";
	$opts{COGEPOS} = qq{<option value="cogepos$i" selected="selected">CoGe Database Position</option>} if $pos;
	push @seq_sub, {%opts}
	  
      }
    #page preferences
    my $pad_gs = $form->param("pad_gs") if $form->param("pad_gs");
    $pad_gs = 0 unless $pad_gs;
    my $prog = get_opt(params=>$prefs, form=>$form, param=>'prog');
    $prog = "blastz" unless $prog;
    my $image_width = get_opt(params=>$prefs, form=>$form, param=>'iw');
    $image_width = 1000 unless $image_width;
    my $image_height = get_opt(params=>$prefs, form=>$form, param=>'ih');
    $image_height = 100 unless $image_height;
    my $feature_height = get_opt(params=>$prefs, form=>$form, param=>'fh');
    $feature_height = 10 unless $feature_height;
    my $padding  = get_opt(params=>$prefs, form=>$form, param=>'padding');
    $padding = 2 unless defined $padding;
    my $gc_color = get_opt(params=>$prefs, form=>$form, param=>'gc');
    $gc_color = 0 unless $gc_color;
    my $nt_color = get_opt(params=>$prefs, form=>$form, param=>'nt');
    $nt_color = 1 unless $nt_color;
    my $auto_adjust_feats = get_opt(params=>$prefs, form=>$form, param=>'overlap');
    $auto_adjust_feats = 0 unless defined $auto_adjust_feats;
    my $hiqual = get_opt(params=>$prefs, form=>$form, param=>'hiqual');
    $hiqual = 0 unless $hiqual;
    my $color_hsp = get_opt(params=>$prefs, form=>$form, param=>'colorhsp');
    $color_hsp = 0 unless $color_hsp;
    my $hsp_label = get_opt(params=>$prefs, form=>$form, param=>'hsplabel');
    $hsp_label = undef unless defined $hsp_label;
    my $hsp_limit = get_opt(params=>$prefs, form=>$form, param=>'hsplim');
    $hsp_limit = 0 unless $hsp_limit;
    my $hsp_limit_num = get_opt(params=>$prefs, form=>$form, param=>'hsplimnum');
    $hsp_limit_num = 20 unless defined $hsp_limit_num;
    my $draw_model = get_opt(params=>$prefs, form=>$form, param=>'draw_model');
    $draw_model = "full" unless $draw_model;
    my $hsp_overlap_limit = get_opt(params=>$prefs, form=>$form, param=>'hsp_overlap_limit');
    $hsp_overlap_limit = 0 unless $hsp_overlap_limit;

    my $template = HTML::Template->new(filename=>'/opt/apache/CoGe/tmpl/GEvo.tmpl');
    $template->param(PAD_GS=>$pad_gs);
    $template->param(IMAGE_WIDTH=>$image_width);
    $template->param(IMAGE_HEIGHT=>$image_height);
    $template->param(FEAT_HEIGHT=>$feature_height);
    $template->param(PADDING=>$padding);
    $template->param(HSP_OVERLAP_LIMIT=>$hsp_overlap_limit);
    if ($gc_color) {$template->param(GC_COLOR_YES=>"checked");}
    else {$template->param(GC_COLOR_NO=>"checked");}
    if ($nt_color) {$template->param(NT_COLOR_YES=>"checked");}
    else {$template->param(NT_COLOR_NO=>"checked");}
    if ($auto_adjust_feats) {$template->param(OVERLAP_YES=>"checked");}
    else {$template->param(OVERLAP_NO=>"checked");}
    if ($hiqual) {$template->param(HIQUAL_YES=>"checked");}
    else {$template->param(HIQUAL_NO=>"checked");}
    if ($color_hsp) {$template->param(COLOR_HSP_YES=>"checked");}
    else {$template->param(COLOR_HSP_NO=>"checked");}
    if ($hsp_label eq "staggered") {$template->param(HSP_LABELS_STAG=>"selected");}
    elsif ($hsp_label eq "linear") {$template->param(HSP_LABELS_LIN=>"selected");}
    else {$template->param(HSP_LABELS_NO=>"selected");}
    if ($hsp_limit) {$template->param(HSP_LIMIT_YES=>"checked");}
    else {$template->param(HSP_LIMIT_NO=>"checked");}
    $template->param(HSP_LIMIT_NUM=>$hsp_limit_num);
    if ($draw_model eq "full") {$template->param(DRAW_MODEL_FULL=>"selected");}
    elsif ($draw_model eq "gene") {$template->param(DRAW_MODEL_GENE=>"selected");}
    elsif ($draw_model eq "mRNA") {$template->param(DRAW_MODEL_mRNA=>"selected");}
    elsif ($draw_model eq "CDS") {$template->param(DRAW_MODEL_CDS=>"selected");}
    elsif ($draw_model eq "RNA") {$template->param(DRAW_MODEL_RNA=>"selected");}
    else {$template->param(DRAW_MODEL_NO=>"selected");}

    my $box = HTML::Template->new(filename=>'/opt/apache/CoGe/tmpl/box.tmpl');
    #generate sequence submission selector
    $template->param(SEQ_SELECT=>1);
    $template->param(SEQ_SELECT_LOOP=>\@seq_sub);
    my $seq_submission = $template->output;
    $template->param(SEQ_SELECT=>0);

    #generate the hsp color option
    my $hsp_colors = gen_hsp_colors(num_seqs=>$num_seqs, prefs=>$prefs);
    my $html;
    my $spike_len = spike_filter_select();

    $template->param(SPIKE_LEN=>$spike_len);
    $template->param(SEQ_RETRIEVAL=>1);
    $template->param(NUM_SEQS=>$num_seqs);
    $message .= "<BR/>" if $message;
    $template->param(MESSAGE=>$message);
    $template->param(SEQ_SUB=>$seq_submission);
    $template->param(HSP_COLOR=>$hsp_colors);
    $template->param(GO_RUN=>gen_go_run($num_seqs));
    $template->param(AUTOSEARCH=>$autosearch_string);
    my $gobe_version = `svnversion /opt/apache/CoGe/gobe/flash`;
    $gobe_version =~ s/\n//g;;
    $template->param(GOBE_VERSION=>$gobe_version);
    $box->param(BOX_NAME=>"Options:");
    $template->param(OPTIONS=>1);
    $template->param(ALIGNMENT_PROGRAMS=>algorithm_list($prog));
    $template->param(SAVE_SETTINGS=>gen_save_settings($num_seqs)) unless !$USER || $USER =~ /public/i;
    $box->param(BODY=>$template->output);
    $html .= $box->output;
    return $html;
  }

sub run
  {
    my %opts = @_;
    my $num_seqs = $opts{num_seqs} || $NUM_SEQS;
    my $spike_len = $opts{spike};
    my $iw = $opts{iw};
    my $ih = $opts{ih};
    my $feat_h = $opts{fh};
    my $show_gc = $opts{gc};
    my $show_nt = $opts{nt};
    my $color_hsp = $opts{colorhsp};
    my $hsp_label = $opts{hsplabel};
    my $draw_model = $opts{draw_model};
    my $hsp_overlap_limit = $opts{hsp_overlap_limit};
    my $overlap_adjustment = $opts{overlap};
    my $hiqual = $opts{hiqual};
    my $hsp_limit = $opts{hsplim};
    my $hsp_limit_num = $opts{hsplimnum};
    my $show_hsps_with_stop_codon = $opts{showallhsps};
    my $padding = $opts{padding};
    my ($analysis_program, $param_string, $parser_opts) = get_algorithm_options(%opts);
    my $basefilename = $opts{basefile};
    my $show_spike = $opts{show_spike} || 0;
    my $pad_gs = $opts{pad_gs} || 0;
    my $color_overlapped_features = $opts{color_overlapped_features};
    my $hsp_overlap_length = $opts{hsp_overlap_length};
    my $message;
    $cogeweb = initialize_basefile(basename=>$basefilename, prog=>"GEvo");
    my @hsp_colors;
    for (my $i = 1; $i <= num_colors($num_seqs); $i++)
      {
	my $r = $opts{"r$i"};
	my $g = $opts{"g$i"};
	my $b = $opts{"b$i"};
	my @tmp;
	foreach my $color ($r, $g, $b)
	  {
	    $color = 0 unless $color && $color =~ /^\d+$/;
	    $color = 0 if $color < 0;
	    $color = 255 if $color > 255;
	    push @tmp, $color;
	  }
	push @hsp_colors,\@tmp;
      }



    my $stagger_label = $hsp_label && $hsp_label =~ /staggered/i ? 1 : 0;
    my $feature_labels = !$hsp_label ? 0 : 1;
    my $form = $FORM;
    my $gevo_link = $form->url."?prog=$analysis_program";
    $gevo_link .= ";spike_len=$spike_len";
    
    my @sets;
    my $html;
    my $t1 = new Benchmark;
    my $spike_seq;
    my $seqcount = 1;
    for (my $i = 1; $i <= $num_seqs; $i++)
      {
	my $skip_seq =$opts{"skip_seq$i"};
	next if $skip_seq;
	my $accn = $opts{"draccn$i"};
	my $featid = $opts{"featid$i"};
	my $feat = $coge->resultset('Feature')->find($featid) if $featid;
	my $dsid;
	$dsid = $feat->dataset_id if $feat;
	$dsid = $opts{"dsid$i"} unless $dsid;
	$dsid = $opts{"posdsid$i"} unless $dsid;
	my $chr;
	$chr = $feat->chromosome if $feat;
	$chr = $opts{"chr$i"} unless $chr;

	my $drup = $opts{"drup$i"};
	my $drdown = $opts{"drdown$i"};
	$drup += $pad_gs;
	$drdown += $pad_gs;
	my $pos = $opts{"pos$i"};
	my $gbaccn = $opts{"gbaccn$i"};
	my $gbstart = $opts{"gbstart$i"};
	$gbstart = 1 unless defined $gbstart;
	my $gblength = $opts{"gblength$i"};

	my $dirseq = $opts{"dirseq$i"};
	my $dirstart = $opts{"dirstart$i"};
	my $dirlength = $opts{"dirlength$i"};

	my $rev = $opts{"rev$i"};
	my $mask_cds_flag = $opts{"maskcds$i"};
	my $mask_ncs_flag = $opts{"maskncs$i"};
	my ($up, $down, $seq);
	my ($file, $obj);
	my $reference_seq =$opts{"ref_seq$i"};
	my $repeat_mask =$opts{"repmask$i"};
	next unless $accn || $featid || $gbaccn || $dirseq|| $pos;
	$gevo_link .= ";accn$seqcount=".CGI::escape($accn) if $accn;
	$gevo_link .= ";x$seqcount=".CGI::escape($pos) if $pos;
	$gevo_link .= ";fid$seqcount=".CGI::escape($featid) if $featid;
	$gevo_link .= ";dsid$seqcount=".CGI::escape($dsid) if $dsid;
	$gevo_link .= ";chr$seqcount=".CGI::escape($chr) if $chr;
	$gevo_link .= ";dr$seqcount"."up=$drup" if defined $drup;
	$gevo_link .= ";dr$seqcount"."down=$drdown" if defined $drdown;
	$gevo_link .= ";gbaccn$seqcount=".CGI::escape($gbaccn) if $gbaccn;
	$gevo_link .= ";gbstart$seqcount=$gbstart" if $gbstart;
	$gevo_link .= ";gblength$seqcount=$gblength" if $gblength;
	$gevo_link .= ";rev$seqcount=1" if $rev;
	$gevo_link .= ";nref$seqcount=1" unless $reference_seq;
	$seqcount++;
#	print STDERR "pos: $pos, dsid: $dsid, chr: $chr\n";
	if ($featid || $pos)
	  {
	    $obj = get_obj_from_genome_db( accn=>$accn, featid=>$featid, pos=>$pos, dsid=>$dsid, rev=>$rev, up=>$drup, down=>$drdown, chr=>$chr );
	    if ($obj)
	      {
		($file, $spike_seq, $seq) = 
		  generate_seq_file(obj=>$obj,
				    mask_cds=>$mask_cds_flag,
				    mask_ncs=>$mask_ncs_flag,
				    spike_len=>$spike_len,
				    seq_num=>$i,
				    repeat_mask=>$repeat_mask,
				   );
		$up = $drup;
		$down = $drdown;
	      }
	    else
	      {
		$message .=  "Unable to generate sequence for sequence $i.  (Skipped.)";
	      }
	  }
 	elsif ($dirseq )
 	  {
 	    ($obj) = generate_obj_from_seq($dirseq, $i, $rev);
	    $dirlength = length($obj->sequence)-$dirstart+1 unless $dirlength;
 	    if ($obj)
	      {
		#add an anchor
		$obj->add_feature(
			      type=>"direct sequence submission",
			      location=>(1-$dirstart+1)."..".(length($obj->sequence)-$dirstart+1),
			      strand=>1,
			      qualifiers=>{
					   type=>"anchor",
					   names=>[$obj->accn],
					  }
			     );

		($file, $spike_seq, $seq) = 
		  generate_seq_file (
				     obj=>$obj,
				     mask_cds=>$mask_cds_flag,
				     mask_ncs=>$mask_ncs_flag,
				     startpos=>$dirstart,
				     length=>$dirlength,
				     spike_len=>$spike_len, 
				     seq_num=>$i,
				     repeat_mask=>$repeat_mask,
				    );
		$obj->start($dirstart);
		$obj->stop($dirstart+length($seq)-1);
		$obj->chromosome(1);
#		$obj->dataset("NA");
		$up = $dirstart;
		$down = $dirlength;
	      }
	    else
	      {
		$message .= "Problem with direct sequence submission\n";
	      }
	  }
	elsif ($gbaccn )
	  {
	    $obj = new CoGe::Accessory::GenBank;
	    $obj->add_gene_models(1); #may want to make a user selectable option
 	    my ($res, $error) = $obj->get_genbank_from_ncbi($gbaccn, $rev);
	    $message .= $error."\n" unless $res;
	    
	    if ($obj->accn)
	      {
		#add an anchor
		$obj->add_feature(
				  type=>"genbank entry",
				  location=>(1-$gbstart+1)."..".(length($obj->sequence)-$dirstart+1),
				  strand=>1,
				  qualifiers=>{
					       type=>"anchor",
					       names=>[$gbaccn],
					      }
				 );

		($file, $spike_seq, $seq) = 
		  generate_seq_file (
				     obj=>$obj,
				     mask_cds=>$mask_cds_flag,
				     mask_ncs=>$mask_ncs_flag,
				     startpos=>$gbstart,
				     length=>$gblength,
				     spike_len=>$spike_len,
				     seq_num=>$i,
				     repeat_mask=>$repeat_mask,
				    );
		$obj->start($gbstart);
		$obj->stop($gbstart + length($seq)-1);
		$obj->chromosome("?") unless $obj->chromosome;
		$up = $gbstart;
		$down = $gblength;
	      }
	    else
	      {
		$message .= "No GenBank entry found for $gbaccn\n";
	      }
 	  }
	next unless $obj;
	unless ($show_spike)
	  {
	    $seq =~ s/N*$spike_seq$//;
	    $seq =~ s/N*$spike_seq$//;
	  }
	$obj->sequence($seq);

	if ($obj && $obj->sequence)
	  {
	    #need to check for duplicate accession names -- sometimes happens and major pain in the ass for other parts of the code
	    my $accn = $obj->accn;
	    my $count = 0;
	    foreach my $accn2 (map {$_->{obj}->accn()} @sets)
	      {
		$accn2 =~ s/\*\*\d+\*\*$//;
		$count++ if $accn eq $accn2;
	      }
	    $accn .= "**".($count+1)."**" if $count;
	    $obj->accn($accn);
	    push @sets, {
			 obj=>$obj,
			 file=>$file,
#			 file_begin=>$file_begin,
#			 file_end=>$file_end,
			 accn=>$accn,
			 rev=>$rev,
			 up=>$up,
			 down=>$down,
			 spike_seq=>$spike_seq,
			 reference_seq=>$reference_seq,
			 seq_num=>$i,
			};
	  }
	else
	  {
#	    push @sets, {seq_num=>$i};
	  }
      }
    $seqcount--;

    $gevo_link .= ";num_seqs=".$seqcount;
    $gevo_link .= ";hsp_overlap_limit=".$hsp_overlap_limit if defined $hsp_overlap_limit;
    unless (@sets >1)
      {
	$message .= "Problem retrieving information.  Please check submissions.\n";
	return '', '', '', '',0,$message;
      }
    my $t2 = new Benchmark;
    # set up output page
    
    # run bl2seq
    
    my $analysis_reports;

    if ($analysis_program eq "blastz")
      {
	$analysis_reports = run_blastz(sets=>\@sets, params=>$param_string, parser_opts=>$parser_opts);
      }
    elsif ($analysis_program eq "LAGAN")
      {
	$analysis_reports = run_lagan (sets=>\@sets, params=>$param_string, parser_opts=>$parser_opts);
      }
    elsif ($analysis_program eq "CHAOS")
      {
	$analysis_reports = run_chaos (sets=>\@sets, params=>$param_string, parser_opts=>$parser_opts);
      }
    elsif ($analysis_program eq "DiAlign_2")
      {
	($analysis_reports, $analysis_program,$message) = run_dialign (sets=>\@sets, params=>$param_string, parser_opts=>$parser_opts);
      }
    else
      {
	$analysis_reports = run_bl2seq(sets=>\@sets, params=>$param_string, parser_opts=>$parser_opts, blast_program=>$analysis_program, spike_seq=>$spike_seq);
      }
    $analysis_reports = [] unless ref($analysis_reports) =~ /ARRAY/i;
    write_log($message, $cogeweb->logfile) if $message;
   #sets => array or data for blast
   #blast_reports => array of arrays (report, accn1, accn2, parsed data)
    my $t3 = new Benchmark;
    my $count = 1;
    my $frame_height;
    initialize_sqlite();
    my @gfxs;
    foreach my $item (@sets)
      {
	my $obj = $item->{obj};
	next unless $obj->sequence;
	my $rev = $item->{rev};
	my $spike_seq = $item->{spike_seq};
	if ($obj)
	  {
	    write_log("generating image ($count/".scalar @sets.")for ".$obj->accn, $cogeweb->logfile);
	    $count++;
#	    my ($image, $map, $mapname, $gfx, $eval_cutoff, $feat_count, $overlap_count) = 
	    my ($gfx, $stats) = 
	      generate_image(
			     set=>$item,
			     start=>1,
			     stop=>length($obj->sequence),
			     data=>$analysis_reports,
			     iw=>$iw,
			     ih=>$ih,
			     fh=>$feat_h,
			     show_gc=>$show_gc,
			     show_nt=>$show_nt,
			     stagger_label=>$stagger_label,
			     overlap_adjustment=>$overlap_adjustment,
			     feature_labels=>$feature_labels,
			     hsp_limit=>$hsp_limit,
			     hsp_limit_num=>$hsp_limit_num,
			     color_hsp=>$color_hsp,
			     hsp_colors=>\@hsp_colors,
			     spike_sequence=>$spike_seq,
			     show_hsps_with_stop_codon => $show_hsps_with_stop_codon,
			     hiqual=>$hiqual,
			     padding=>$padding,
#			     seq_num=>$item->{seq_num},
			     reverse_image=>$rev,
			     draw_model=>$draw_model,
			     hsp_overlap_limit=>$hsp_overlap_limit,
			     color_overlapped_features=>$color_overlapped_features,
			     hsp_overlap_length=>$hsp_overlap_length,
			    );
#	    $item->{feat_count} = $feat_count;
#	    $item->{overlap_count} = $overlap_count;
	    my $filename = $cogeweb->basefile."_".$item->{seq_num}.".png";
	    $filename = check_filename_taint($filename);
	    my $image = basename($filename);
	    $item->{image}=$image;
	    $item->{stats}=$stats;
	    $gfx->generate_png(file=>$filename);
	    generate_image_db(set=>$item, gfx=>$gfx);
	    $frame_height += $gfx->ih + $gfx->ih*.1;
	  }
      }
    my $t3_5 = new Benchmark;
    image_db_create_hsp_pairs();

    my $t4 = new Benchmark;
    $html .= qq{<DIV id=flash_viewer></DIV>};
    $html .= qq{<table>};
    $html .= qq{<tr valign=top><td class = small>Alignment reports};
    my $stats_file = $cogeweb->basefile."_stats.txt";
    $stats_file = check_filename_taint($stats_file);

    if ($analysis_reports && @$analysis_reports)
      {
	foreach my $item (@$analysis_reports)
	  {
	    my $report = $item->[0];
	    my $accn1 = $item->[1];
	    my $accn2 = $item->[2];
	    my $basereportname = basename( $report );
	    $basereportname = $TEMPURL . "/$basereportname\n";
	    $html .= "<div><font class=small><A HREF=\"$basereportname\" target=_new>View alignment output for $accn1 versus $accn2</A></font></DIV>\n";
	  }
      }
    else
      {
	$html .= "<div class=small>No alignment reports were generated</DIV>\n";
      }
    $html .= qq{<td class = small>Fasta files};
    foreach my $item (@sets)
      {
	next unless $item->{file};
	my $basename = $TEMPURL."/".basename ($item->{file});
	print STDERR "basename is undefined: $basename\n" if $basename =~ /defined/i;
	my $accn = $item->{accn};
	$html .= "<div><font class=small><A HREF=\"$basename\" target=_new>Fasta file for $accn</A></font></DIV>\n";
      }
    $html .= qq{<td class = small><a href = "http://baboon.math.berkeley.edu/mavid/gaf.html">GAF</a> annotation files};
    foreach my $item (@sets)
      {
	my $anno_file = generate_annotation(%$item);
	next unless $anno_file;
	my $basename = $TEMPURL."/".basename ($anno_file);
	my $accn = $item->{accn};
	$html .= "<div><font class=small><A HREF=\"$basename\" target=_new>Annotation file for $accn</A></font></DIV>\n";
     }
    $html .= qq{<td class = small>SQLite db};
    my $dbname = $TEMPURL."/".basename($cogeweb->sqlitefile);
    
    $html .= "<div class=small><A HREF=\"$dbname\" target=_new>SQLite DB file</A></DIV>\n";
    $html .= qq{<td class = small>Log File};
    my $logfile = $TEMPURL."/".basename($cogeweb->logfile);
    $html .= "<div class=small><A HREF=\"$logfile\" target=_new>Log</A></DIV>\n";
    my $tiny = get("http://tinyurl.com/create.php?url=$gevo_link");
    ($tiny) = $tiny =~ /<b>(http:\/\/tinyurl.com\/\w+)<\/b>/;
    $html .= qq{<td class = small>GEvo Link<div class=small><a href=$tiny target=_new>$tiny<br>(See log file for full link)</a></div>};
    $html .=qq{<a href="GEvo_direct.pl?name=$basefilename" target=_new>Results only</a>};
    $html .= qq{<td class = small>Gene Homology Stats:<table>};


    write_log("Summary of overlap data for genes only:", $stats_file);
    write_log(join ("\t", qw(REGION1 REGION2 LENGTH OVERLAP/TOTAL PERCENT_OVERLAP)),$stats_file);
    foreach my $item (@sets)
      {
	$html .= "<tr class=small>";
	my @data = ("<td>".$item->{obj}->accn,commify(length $item->{obj}->sequence)." bp",$item->{stats}{feat_counts}{gene}{overlap}."/".$item->{stats}{feat_counts}{gene}{count});
	push @data, sprintf("%.2f", $item->{stats}{feat_counts}{gene}{overlap}/$item->{stats}{feat_counts}{gene}{count}*100)."%" if $item->{stats}{feat_counts}{gene}{count};
	$html .= join ("<td nowrap>",@data);
	shift @data;
	write_log(join ("\t",$item->{obj}->accn,"all",@data), $stats_file);
      }

    write_log("\n",$stats_file);
    write_log("Features overlapped by HSPs:", $stats_file);
    write_log(join ("\t", qw(REGION1 REGION2 TYPE OVERLAP_COUNT TOTAL PERCENT_W_OVERLAP)),$stats_file);
    foreach my $item (@sets)
      {
	my $stats = $item->{stats};
	foreach my $type (sort keys %{$stats->{feat_counts}})
	  {
	    foreach my $accn1 (sort keys %{$stats->{data}})
	      {
		foreach my $accn2 (sort keys %{$stats->{data}{$accn1}})
		  {
		    next if $type eq "hsp";
		    my $overlap_count = $stats->{data}{$accn1}{$accn2}{$type}{overlapped_hsps} || 0;
		    my $percent = sprintf("%.2f",$overlap_count/$stats->{feat_counts}{$type}{count}*100) if $stats->{feat_counts}{$type}{count};
		    write_log(join ("\t", $accn1, $accn2, $type, $overlap_count, $stats->{feat_counts}{$type}{count}, $percent), $stats_file);
		  }
	      }
	  }
      }


    write_log("\n",$stats_file);
    write_log("Average percent identity for all HSPs (normalized to HSP length):", $stats_file);
    write_log(join ("\t", qw(REGION1 REGION2 AVG_PERCENT_ID)),$stats_file);
    foreach my $item (@sets)
      {
	my $stats = $item->{stats};
	foreach my $accn1 (sort keys %{$stats->{data}})
	  {
	    foreach my $accn2 (sort keys %{$stats->{data}{$accn1}})
	      {
		my $len= 0;
		map {$len+=$_->[1]} @{$stats->{data}{$accn1}{$accn2}{hsp}};
		my $pid = 0;
		map {$pid+=$_->[1]*$_->[0]} @{$stats->{data}{$accn1}{$accn2}{hsp}};
		write_log("$accn1\t$accn2\t".sprintf("%.2f",$pid/$len), $stats_file);
	      }
	  }
      }
    my $stats_url = $TEMPURL."/".basename($stats_file);
    $html .= qq{<tr class=small><td><a href=$stats_url target=_new>Stats file</a>};
    $html .= qq{</table>};
    $html .= qq{</table>};

    my $template = HTML::Template->new(filename=>'/opt/apache/CoGe/tmpl/box.tmpl');
    my $results_name = "Results: $analysis_program";
    $results_name .= qq{ <span class="small">(spike sequence filter length: $spike_len)</span>} if $spike_len;
    $template->param(BOX_NAME=>$results_name);
    $template->param(BODY=>$html);
    my $outhtml = $template->output;
    my $t5 = new Benchmark;
    my $db_time = timestr(timediff($t2,$t1));
    my $blast_time = timestr(timediff($t3,$t2));
    my $image_time = timestr(timediff($t3_5,$t3));
    my $db_hsp_time = timestr(timediff($t4,$t3_5));
    my $html_time = timestr(timediff($t5,$t4));
    my $bench =  qq{
GEvo Benchmark: $DATE
Time to get sequence                              : $db_time
Time to run $analysis_program                     : $blast_time
Time to generate images, maps, and sqlite database: $image_time
Time to find and update sqlite database for HSPs  : $db_hsp_time
Time to process html                              : $html_time
};
    print STDERR $bench if $BENCHMARK;
    write_log($bench, $cogeweb->logfile);
    write_log("Finished!", $cogeweb->logfile);
    write_log("\nGEvo link: $gevo_link\n", $cogeweb->logfile);
    write_log("Tiny url: $tiny", $cogeweb->logfile);
    $count--;
    return $outhtml, $iw+400, $frame_height, $cogeweb->basefilename,$count,$message;

}


sub generate_image
  {
    my %opts = @_;
    my $set = $opts{set};
    my $gbobj = $set->{obj};
    my $start = $opts{start};
    my $stop = $opts{stop};
    my $data = $opts{data};
    my $masked_exons = $opts{mask};
    my $masked_ncs = $opts{mask_ncs};
    my $spike_seq = $opts{spike_sequence};
    my $iw = $opts{iw} || 1600;
    my $ih = $opts{ih} || 200;
    my $fh = $opts{fh} || 25;
    my $show_gc = $opts{show_gc};
    my $show_nt = $opts{show_nt};
    my $reverse_image = $opts{reverse_image};
    my $stagger_label = $opts{stagger_label};
    my $overlap_adjustment = $opts{overlap_adjustment};
    my $feature_labels = $opts{feature_labels};
    my $hsp_limit = $opts{hsp_limit};
    my $color_hsp = $opts{color_hsp};
    my $hsp_limit_num = $opts{hsp_limit_num};
    my $eval_cutoff = $opts{eval_cutoff};
    my $hsp_colors = $opts{hsp_colors};
    my $show_hsps_with_stop_codon = $opts{show_hsps_with_stop_codon};
    my $hiqual = $opts{hiqual};
    my $padding = $opts{padding} || 5;
#    my $seq_num = $opts{seq_num};
    my $draw_model = $opts{draw_model};
    my $hsp_overlap_limit = $opts{hsp_overlap_limit};
    my $color_overlapped_features = $opts{color_overlapped_features};
    my $hsp_overlap_length = $opts{hsp_overlap_length};
    my $graphic = new CoGe::Graphics;
    my $gfx = new CoGe::Graphics::Chromosome;
#    print STDERR "$start"."::"."$stop"." -- ".length($gbobj->sequence)."\n";
    $graphic->initialize_c (
			    c=>$gfx,
			    iw=>$iw,
			    start=> $start,
			    stop => $stop,
			    draw_chr=>1,
			    draw_ruler=>1,
			    draw_chr_end=>0,
			    chr_start_height=>$ih,
			    chr_mag_height=>5,
			    feature_start_height=>$fh,
			    mag=>0,
			    mag_off=>1,
			    chr_length => length($gbobj->sequence),
			    fill_labels=>1,
			    forcefit=>1,
			    minor_tick_labels=>1,
			    feature_labels=>$feature_labels,
			    draw_hi_qual=>$hiqual,
			    padding=>$padding,
			   );
 #   print STDERR $gfx->_region_start ."--".$gfx->_region_stop."\n";
    $gfx->overlap_adjustment(1);
    $gfx->top_padding(15);
    $gfx->skip_duplicate_features(1);
    $gfx->DEBUG(0);
    $gfx->major_tick_labels(0);
    $graphic->process_nucleotides(c=>$gfx, seq=>$gbobj->sequence, layers=>{gc=>$show_gc, nt=>$show_nt});
    my $stats = process_hsps(
			  c=>$gfx, 
			  data=>$data, 
			  accn=>$gbobj->accn, 
			  rev=>$reverse_image, 
			  seq_length=> length($gbobj->sequence), 
			  stagger_label=>$stagger_label, 
			  hsp_limit=>$hsp_limit,
			  hsp_limit_num=>$hsp_limit_num,
			  gbobj=>$gbobj,
			  spike_seq=>$spike_seq,
			  color_hsp=>$color_hsp,
			  colors=>$hsp_colors,
			  show_hsps_with_stop_codon=>$show_hsps_with_stop_codon,
			  hsp_overlap_limit=>$hsp_overlap_limit,
			  hsp_overlap_length=>$hsp_overlap_length,
			 );
    my ($feat_counts) = process_features(c=>$gfx, obj=>$gbobj, start=>$start, stop=>$stop, overlap_adjustment=>$overlap_adjustment, draw_model=>$draw_model, color_overlapped_features=>$color_overlapped_features);
    $stats->{feat_counts} = $feat_counts;
    return ($gfx, $stats);
  }

sub initialize_sqlite
  {
    my %opts = @_;
    my $dbfile = $cogeweb->sqlitefile;
    return if -r $dbfile;
    my $dbh = DBI->connect("dbi:SQLite:dbname=$dbfile","","");
    my $create = qq{
CREATE TABLE image_data
(
id INTEGER PRIMARY KEY,
name varchar,
type varchar,
xmin integer,
xmax integer,
ymin integer,
ymax integer,
bpmin integer,
bpmax integer,
image_track integer,
image_id integer,
pair_id integer,
link varchar,
annotation text,
color varchar
)
};
    $dbh->do($create);
    $dbh->do('CREATE INDEX type ON image_data (type)');
    $dbh->do('CREATE INDEX image_id ON image_data (image_id)');
    $dbh->do('CREATE INDEX x_idx ON image_data (xmin,xmax)');
    $dbh->do('CREATE INDEX y_idx ON image_data (ymin,ymax)');
    $dbh->do('CREATE INDEX image_track ON image_data (image_track)');
    $dbh->do('CREATE INDEX pair_id ON image_data (pair_id)');

#id INTEGER PRIMARY KEY AUTOINCREMENT,
    $create = qq{
CREATE TABLE image_info
(
id INTEGER,
iname varchar,
title varchar,
px_width integer,
bpmin integer,
bpmax integer,
dsid integer,
chromosome varchar
)
};
#TODO: make sure to populate the bpmin and pbmax and image width!!!!!!!!!!!!!!!!!
    $dbh->do($create);
    $dbh->do('CREATE INDEX id ON image_info (id)');
    system "chmod +rw $dbfile";
  }

sub generate_image_db
  {
    my %args = @_;
    my $gfx = $args{gfx};
    my $set = $args{set};
    next unless $gfx;
    my $dbh = DBI->connect("dbi:SQLite:dbname=".$cogeweb->sqlitefile,"","");
    my $image = $set->{image};
    my $accn = $set->{accn};
    my $title;
    $title = $accn;
    $title .= " ".$set->{obj}->organism() if $set->{obj}->organism();
    $title .= "(chr: ".$set->{obj}->chromosome." ". $set->{obj}->start."-".$set->{obj}->stop.")" if defined $set->{up};
    $title .= qq! Reverse Complement! if $set->{rev};
    my $width = $gfx->image_width;
    my $dsid = $set->{obj}->dataset;
    my ($chr) =  $set->{obj}->chromosome =~ /(\d+)/;
    $chr = "NULL" unless $chr;
    $dsid = "NULL" unless $dsid;
    my $image_start = $set->{obj}->start;
    my $image_stop = $set->{obj}->stop;
    my $image_id = $set->{seq_num};
    my $statement = qq{
INSERT INTO image_info (id, iname, title, px_width,dsid, chromosome, bpmin, bpmax) values ($image_id, "$image", "$title", $width, "$dsid", $chr, $image_start, $image_stop)
};
    print STDERR $statement unless $dbh->do($statement);
    foreach my $feat ($gfx->get_feats)
      {
	if ($feat->fill)
	  {
	    next unless $feat->type eq "anchor";
	  }
	next unless $feat->image_coordinates;
	my $type = $feat->type;
	my $pair_id = "-99";
	my $coords = $feat->image_coordinates;
	$coords =~ s/\s//g;
	next if $feat->type eq "unknown";
	my $name = $feat->type =~ /HSP/i ? $feat->alt : $feat->label;
	$name .= "_".$feat->type;# unless $name;

	my $color = "NULL";
	if ($feat->type =~ /HSP/)
	  {
	    $color ="#";
	    foreach my $c (@{$feat->color})
	      {
		$color .= sprintf("%X",$c);
	      }
	  }
	#generate link
	my $link = $feat->link;
	$link = " " unless $link;
	$link =~ s/'//g if $link;
	#generate image track
	my $image_track = $feat->track;
	$image_track = "-".$image_track if $feat->strand =~ /-/;
	
	my ($xmin, $ymin, $xmax, $ymax) = split /,/, $coords;
	$xmin++;
	$xmax++;
	my $anno = $feat->description;
	#	    $anno =~ s/'|"//g;
	$anno =~ s/'//g if $anno;
	$anno =~ s/<br\/?>/&#10;/ig if $anno;
	$anno =~ s/\n/&#10;/g if $anno;
	$anno =~ s/[\[\]\(\)]/ /g if $anno;
	$anno = " " unless $anno;
	my $start = $feat->start;
	my $stop = $feat->stop;
	my $length_nt = $stop-$start+1;
	my $length_pix = $xmax-$xmin;
	next if !$feat->{anchor} && ($length_nt == 0 );
	$type = "anchor" if $feat->{anchor};
	my $bpmin = $set->{obj}->start+$start-1;
	my $bpmax = $set->{obj}->start+$stop-1;
	$statement = qq{
INSERT INTO image_data (name, type, xmin, xmax, ymin, ymax, bpmin,bpmax,image_id, image_track,pair_id, link, annotation, color) values ("$name", "$type", $xmin, $xmax, $ymin, $ymax, $bpmin, $bpmax,$image_id, "$image_track",$pair_id, '$link', '$anno', '$color')
};
	print STDERR $statement unless $dbh->do($statement);
      }
  }

sub image_db_create_hsp_pairs
  {
    my $dbh = DBI->connect("dbi:SQLite:dbname=".$cogeweb->sqlitefile,"","");
    my $query = qq{select id,name from image_data where type = "HSP"};
    my %ids;
    foreach my $item (@{$dbh->selectall_arrayref($query)})
      {
	push @{$ids{$item->[1]}},$item->[0];
      }
    foreach my $name (keys %ids)
      {
	if (@{$ids{$name}} != 2)
	  {
	    print STDERR "Problem with $name.  Retrieved", Dumper $ids{$name} if $DEBUG;
	    next;
	  }
	my ($id1, $id2) = @{$ids{$name}};
	my $statement = "update image_data set pair_id = $id1 where id = $id2";
	$dbh->do($statement);
	$statement = "update image_data set pair_id = $id2 where id = $id1";
	$dbh->do($statement);
      }
  }


sub process_features
  {
    #process features
    my %opts = @_;
    my $c = $opts{c};
    my $obj=$opts{obj};
    my $start=$opts{start};
    my $stop = $opts{stop};
    my $overlap = $opts{overlap_adjustment};
    my $draw_model = $opts{draw_model};
    my $color_overlapped_features = $opts{color_overlapped_features};
#    return unless $draw_model;
    my $accn = $obj->accn;
    my $track = 1;
    my %feat_counts;
    my @opts = ($start, $stop) if $start && $stop;
    unless (ref $obj)
      {
	warn "Possible problem with the object in process_features.  Returning";
	return 0;
      }
    my %prior_feat;
    foreach my $feat (sort {$a->start <=> $b->start} $obj->get_features())
      {
        my $f;
	my $type = $feat->type;
	my ($name) = sort { length ($b) <=> length ($a) || $a cmp $b} @{$feat->qualifiers->{names}} if ref ($feat->qualifiers) =~ /hash/i;
	my $anchor = $feat if ref ($feat->qualifiers) =~ /hash/i && $feat->qualifiers->{type} && $feat->qualifiers->{type} eq "anchor";
        if ($type =~ /pseudogene/i)
          {
	    next unless $draw_model eq "full";
	    $f = CoGe::Graphics::Feature::Gene->new();
	    $f->color([255,33,0,50]);
	    $f->order($track);
	    $f->overlay(1);
	    $f->mag(0.5);
          }
        elsif ($type =~ /Gene/i)
          {
	    next unless $draw_model eq "full" || $draw_model eq "gene" || $anchor;
	    $f = CoGe::Graphics::Feature::Gene->new();
	    #$f->color([255,0,0,50]);
        $f->color([219,219,219,50]);
	    $f->order($track);
	    $f->overlay(1);
	    $f->mag(0.5);
          }
        elsif ($type =~ /CDS/i)
          {
	    next unless $draw_model eq "full" || $draw_model eq "CDS";
	    $f = CoGe::Graphics::Feature::Gene->new();
	    $f->color([0,255,0, 50]);
	    $f->order($track);
	    $f->overlay(3);
	    if ($accn)
	      {
		foreach my $name (@{$feat->qualifiers->{names}})
		  {
            my $cleaned_name = $name;
            $cleaned_name =~ s/[\(\)]//g;
		    if ($accn =~ /^$cleaned_name\(?\d*\)?$/i)
		      {
			$f->color([255,255,0]) ;
			$f->label($name);
		      }
		  }
	      }
	    
		
          }
        elsif ($type =~ /mrna/i)
          {
	    next unless $draw_model eq "full" || $draw_model eq "mRNA";
        	$f = CoGe::Graphics::Feature::Gene->new();
        	$f->color([0,0,255, 50]);
        	$f->order($track);
		$f->overlay(2);
		$f->mag(0.75);
          }
        elsif ($type =~ /rna/i)
          {
	    next unless $draw_model eq "full" || $draw_model eq "RNA";
        	$f = CoGe::Graphics::Feature::Gene->new();
        	$f->color([200,200,200, 50]);
        	$f->order($track);
		$f->overlay(2);
		if ($accn)
		  {
		    foreach my $name (@{$feat->qualifiers->{names}})
		      {
			if ($accn =~ /^$name\(?\d*\)?$/i)
			  {
			    $f->color([255,255,0]) ;
			    $f->label($name);
			  }
		      }
		  }

          }
	elsif ($type =~ /anchor/)
	  {
	    $f = CoGe::Graphics::Feature::NucTide->new({type=>'anchor',start=>$feat->blocks->[0][0], strand=>1});
	    $f->color([0,0,255]);
	    $f->type($type);
	    $f->description($feat->annotation);
	    $c->add_feature($f);
	    $f = CoGe::Graphics::Feature::NucTide->new({type=>'anchor',start=>$feat->blocks->[0][0], strand=>-1});
	    $f->color([0,0,255]);
	    $f->type($type);
	    $f->description($feat->annotation);
	    $c->add_feature($f);
	    next;
	  }
	elsif ($type =~ /cns/i)
	  {
	    $f = CoGe::Graphics::Feature::HSP->new({start=>$feat->blocks->[0][0], stop=>$feat->blocks->[0][1]});
	    $f->color([255,51,0]);
	    $f->order($track);
	    $f->overlay(4);
	    $f->type($type);
	    $f->force_draw(1);
	    $c->add_feature($f);
	    next;
	  }
	#need to create an anchor if this feature is an anchor, but not to be drawn
	if ($anchor && !$f)
	  {
	    $f = CoGe::Graphics::Feature->new({type=>'anchor',start=>$feat->blocks->[0][0], stop=>$feat->blocks->[0][1], strand=>1});
	    $f->{anchor}=1;
	    $f->order(0);
	    $f->skip_overlap_search(1);
	    $f->description("auto generated anchor for GEvo");
	    $c->add_feature($f);
	    next;
	  }
        next unless $f;
	$f->color([255,0,255]) if $color_overlapped_features && $feat->qualifiers->{overlapped_hsp};
	my $strand = 1;
 	$strand = -1 if $feat->location =~ /complement/;
	foreach my $block (@{$feat->blocks})
	  {
	    $block->[0] =1 unless $block->[0]; #in case $block is set to 0
	    $f->add_segment(start=>$block->[0], stop=>$block->[1]);
	    $f->strand($strand);
	    print STDERR "\t", join ("-", @$block),"\n" if $DEBUG;
	  }

	print STDERR $name,"\n\n" if $DEBUG;
        $f->type($type);
	$f->description($feat->annotation);
	$f->link("FeatList.pl?fid=".$feat->qualifiers->{id}) if $feat->qualifiers->{id};
	my $foverlap = $overlap ? 0 : 1; #I think I need this to get this to work as expected
	$f->skip_overlap_search($foverlap);
	$f->{anchor}=1 if $anchor;
        $c->add_feature($f);

	if ($prior_feat{$feat->type}{$strand})
	  {
	    if ($feat->start < $prior_feat{$feat->type}{$strand}->stop)
	      {
		$prior_feat{$feat->type}{$strand} = $feat;
		next;
	      }
	  }
	next if $feat->start > $stop;
	next if $feat->stop < $start;
	$feat_counts{$feat->type}{count}++;
	$feat_counts{$feat->type}{overlap}++ if $feat->qualifiers->{overlapped_hsp};
	$prior_feat{$feat->type}{$strand} = $feat;
      }
    return (\%feat_counts);
  }

sub process_hsps
  {
    my %opts = @_;
    my $c = $opts{c};
    my $data = $opts{data};
    my $accn = $opts{accn};
    my $reverse = $opts{rev};
    my $seq_len = $opts{seq_length};
    my $stagger_label = $opts{stagger_label};
{'fraction_non_alphanumeric': 0.20222558844092287, 'fraction_numerical': 0.00908879049172687, 'mean_word_length': 3.597717775634844, 'pattern_counts': {'":': 2, '<': 169, '<?xml version=': 0, '>': 1578, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1007056', 'n_tokens_mistral': 35861, 'n_tokens_neox': 33589, 'n_words': 6229}