#!/usr/bin/perl -Wall

use strict;
use Data::Dumper;

# LOAD CGDB->RegulonDb ID mapping
my $OPs;
open(F,'regid_cgdbid_utid.txt');
while (<F>) {
	chomp;
	my ($r,$c,$u) = split '\t';
	if (defined $OPs->{$c}) {
		$OPs->{$c} .= ','.$u; # because of stuff like
		#~ zraR    ZRAR    zraR
		#~ zraR    ZRAR    zraSR
	}
	else {
		$OPs->{$c} = $u;
	}
}
close(F);

# Load regulondb operon ids
my $ref;
open(F,'EcolA.grn.sif');
while (<F>) {
	chomp;
	my ($from, $to) = $_ =~ /^(\S+)\s\S+\s(\S+)$/;
	$ref->{$from} = 1;
	$ref->{$to} = 1;
}
close(F);

# Load profiles
my $exprs;
open(F,'EcolA.GPL3154.rma.cgdb.txt');
while (<F>) {
	chomp;
	# header
	print $_ if $. == 1; # keep header
	next if $. == 1;  # then do not process it
	my ($gene, @expr) = split '\t';
	# filter unrecognized gene profiles
	next if ! defined $OPs->{$gene};
	# store array of profiles for operon
	foreach (split ',',$OPs->{$gene}) {
		push @{$exprs->{$_}}, \@expr;
	}
}
close(F);

# Average profiles / Operon
foreach my $operon (sort keys %$exprs) {
	next if ! defined $ref->{$operon}; # filter gene profiles not in RegulonDb
	if (scalar @{$exprs->{$operon}} == 1) { # only one gene -> ouput as is
		my $profile = $exprs->{$operon}->[0];
		print "$operon\t".join("\t", @$profile);
	}
	else { # average profiles
		my @profiles = @{$exprs->{$operon}};
		my $output = $operon; # intialize line to be printed to operonId
		my $nb = scalar @{$profiles[0]}; # number of profiles columns (number of elms of first profile)
		for (my $i=0;$i<$nb;$i++) { # foreach hybridization/column
			my $sum = 0;
			for (my $j=0;$j<scalar @profiles;$j++) { # foreach gene in operon
				$sum += $profiles[$j]->[$i]; # gene j, column i
			}
			my $avg = $sum / scalar @profiles;
			$output.= "\t$avg";
		}
		print $output;
	}
}

