genome_canvas

#!/usr/bin/env perl

use strict;
use warnings;

### genome_canvas

use GenomeCanvas;
use GenomeCanvas::Band::TilingPath;
use GenomeCanvas::Band::Plot;
use GenomeCanvas::DensityBand::RepeatFeature;
#use GenomeCanvas::DensityBand::Gene;
use GenomeCanvas::Band::Gene;
use GenomeCanvas::Band::SNP;
use GenomeCanvas::Band::CpG;
use GenomeCanvas::Band::GeneticMap;
use GenomeCanvas::Band::FileFeatures;
use GenomeCanvas::Band::HomolFeature;
use GenomeCanvas::Band::GeneCluster;
use GenomeCanvas::Band::Scale;

use Bio::Otter::Lace::Defaults;
use Bio::Vega::DBSQL::DBAdaptor;

use Hum::Submission 'prepare_statement'; # For Sanger sequence names
use Hum::EnsCmdLineDB;

{
    my $scale           = 2000;
    my $v_tiles         = 0;
    my $h_tiles         = 0;
    my $chr_name        = '';
    my $chr_start       = 0;
    my $chr_end         = 0;
    my $font_size       = 12;
    my $sanger_names    = 0;
    my $landscape       = 1;
    my $portrait        = 0;
    my $gene_labels     = 1;

    my $show_gc_content = 1;
    my $show_key        = 1;
    my $show_genes      = 1;
    my $show_repeats    = 1;

    my $show_cpg        = 0;
    my $show_ruler      = 0;

    my $homologies      = 0;
    my $markers         = 0;
    my $clusters        = 0;
    my $snps            = 0;
    my $psdir           = '.';
    my $big_picture     = 0;
    my $otter           = 0;
    my $eps_file        = undef;
    my $span_file       = undef;
    my $dataset_name    = undef;
    my $relative_coords = 0;
    
    # Otter/Ensembl switch
    my $command_line_sub = \&Hum::EnsCmdLineDB::do_getopt;
    for (my $i = 0; $i < @ARGV; $i++) {
        if ($ARGV[$i] eq '-otter') {
            $otter = 1;
            splice(@ARGV, $i, 1);
            require Bio::Otter::Lace::Defaults;
            $command_line_sub = \&Bio::Otter::Lace::Defaults::do_getopt;
            last;
        }
    }   
    
    $command_line_sub->(
        'scale=s'           => \$scale,
        'htiles=i'          => \$h_tiles,
        'vtiles=i'          => \$v_tiles,
        'start=s'           => \$chr_start,
        'end=s'             => \$chr_end,
        'font_size=i'       => \$font_size,
        'sanger!'           => \$sanger_names,
        'landscape!'        => \$landscape,
        'portrait!'         => \$portrait,
        'gene_labels!'      => \$gene_labels,
        'genes!'            => \$show_genes,
        'gene_spans=s'      => \$span_file,
        'repeats!'          => \$show_repeats,
        'cpg!'              => \$show_cpg,
        'gc_content!'       => \$show_gc_content,
        'ruler!'            => \$show_ruler,
        'key!'              => \$show_key,
	    'homologies=s'      => \$homologies,
	    'markers=s'         => \$markers,
	    'clusters=s'        => \$clusters,
        'snps=s'            => \$snps,
        'psdir=s'           => \$psdir,
        'big!'              => \$big_picture,
        'eps=s'             => \$eps_file,
        'dataset=s'         => \$dataset_name,
        'chr|set=s'         => \$chr_name,
        'relative!'         => \$relative_coords,
        );
    
    my( $db, $ds, $cl );
    if ($otter) {
        unless ($dataset_name and $chr_name) {
            die "otter connection requires dataset and set";
        }
        $cl = Bio::Otter::Lace::Defaults::make_Client();
        $ds = $cl->get_DataSet_by_name($dataset_name);

        $db = $ds->get_cached_DBAdaptor;
    } else {
        $db = Hum::EnsCmdLineDB::connect();
    }

    if ($eps_file) {
        $eps_file .= '.eps' unless $eps_file =~ /\.eps$/i;
    }

    Bio::Otter::Lace::Defaults::show_help() unless $chr_name;

    if ($chr_start) {
        my $chr_start_num = parse_num($chr_start)
            or die "Bad chr_start '$chr_start'";
        $chr_start = $chr_start_num;
    }
    if ($chr_end) {
        my $chr_end_num = parse_num($chr_end)
            or die "Bad chr_end '$chr_end'";
        $chr_end = $chr_end_num;
    }
    
    my $mw = CanvasWindow::MainWindow->new;
    $mw->title($chr_name);
    my $gc = GenomeCanvas->new($mw);

    $gc->font_size( $font_size );
    
    if ($h_tiles and $v_tiles) {
        die "Can't set both htiles and vtiles";
    }
    elsif ($h_tiles) {
        $gc->horizontal_tile($h_tiles);
    }
    elsif ($v_tiles) {
        $gc->vertical_tile($v_tiles);
    }
    
    $gc->landscape($landscape);
    if ($portrait) {
        $gc->landscape(0);
    }

    my $canvas = $gc->canvas;

    {

        if ($chr_start and ! $chr_end) {
            $chr_end = get_chr_end($db, $chr_name);
        }
        elsif ($chr_end and ! $chr_start) {
            $chr_start = 1;
        }

	    my $stadp = $db->get_SliceAdaptor();

        my( $vc );
        if ($chr_start) {
            print STDERR "Fetching region of chr '$chr_name' from $chr_start to $chr_end\n";
	        $vc = $stadp->fetch_by_region(
                'chromosome',
                $chr_name,
                $chr_start,
                $chr_end,
                );
        } else {
            print STDERR "Fetching all of chromosome '$chr_name'\n";
	        $vc = $stadp->fetch_by_region('chromosome', $chr_name);
        }
        
        my $tiling = GenomeCanvas::Band::TilingPath->new($vc);
        $tiling->title("Contig\nTiling\nPath");
        if ($sanger_names) {
            $tiling->name_morpher(\&get_sanger_name);
        } else {
            $tiling->name_morpher(\&trim_name);
        }
        $tiling->tiling_direction(-1);
        $gc->add_Band($tiling);

	    # gentic markers
	    if ($markers) {
	        my $m = GenomeCanvas::Band::GeneticMap->new;
	        $m->virtual_contig($vc);
	        $m->results_file($markers);
	        $gc->add_Band($m);
	    }

        if ($show_repeats) {
            my $repeat = GenomeCanvas::DensityBand::RepeatFeature->new;
            $repeat->virtual_contig($vc);
            if ($otter) {
                warn "Attaching Client and Dataset to band";
                $repeat->Client($cl);
                $repeat->DataSet($ds);
            }
            $repeat->title('Repeats');
            $repeat->repeat_classifier(\&repeat_class);
            $gc->add_Band($repeat);
        }

        if ($show_gc_content) {
            my $plot = GenomeCanvas::Band::Plot->new;
            $plot->title('C+G Content');
            $plot->virtual_contig($vc);
            $plot->range(0.3, 0.7);
            $plot->x_major(0.2);
	        $plot->show_horizontal_scale(1);
            $plot->relative_coords($relative_coords);
            $gc->add_Band($plot);
        }
        
        if ($show_cpg) {
            my $cpg_plot = GenomeCanvas::Band::Plot->new;
            $cpg_plot->title('CpG Content');
            $cpg_plot->virtual_contig($vc);
            $cpg_plot->plot_method('cpg_profile');
            $cpg_plot->range(0.0,0.08);
            $cpg_plot->x_major(0.02);
            $cpg_plot->show_horizontal_scale(0);
            $gc->add_Band($cpg_plot);

            my $cpg = GenomeCanvas::Band::CpG->new;
            $cpg->virtual_contig($vc);
            $gc->add_Band($cpg);
        }

        if ($snps) {
            my $snp = GenomeCanvas::Band::SNP->new;
            $snp->virtual_contig($vc);
            $snp->tiling_direction(-1);
	        $snp->chr_snp_file($snps);
	        $snp->max_snp_density(30);
	        $snp->title("Total SNPs");
            $gc->add_Band($snp);
        }

	    if ($homologies) {
	        my $track = GenomeCanvas::Band::HomolFeature->new();
	        $track->feature_file( $homologies );
	        $track->type_start_end_column_indices(2, 3, 4);
	        $track->virtual_contig($vc);
	        $track->title("Similarity");
	        $gc->add_Band($track);
	    }

        if ($show_ruler) {
	        my $scale_track = GenomeCanvas::Band::Scale->new;
	        $scale_track->title("scale");
	        $scale_track->virtual_contig($vc);
	        $gc->add_Band($scale_track);
        }

	    if ($clusters) {
	        my $track = GenomeCanvas::Band::GeneCluster->new();
	        $track->cluster_file($clusters);
	        $track->span_file($span_file);
	        $track->virtual_contig($vc);
	        $track->title("Clusters");
                $track->tiling_direction(1);
	        $gc->add_Band($track);

	    }

        if ($show_genes) {
            my $gene_names = GenomeCanvas::Band::Gene->new;
            $gene_names->show_labels($gene_labels);
            $gene_names->label_type_color_list(
                [ 'Known',                  'Known_CDS',                '#4169e7' ],
                [ 'Novel CDS',              'Novel_CDS',                '#86ceff' ],
                [ 'Novel Transcript',       'Novel_Transcript',         '#ffa600' ],
                [ 'Putative',               'Putative',                 '#cccccc' ],
                #[ 'Processed Pseudogene',   'Processed_pseudogene',     '#31ce31' ],
                #[ 'Unprocessed Pseudogene', 'Unprocessed_pseudogene',   '#006400' ],

                #[ 'Known',              'Known',             '#4169e7' ],
                #[ 'Novel CDS',          'Novel_CDS',         '#86ceff' ],
                #[ 'Novel Transcript',   'Novel_Transcript',  '#ffa600' ],
                [ 'Pseudogene',         'Pseudogene',        '#31ce31' ],
                );
            $gene_names->virtual_contig($vc);
            $gene_names->span_file($span_file);
            $gene_names->tiling_direction(1);
            $gene_names->show_key($show_key);
            $gc->add_Band($gene_names);
        }

        print STDERR "Rendering canvas\n";
        $gc->residues_per_pixel($scale);
        $gc->render;
    }
    
    $gc->set_scroll_region;

    my $button_frame = $mw->Frame;
    $button_frame->pack(
        -side => 'top',
        );

    make_button($button_frame, 'Exit',
        #sub{ kill 'INT', $$ } );
        sub{ Tk::exit() } );

    make_button($button_frame, 'Zoom in', sub{ $gc->zoom(2) } );
    make_button($button_frame, 'Zoom out', sub{ $gc->zoom(-2) } );
    make_button($button_frame, 'Redraw', sub{
        $gc->render;
        $gc->fix_window_min_max_sizes
        } );
    make_button($button_frame, 'Show outline', sub{
            foreach my $obj ($canvas->find('withtag', 'all')) {
                if ($canvas->type($obj) eq 'rectangle') {
                    $canvas->itemconfigure($obj, -outline => 'red');
                }
            }
        } );

    if ($eps_file) {
        # For making figures of a fixed scale for posters
        make_button($button_frame, 'PostScript',
            sub {
                warn "Printed to file '$eps_file'\n";

                my $canvas = $gc->canvas;
                my ($x1, $y1, $x2, $y2) = $canvas->bbox('all');
                my $border = 10;
                my $width  = $x2 - $x1 + (2 * $border);
                my $height = $y2 - $y1 + (2 * $border);
                my $x = $x1 - $border;
                my $y = $y1 - $border;
                $canvas->postscript(
                    -file       => $eps_file,
                    -x          => $x,
                    -y          => $y,
                    -width      => $width,
                    -height     => $height,
                    -pageanchor => 'nw',
                    -pagex      => $border,
                    -pagey      => $border,
                    -pagewidth  => $width / 2,
                    -rotate     => 0,
                    );
                system('./patch_strings', $eps_file) == 0 or die "Error running './patch_strings' on '$eps_file'";
                if ($^O eq 'darwin') {
                    system('open', $eps_file);
                }
            });
    } else {
        make_button($button_frame, 'PostScript',
            sub {
                warn "Printed to files:\n",
                    map "  $_\n", $gc->print_postscript("$psdir/$chr_name");

                if ($big_picture) {
                    my $htile = $gc->horizontal_tile || 1;
                    my $vtile = $gc->vertical_tile   || 1;
                    my( $pwidth, $pheight );
                    if ($gc->landscape) {
                        $pwidth  = $gc->page_width  * $htile;
                        $pheight = $gc->page_height * $vtile;
                    } else {
                        $pwidth  = $gc->page_width  * $vtile;
                        $pheight = $gc->page_height * $htile;
                    }
                    if ($pwidth > $pheight) {
                        ($pwidth, $pheight) = ($pheight, $pwidth);
                    }
                    $gc->page_width ($pwidth);
                    $gc->page_height($pheight);
                    $gc->horizontal_tile(0);
                    $gc->vertical_tile(0);

                    warn "Whole diagram printed to:\n",
                        map "  $_\n", $gc->print_postscript($chr_name);
                }
            });
    }
    $gc->fix_window_min_max_sizes;
}


                            ###############
                            Tk::MainLoop();
                            ###############


sub get_chr_end {
    my( $db, $chr_name ) = @_;
    
    my $get_max = $db->prepare(q{
        SELECT length
        FROM seq_region
        WHERE name = ?
        });
    $get_max->execute($chr_name);
    my ($max) = $get_max->fetchrow;
    die "No max for chr_name = '$chr_name'" unless $max;
    return $max;
}

sub parse_num {
    my( $string ) = @_;
    
    my ($int, $mod) = $string =~ /^([\d\.]+)(.?)/ or return;
    if ($mod) {
        $mod = lc $mod;
        if ($mod eq 'k') {
            $int *= 1e3;
        }
        elsif ($mod eq 'm') {
            $int *= 1e6
        }
        else {
            return;
        }
    }
    return $int;
}

sub get_sanger_name {
    my( $name ) = @_;
    
    my $acc = trim_name($name);
    return seq_name_from_acc($acc) || "Em:$acc";
}

{
    my( $sth );
    
    sub seq_name_from_acc {
        my( $acc ) = @_;
        
        $sth ||= prepare_statement(q{
            SELECT s.sequence_name
            FROM project_acc a
              , project_dump d
              , sequence s
            WHERE a.sanger_id = d.sanger_id
              AND d.seq_id = s.seq_id
              AND d.is_current = 'Y'
              AND a.accession = ?
            });
        $sth->execute($acc);
        my ($seq_name) = $sth->fetchrow;
        return $seq_name;
    }
}

sub trim_name {
    my( $name ) = @_;
    
    my ($trim) = $name =~ /^([^\.]+)/
        or confess("Can't get trimmed name from '$name'");
    return $trim;
}


sub make_button {
    my( $parent, $label, $command ) = @_;
    
    my $button = $parent->Button(
        -text => $label,
        -command => $command,
        );
    $button->pack(
        -side => 'left',
        );
    return $button;
}

sub add_contig {
    my( $canvas, $name, $length, $phase, $gap ) = @_;

    my $fill_color = $phase == 3 ? 'black' : 'grey';

    my ($right_most) = ($canvas->bbox('contig'))[2] || 0;

    my $size = $length / 2000;  # 1 pixel = 2 kbp
    my $x1 = $right_most - 1;
    $x1 += $gap * 7.5;
    my $rec = $canvas->createRectangle(
        $x1, 0, $x1 + $size, 8,
        -fill => $fill_color,
        -outline => undef,
        -tags => ['contig', $name],
        );

    my $label = $canvas->createText(
        $x1, -1,
        -text => $name,
        #-font => '-*-helvetica-medium-r-*-*-24-*-*-*-*-*-*-*',
        -font => ['helvetica', 12],
        -anchor => 'sw',
        -tags => ['contig_label', $name],
        );

    my @bbox = $canvas->bbox($name);

    my $sp = 5;
    expand_bbox(\@bbox, $sp);
    my $bkgd = $canvas->createRectangle(
        @bbox,
        -fill    => undef,
        -outline => undef,
        -tags => ['contig_bkgd', $name],
        );
    $canvas->lower($bkgd, $rec);
    nudge_into_free_space($canvas, $name);
}

sub expand_bbox {
    my( $bbox, $pad ) = @_;
    
    $bbox->[0] -= $pad;
    $bbox->[1] -= $pad;
    $bbox->[2] += $pad;
    $bbox->[3] += $pad;
}

sub nudge_into_free_space {
    my( $canvas, $ac ) = @_;
    
    my %self = map {$_, 1} $canvas->find('withtag', $ac);
    while (grep ! $self{$_}, $canvas->find('overlapping', $canvas->bbox($ac))) {
        $canvas->move($ac, 0, -5);
    }
}

{
    my %class_hash;

    sub repeat_class {
        my( $band, $repeat_feature ) = @_;
        
	    if (not %class_hash) {
	        %class_hash = map { $_ , 1 } $band->repeat_classes;
	    }

	    my $name = $repeat_feature->repeat_consensus->name;
	    my $class = $repeat_feature->repeat_consensus->repeat_class;
	    $class =~ s{/.+}{};
	    $class = "Other" if not exists $class_hash{$class};

        return $class;
    }
    
}



__END__

=head1 NAME - genome_canvas

=head1 SYNOPSIS

  # Show the region of chromosome 6 from position 1 to 2Mbp
  # but don't show the repeats
  genome_canvas -chr chr6 -start 1 -end 2m -norepeats

=head1 AUTHOR

James Gilbert B<email> jgrg@sanger.ac.uk