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

View of /FigKernelScripts/ex_make_initial_atomic_regulons.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (download) (as text) (annotate)
Fri Nov 22 22:08:59 2013 UTC (5 years, 11 months ago) by overbeek
Branch: MAIN
CVS Tags: rast_rel_2014_0729, rast_rel_2014_0912, HEAD
Changes since 1.1: +2 -2 lines
two embarrassing bugs in same_profile

use Data::Dumper;
use File::Basename;
use Getopt::Long;
use SeedEnv;
use strict;

my $usage = "usage: ex_make_initial_atomic_regulons -d DataDir\n";
my($i,$j);
my $dataD;
my $rc  = GetOptions('d=s' => \$dataD);

if ((! $rc) || (! -d $dataD))
{ 
    print STDERR $usage; exit ;
}

my %to_peg = map { ($_ =~ /(fig\|\d+\.\d+\.peg\.\d+).*\t(\S[^\s]+\S)$/) ? ($2 => $1) : () }
             `cat $dataD/aliases`;

open(CONN,"| cluster_objects | tabs2rel > $dataD/initial.atomic.regulons") 
    || die "could not open $dataD/initial.atomic.regulons";
foreach my $line (`cat $dataD/coreg.based.on*`)
{
    if ($line =~ /^(\S+)/)
    {
	my @pegs = split(/,/,$1);
	for ($i=0; ($i < $#pegs); $i++)
	{
	    for ($j=$i+1; ($j < @pegs); $j++)
	    {
		print CONN join("\t",($pegs[$i],$pegs[$j])),"\n";
	    }
	}
    }
}

open(CALLS,"<$dataD/peg.on.off.calls") 
    || die "could not open $dataD/peg.on.off.calls";
my %pegs;
my %exp_peg_call;
my $line;
while (defined($line = <CALLS>))
{
    if (($line =~ /^(\S+)\t(\S+)\t(\S+)$/) && ($_ = $to_peg{$2}))
    {
	$exp_peg_call{$1}->{$_} = $3;
	$pegs{$_} = 1;
    }
}
close(CALLS);

my @exp = sort keys(%exp_peg_call);
my @pegs = keys(%pegs);
my $pvecs = &to_pvec(\@exp,\@pegs,\%exp_peg_call);
for ($i=0; ($i < $#pegs); $i++)
{
    my $pvecsI = $pvecs->{$pegs[$i]};
    print STDERR ".";
    for ($j=$i+1; ($j < @pegs); $j++)
    {
	my $pvecsJ = $pvecs->{$pegs[$j]};
	if (&same_profile($pvecsI,$pvecsJ))
	{
	    print CONN join("\t",($pegs[$i],$pegs[$j])),"\n";
	}
    }
}
close(CONN);

sub to_pvec {
    my($exps,$pegs,$exp_peg_call) = @_;

    my $pvecs = {};
    foreach my $peg (@$pegs)
    {
	foreach my $exp (@$exps)
	{
	    push(@{$pvecs->{$peg}},$exp_peg_call->{$exp}->{$peg});
	}
    }
    return $pvecs;
}

sub same_profile {
    my($pvecsI,$pvecsJ) = @_;

    my $i;
    my $solid1_on = 0;
    my $solid1_off = 0;
    my $solid2_off = 0;
    my $solid2_on = 0;

    my $ok = 1;
    for ($i=0; ($i < @$pvecsI) && $ok; $i++) 
    {	
	my $v1 = $pvecsI->[$i];
	my $v2 = $pvecsJ->[$i];
	$ok = &ok($v1,$v2);
	if    ($v1 == 1)    { $solid1_on++ }
	elsif ($v1 == -1)   { $solid1_off++ }
	if    ($v2 ==  1)   { $solid2_on++ }
	elsif ($v2 == -1)   { $solid2_off++ }
    }
    return (($i == @$pvecsI) && 
	    ($solid1_on  > (@$pvecsI * 0.2)) &&
	    ($solid1_off > (@$pvecsI * 0.2)) &&
	    ($solid1_on  > (@$pvecsI * 0.2)) &&
	    ($solid2_off > (@$pvecsJ * 0.2)));
}

sub ok {
    my($x,$y) = @_;

    return (($x == $y) || ($x == 0) || ($y == 0));
}

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3