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

View of /FigKernelScripts/ex_fillout_ar_data.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (download) (as text) (annotate)
Mon Sep 2 02:37:59 2013 UTC (6 years, 2 months ago) by overbeek
Branch: MAIN
CVS Tags: rast_rel_2014_0729, rast_rel_2014_0912, HEAD
added steps

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

my $usage = "usage: ex_fillout_ar_data -d DataDir\n";

my $dataD;
my $rc  = GetOptions('d=s' => \$dataD);

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

my %peg_to_alias;
my %alias_to_peg;

foreach $_ (`cat $dataD/aliases`)
{
    if ($_ =~ /^(fig\S+)\t\S*\t(\S+)$/)
    {
	$peg_to_alias{$1} = $2;
	$alias_to_peg{$2} = $1;
    }
}

my %ar_counts;
foreach $_ (`cat $dataD/atomic.regulon.on.off.calls`)
{
    if ($_ =~ /^\S+\t(\S+)\t(\S+)$/)
    {
	my $ar = $1;
	my $v = $2;
	if (! defined($ar_counts{$ar})) { $ar_counts{$ar} = [0,0,0] }
	$ar_counts{$ar}->[$v+1] += 1;
    }
}

my %peg_counts;
foreach $_ (`cat $dataD/peg.on.off.calls`)
{
    if ($_ =~ /^\S+\t(\S+)\t(\S+)$/)
    {
	my $peg = $1;
	my $v = $2;
	if (! defined($peg_counts{$peg})) { $peg_counts{$peg} = [0,0,0] }
	$peg_counts{$peg}->[$v+1] += 1;
    }
}

my %pcc = map { ($_ =~ /^(\S+\t\S+)\t(\S+)/) ? ($1 => $2) : () } `cat $dataD/locus.pccs`;

open(EXP,">$dataD/atomic.regulons.expanded") || die "could not open $dataD/atomic.regulons.expanded";
open(AR,"<$dataD/atomic.regulons") || die "could not open $dataD/atomic.regulons";
my $last = <AR>;
while ($last && ($last =~ /^(\S+)/))
{
    my $ar = $1;
    my @set;
    while ($last && ($last =~ /^(\S+)\t(\S+)/) && ($1 eq $ar))
    {
	push(@set,$2);
	$last = <AR>;
    }
    &process($dataD,$ar,\@set,\%peg_to_alias,\%alias_to_peg,\%ar_counts,\%peg_counts,\%pcc,\*EXP);
}
close(EXP);

sub process {
    my($dataD,$ar,$set,$peg_to_alias,$alias_to_peg,$ar_counts,$peg_counts,$pccs,$exp_fh) = @_;

    my $tot = 0;
    my @locus_ids = grep { $_ =~ /^\S/ } map { $peg_to_alias->{$_} } @$set;
    my $ar_on_off = $ar_counts->{$ar};
    if (! $ar_on_off) { $ar_on_off = [0,0,0] }
    foreach my $peg (sort { &SeedUtils::by_fig_id($a,$b) } @$set)
    {
	my $alias = $peg_to_alias->{$peg};
	if (! $alias) { $alias = '' }
	my $peg_on_off = $peg_counts->{$peg};
	if (! $peg_on_off) { $peg_on_off = [0,0,0] }
	my @others = grep { $_ && ($_ ne $alias) } map { $peg_to_alias->{$_} } @$set;
	my $tot = 0;
	foreach my $other (@others)
	{
	    if (($_ = $pccs->{"$alias\t$other"}) || ($_  = $pccs->{"$other\t$alias"}))
	    {
		$tot += $_;
	    }
	}
	my $pcc = (@others > 0) ? sprintf("%0.3f",($tot/@others)) : 0 ;
	print $exp_fh join("\t",($ar,$peg,$alias,$pcc,@$ar_on_off,@$peg_on_off)),"\n";
    }
}

			     
			     

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3