[Bio] / FigKernelScripts / compare_two_sets_of_calls.pl Repository:
ViewVC logotype

View of /FigKernelScripts/compare_two_sets_of_calls.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.3 - (download) (as text) (annotate)
Sat Apr 11 22:16:59 2009 UTC (10 years, 10 months ago) by overbeek
Branch: MAIN
CVS Tags: mgrast_dev_08112011, rast_rel_2009_05_18, mgrast_dev_08022011, rast_rel_2014_0912, myrast_rel40, mgrast_dev_05262011, mgrast_dev_04082011, rast_rel_2010_0928, mgrast_version_3_2, mgrast_dev_12152011, mgrast_dev_06072011, rast_rel_2009_0925, rast_rel_2010_0526, rast_rel_2014_0729, mgrast_dev_02212011, rast_rel_2010_1206, mgrast_release_3_0, mgrast_dev_03252011, rast_rel_2010_0118, rast_rel_2011_0119, mgrast_release_3_0_4, mgrast_release_3_0_2, mgrast_release_3_0_3, mgrast_release_3_0_1, mgrast_dev_03312011, mgrast_release_3_1_2, mgrast_release_3_1_1, mgrast_release_3_1_0, mgrast_dev_04132011, mgrast_dev_04012011, rast_rel_2009_07_09, rast_rel_2010_0827, myrast_33, rast_rel_2011_0928, mgrast_dev_04052011, mgrast_dev_02222011, mgrast_dev_10262011, HEAD
Changes since 1.2: +37 -13 lines
stuff for gene calling

use TBLstuff;
use FIG;
my $fig = new FIG;
use Data::Dumper;
use Carp;
use strict;
use compare_coding;

my $usage = "usage: compare_two_sets_of_calls [-h HTMLdir] [-g Genome] [-p PegTbl -r RNAtbl]";

my($html_dir,$genome,$peg_tbl,$rna_tbl);

while ((@ARGV > 0) && ($ARGV[0] =~ s/^-//))
{
    $_ = shift;
    if     (s/^h//)         { $html_dir   = $_ || shift }
    elsif  (s/^g//)         { $genome     = $_ || shift }
    elsif  (s/^p//)         { $peg_tbl    = $_ || shift }
    elsif  (s/^r//)         { $rna_tbl    = $_ || shift }
    else
    {
	die "Invalid Flag '$_':\n$usage";
    }
}

my @features2;

if ($genome)
{
    @features2 = &TBLstuff::current_features_of($fig,$genome);
}
else
{
    ($peg_tbl && $rna_tbl) || die "you need to give either a genome id or a pair of peg/rna tbls\n$usage";
    @features2 = &TBLstuff::flat_tbl_to_features($peg_tbl,$rna_tbl);
}

my @features1 = map { chomp; 
                      my($type,$loc,$sc,%extras) = split(/\t/,$_); 
		      $loc =~ s/^[^:]+://;
		      $extras{score} = $sc;
		      [$type,$loc,\%extras] 
		    }
                <STDIN>;

my @merged = &TBLstuff::compare_sets_of_features($fig,\@features1,\@features2);
#
# [relation,contig,midpt,call-set-1-entry,call-set-2-entry]
#
#     relation: { same, diff-start, '' means unique }
#

@merged = sort { ($a->[1] cmp $b->[1]) or ($a->[2] <=> $b->[2]) } @merged;  # sort by position
@merged = map { [$_->[0],$_->[3]->[1],$_->[4]->[1]] } @merged;              # [relation,call-set-1,call-set2]
                                                                            # call-set-entry is
                                                                            #   [feature-type,loc,\%extras] 
                                                                            # useful extras:  score => $score, query => $queryID
my $html = &to_html(\@merged,$html_dir);
print $html;

sub to_html {
    my($merged,$html_dir) = @_;

    my @html = ();
    push @html, <<"End_of_head";
<HTML>
<HEAD>
<TITLE>Comparision of new to old gene calls</TITLE>
</HEAD>
<BODY>
<TABLE>
End_of_head
;
    foreach ( @$merged )
    {
	my ( $rel, $item1, $item2 ) = @$_;
	my $type   = $item1 ? $item1->[0] : $item2->[0];

	my $extra1 = $item1 ? $item1->[2] : {};
	my $loc1   = $item1 ? $item1->[1] : '';
	my $score  = $item1 ? $extra1->{ score } : '';

	my $extra2 = $item2 ? $item2->[2] : {};
	my $loc2   = $item2 ? $item2->[1] : '';
	my $id2    = $item2 ? $extra2->{ id } : '';
	my $query  = $item1 ? query_link( $extra1->{ query }, $extra1->{ orig_id },
					  $html_dir) 
	                    : '';
	if (($rel =~ /^diff-/) && $loc1 && $loc2)
	{
	    my($b1,$b2,$e1,$e2,@locs1,@locs2);
	    @locs1 = split(/,/,$loc1);
	    @locs2 = split(/,/,$loc2);
	    
	    if (($locs1[0]  =~ /^\S+_(\d+)_\d+$/) && ($b1 = $1) &&
		($locs1[-1] =~ /^\S+_\d+_(\d+)$/) && ($e1 = $1) &&
		($locs2[0]  =~ /^\S+_(\d+)_\d+$/) && ($b2 = $1) &&
		($locs2[-1] =~ /^\S+_\d+_(\d+)$/) && ($e2 = $1))
	    {
		my $db  = ($b1 > $e1) ? $b1-$b2 : $b2-$b1;
		my $de  = ($b1 < $e1) ? $e1-$e2 : $e2-$e1;
		$rel = sprintf("%+d:%+d",$db,$de);
	    }

	}
        push @html, "  <TR>",
                    "    <TD>$rel</TD>",
	            "    <TD>$type</TD>",
                    "    <TD>$query</TD>",
                    "    <TD align=right>$score&nbsp;&nbsp;</TD>",
                    "    <TD>$loc1</TD>",
                    "    <TD>$loc2</TD>",
                    "    <TD>$id2</TD>",
	            "  </TD>";
    }
    push @html, <<"End_of_tail";
</TABLE>
</BODY>
</HTML>
End_of_tail

    return join( "\n", @html );
}


sub query_link
{
    my ( $query,$orig_id,$html_dir) = @_;
    return '' if ! $query;
    return $query if ! $html_dir;

    return "<A HRef=\"$html_dir/$query.html\">$query</A>";
}
    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3