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

View of /FigKernelScripts/ex_make_on_off_calls.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (download) (as text) (annotate)
Sun Jul 28 15:56:13 2013 UTC (6 years, 4 months ago) by overbeek
Branch: MAIN
Changes since 1.1: +1 -1 lines
fix on arg order - ids

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

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

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

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

my $initial_calls = &make_initial_calls($dataD);
open(CALLS,">$dataD/peg.on.off.calls")
    || die " could not open $dataD/initial.peg.on.off.calls";
foreach my $exp (sort keys(%$initial_calls))
{
    my $locusH = $initial_calls->{$exp};
    foreach my $locus (keys(%$locusH))
    {
	my $v = $locusH->{$locus};
	print CALLS join("\t",($exp,$locus,$v)),"\n";
    }
}
close(CALLS);

sub make_initial_calls {
    my($dataD) = @_;

    my %roles_always_on = map { chop; ($_ => 1) } `cat $dataD/roles.always.on`;
    my %pegs_always_on = map { chop; 
			       my(undef,$locus,undef,$func) = split(/\t/,$_); 
			       $roles_always_on{$func} ? ($locus => 1) : () }
                         `cat $dataD/aliases.with.function`;
    my %exp_peg_val;
    open(EPV,"<$dataD/gene.locus.exp.val.tuples") 
	|| die "could not open $dataD/gene.locus.exp.val.tuples";
    while (defined($_ = <EPV>))
    {
	chop;
	my($gene,$locus,$exp,$v) = split(/\t/,$_);
	$exp_peg_val{$exp}->{$locus} = $v;
    }
    close(EPV);
    my %on_for_exp;
    my %off_for_exp;
    my %diff_by_exp;
    foreach my $exp (keys(%exp_peg_val))
    {
	my $locusH = $exp_peg_val{$exp};
	my @on_values = sort { $a <=> $b } 
	                map { $pegs_always_on{$_} ? $locusH->{$_} : () } keys(%$locusH);
	my $I_10 = int(@on_values * 0.1);
	my $on_threshold = $on_values[$I_10];
	$on_for_exp{$exp} = $on_threshold;
	my @not_on_values = sort { $a <=> $b } 
            	            map { ((! $pegs_always_on{$_} ) && ($locusH->{$_} < $on_threshold)) ? $locusH->{$_} : () } 
	                    keys(%$locusH);
	my $I_80  = int(@not_on_values * 0.8);
	my $off_threshold = $not_on_values[$I_80];
	$off_for_exp{$exp} = $off_threshold;
	$diff_by_exp{$exp} = $on_threshold - $off_threshold;
    }
    my @diffs = sort { $a <=> $b } map { $diff_by_exp{$_} } keys(%diff_by_exp);
    my $I_25  = int(@diffs * 0.25);
    my $diff_threshold = $diffs[$I_25];
    foreach my $exp (keys(%off_for_exp))
    {
	if (($on_for_exp{$exp} - $off_for_exp{$exp}) < $diff_threshold)
	{
	    $off_for_exp{$exp} = $on_for_exp{$exp} - $diff_threshold;
	}
    }
    my %peg_on_off_by_exp;
    foreach my $exp (keys(%exp_peg_val))
    {
	my $on_threshold = $on_for_exp{$exp};
	my $off_threshold = $off_for_exp{$exp};

	my $locusH = $exp_peg_val{$exp};
	foreach my $peg (keys(%$locusH))
	{
	    my $v = $locusH->{$peg};
	    if ($v >= $on_threshold)
	    {
		$peg_on_off_by_exp{$exp}->{$peg} = 1;
	    }
	    elsif ($v <= $off_threshold)
	    {
		$peg_on_off_by_exp{$exp}->{$peg} = -1;
	    }
	    else
	    {
		$peg_on_off_by_exp{$exp}->{$peg} = 0;
	    }
	}
    }
    return \%peg_on_off_by_exp;
}

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3