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

View of /FigKernelScripts/pg_intersection.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (download) (as text) (annotate)
Mon Jun 3 20:20:04 2013 UTC (6 years, 5 months ago) by overbeek
Branch: MAIN
CVS Tags: rast_rel_2014_0729, rast_rel_2014_0912, HEAD
add Jose stuff

use strict;
use Data::Dumper;
use Getopt::Long;


#########
my $usage = "usage: pg_intersection -d Data\n";
my $dataD;

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

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

###################################################################
opendir(SETS,"$dataD/ReactionSets")  || die "No subsets defined yet";
my @sets = grep { $_ !~ /^\./ } readdir(SETS);
closedir(SETS);
my $dir = "$dataD/ReactionSetIntersections";
mkdir($dir,0777);

foreach my $set1 (@sets)
{
    foreach my $set2 (grep { $_ lt $set1 } @sets)
    {
	my @in1 = map { ($_ =~ /(rxn\d+)/) ? $1 : () } `cat $dataD/ReactionSets/$set1/core_reactions.txt`;
	my @in2 = map { ($_ =~ /(rxn\d+)/) ? $1 : () } `cat $dataD/ReactionSets/$set2/core_reactions.txt`;
	my %in2 = map { $_ => 1 } @in2;
	my @inter = grep { $in2{$_} } @in1;
	my %inter = map { $_ => 1 } @inter;
	my @just1 = grep { ! $inter{$_} } @in1;
	my @just2 = grep { ! $inter{$_} } @in2;
	open(I,">$dir/$set1.and.$set2") || die "could not open $dir/$set1.and.$set2";
	foreach my $r (@inter) { print I "$r\n"; }
	close(I);
	open(D1,">$dir/$set1.not.$set2") || die "could not open $dir/$set1.not.$set2";
	foreach my $r (@just1) { print D1 "$r\n"; }
	close(D1);
	open(D2,">$dir/$set2.not.$set1") || die "could not open $dir/$set2.not.$set1";
	foreach my $r (@just2) { print D2 "$r\n"; }
	close(D2);
    }
    my %reactions;
    my @tuples = map { ($_ =~ /^(\S+)\s.*(rxn\d+)/) ? [$2,$1] : () } `cat $dataD/ReactionSets/$set1/non_core_reactions.txt`;
    foreach my $tuple (@tuples)
    {
	my($reaction,$genome) = @$tuple;
	$reactions{$reaction}->{$genome} = 1;
    }
    open(F,">$dir/$set1.non-core.genome-counts") || die "could not open $dir/$set1.non-core.genome-counts";
    foreach my $r (sort keys(%reactions))
    {
	my $gH = $reactions{$r};
	my @gs = keys(%$gH);
	print F join("\t",($r,scalar @gs,join(",",@gs))),"\n";
    }
    close(F);
}


MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3