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

View of /FigKernelScripts/update_FIGfams.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.12 - (download) (as text) (annotate)
Thu Oct 8 18:51:57 2009 UTC (10 years, 1 month ago) by arodri7
Branch: MAIN
CVS Tags: mgrast_dev_08112011, mgrast_dev_08022011, rast_rel_2014_0912, myrast_rel40, mgrast_dev_05262011, mgrast_dev_04082011, rast_rel_2010_0928, mgrast_version_3_2, mgrast_dev_12152011, mgrast_dev_06072011, rast_rel_2010_0526, rast_rel_2014_0729, mgrast_dev_02212011, rast_rel_2010_1206, mgrast_release_3_0, mgrast_dev_03252011, rast_rel_2010_0118, rast_rel_2011_0119, mgrast_release_3_0_4, mgrast_release_3_0_2, mgrast_release_3_0_3, mgrast_release_3_0_1, mgrast_dev_03312011, mgrast_release_3_1_2, mgrast_release_3_1_1, mgrast_release_3_1_0, mgrast_dev_04132011, mgrast_dev_04012011, rast_rel_2010_0827, myrast_33, rast_rel_2011_0928, mgrast_dev_04052011, mgrast_dev_02222011, mgrast_dev_10262011, HEAD
Changes since 1.11: +233 -34 lines
commit changes

########################################################################
# 
# FIGfams are designed to automatically update in response to changes in
# the collection of subsystems.  There should be no need to do anything
# but work on decision procedures.  The basic operations implemented in
# this code are as follows:
# 
# 1. When a FIGfam is found to contain entries from subsystems that have
#    different functions (not counting comments), the FIGfam must be
#    split.  
# 
# 2. When distinct subsystems exist, each containing PEGs with identical
#    functions and they are globally similar, the FIGfams
#    must be merged.
# 
# 3. When new PEGs are found to be in subsystems, but not in the
#    appropriate FIGfam, they need to be added.
# 
# 4. When new PEGs are found to be in subsystems, but there is no
#    appropriate FIGfam, a new FIGfam must be created
# 
# This last case is a bit unusual, since it requires the generation of a
# new "partition".
# 
# In any event, once the appropriate changes have been made, the
# functions need to be recomputed.  This all is done on a weekly basis,
# so there is no meaningful notion of "release" other than the first and
# the current state.
# 
###############################

my $usage = "usage: update_FIGfams [-f File] [-t tmpdir] NewFamilies2c NewFunctions";
my($new2c,$new_func,$ss_file);

my $tmpdir = "$FIG_Config::temp/update.FIGfams.$$";
while ( $ARGV[0] =~ /^-/ )
{
    $_ = shift @ARGV;
    if    ($_ =~ s/^-f//) { $ss_file       = ($_ || shift @ARGV) }
    elsif ($_ =~ s/^-t//) { $tmpdir        = ($_ || shift @ARGV) }
    else                  { print STDERR  "Bad flag: '$_'\n$usage"; exit 1 }
}

(
 ($new2c    = shift @ARGV) &&
 ($new_func = shift @ARGV)
)
    || die $usage;

my $oldFF = shift @ARGV;

use FIG;
my $fig = new FIG;

use strict;
mkdir($tmpdir,0777) || die "could not make $tmpdir";

my($peg_to_fam,$members_in_fam,$next_famP);
#my($peg_to_fam,$members_in_fam,$next_famP) = &existing_families($oldFF,$fig);
my $peg_to_ss_fam;
#my $peg_to_ss_fam                          = &get_current_subsys_based_families($fig,$tmpdir,$ss_file);

my $n;
#my $n = keys(%$peg_to_fam);      print STDERR "$n pegs in existing families\n";
#$n    = keys(%$members_in_fam);  print STDERR "$n existing families\n";
#$n    = keys(%$peg_to_ss_fam);   print STDERR "$n pegs in new ss_fams\n";

# after 8 hours of $fig being open in a process the hash will expire, and if it does it generates an error and program dies ... regenerate $fig hash
undef $fig;
my $fig = new FIG;
my $state = [$fig,$tmpdir,$peg_to_fam,$members_in_fam,$peg_to_ss_fam,$next_famP];
print STDERR "created initial state\n";

#&delete_if_necessary($state);
#&split_if_necessary($state);
#&add_new_if_necessary($state);
#&merge_if_possible($state);

#&dump_families_2c($fig,$peg_to_fam,$new2c);
#&dump_md5_figfams($fig,$new2c);
&dump_family_functions2($fig,$new2c,$new_func);

sub dump_families_2c {
    my($fig,$peg_to_fam,$new2c) = @_;

    open(FAM2C,"| sort > $new2c") || die "could not open $new2c";
    my($peg,$fam);
    my $prev_fam ="";
    my $fam_pegs=[];
    while (($peg,$fam) = each(%$peg_to_fam))
    {
	#print FAM2C "$fam\t$peg\n";
	if ($prev_fam && $fam ne $prev_fam)
	{
	    if (scalar @$fam_pegs > 1)
	    {
		print FAM2C join ("", @$fam_pegs);
		$fam_pegs=[];
	    }
	}
	push @$fam_pegs, "$fam\t$peg\n";
	$prev_fam = $fam;

	# change the role of the peg and subsystems if the figfam id in the function does not match the new figfam id
#	my $role = $fig->function_of($peg);
#	if ($role =~ /^(FIG\d{6})(:.*)/)
#	{
#	    my $ff = $1;
#	    my $real_func = $2;
#	    if ($ff ne $fam)
#	    {
#		my $newname = $fam . "$real_func";
#		my $seeduser = "FIGfam_update";    ###### what user should I use here? should I just put FIGfam_update?
#		#$fig->change_funcrole($role,$newname,$seeduser);
#		$fig->assign_function($peg,$seeduser,$newname);
#		print STDERR "FUNCLOG\t$peg\t$seeduser\t$newname\n";
#		#print STDERR "changing $ff: funcrole $role to $newname\n";
#	    }
#	}
#	elsif ($role =~ /^(FIG\d{6})$/)
#	{
#	    my $newname = "hypothetical protein";
#	    my $seeduser = "FIGfam_update";    ###### what user should I use here? should I just put FIGfam_update?
#	    $fig->change_funcrole($role,$newname,$seeduser);
#	    print STDERR "FUNCLOG\t$role\t$seeduser\t$newname\n";
#	    #print STDERR "changing funcrole $role to $newname\n";
#	}
    }

    if  (scalar @$fam_pegs > 1)
    {   
	print FAM2C join ("", @$fam_pegs);
	$fam_pegs=[];
    }

    close(FAM2C);
}

sub dump_md5_figfams {
    my ($fig,$new2c) = @_;

    my @paths = split(/\//, $new2c);
    pop @paths;
    my $dir = join("/", @paths);

    if (-e "$dir/md5.figfams") { system "rm $dir/md5.figfams*" }
    open(OUT,">$dir/md5.figfams") || die "could not open $dir/md5.figfams";

    my $md5Hash = {};
    open(TMP,"<$new2c")
        || die "could not open $new2c";
    while (defined($_ = <TMP>))
    {
        if ($_ =~ /^(FIG\d{6})\t(\S+)/)
        {
            my $fam  = $1;
            my $peg  = $2;
            my $md5 = $fig->md5_of_peg($peg);
            push (@{$md5Hash->{$md5}}, $fam);
        }
    }
    close(TMP);

    foreach my $md5 (sort keys %{$md5Hash}){
        my %saw;
        @saw{@{$md5Hash->{$md5}}} = ();
        my @array = sort keys %saw;  # remove sort if undesired

        print OUT join(",", @array) . "\t$md5\n";
    }
    close OUT;

}

sub dump_family_functions2 {
    my ($fig,$new2c,$new_func) = @_;

    my $family_names ={};
    open(FAMFUNC,">$new_func") || die "could not open $new_func";
    open(TMP,"<$new2c") || die "could not open $new2c";
    my $prev_fam;
    my $fam_pegs=[];
    my $fam_members={};
    my ($fam, $peg);

    while (defined($_ = <TMP>))
    {
	if ($_ =~ /^(FIG\d{6})\t(\S+)/)
        {
            $fam  = $1;
            $peg  = $2;

	    if ($prev_fam && $fam ne $prev_fam)
	    {
		$fam_members->{$fam} = $fam_pegs;
		my $func = &function_of_family($fig,$fam,$fam_members);
		push (@{$family_names->{$func}}, $prev_fam);

		# clean out the fam_pegs array
		$fam_pegs=[];
            }
	    push @$fam_pegs, $peg;
	    $prev_fam = $fam;
	}
    }
    close TMP;

    if (scalar @$fam_pegs > 1)
    {
	$fam_members->{$prev_fam} = $fam_pegs;
	my $func = &function_of_family($fig,$prev_fam,$fam_members);
	push (@{$family_names->{$func}}, $prev_fam);
    }
    
    # count the number of times the function occurs, so we can disambiguate FIGfams of same function name
    my $new_fam_funcs={};
    foreach my $func (sort {$family_names->{$b} <=> $family_names->{$a}} keys %$family_names)
    {
        if (scalar @{$family_names->{$func}} > 1 && $func !~ /not subsystem-based/ && $func !~ /^FIG\d\d\d\d\d\d/)
        {
            foreach my $fam (@{$family_names->{$func}})
            {
                $new_fam_funcs->{$fam} = "$fam: $func";
            }
        }
        else
        {
            foreach my $fam (@{$family_names->{$func}})
            {
                $new_fam_funcs->{$fam} = $func;
            }
        }
    }

    foreach my $fam (sort {$a cmp $b} keys %$new_fam_funcs)
    {
        print FAMFUNC "$fam\t".$new_fam_funcs->{$fam}."\n";
        #print FAMFUNC "$fam\t$func\n";
        #print FAMFUNC "$fam\t$func\t$qty\n";
    }
    close(FAMFUNC);
    
}

sub dump_family_functions {
    my($fig,$members_in_fam,$new_func) = @_;

    my $family_names={};
    open(FAMFUNC,">$new_func") || die "could not open $new_func";
    foreach my $fam (sort keys(%$members_in_fam))
    {
	my $func = &function_of_family($fig,$fam,$members_in_fam);
	push (@{$family_names->{$func}}, $fam);
    }

    # count the number of times the function occurs, so we can disambiguate FIGfams of same function name
    my $new_fam_funcs={};
    foreach my $func (sort {$family_names->{$b} <=> $family_names->{$a}} keys %$family_names)
    {
	if (scalar @{$family_names->{$func}} > 1 && $func !~ /not subsystem-based/ && $func !~ /^FIG\d\d\d\d\d\d/)
	{
	    foreach my $fam (@{$family_names->{$func}})
	    {
		$new_fam_funcs->{$fam} = "$fam: $func";
	    }
	}
	else
	{
	    foreach my $fam (@{$family_names->{$func}})
	    {
		$new_fam_funcs->{$fam} = $func;
	    }
	}
    }

    foreach my $fam (sort {$a cmp $b} keys %$new_fam_funcs)
    {
	print FAMFUNC "$fam\t".$new_fam_funcs->{$fam}."\n";
	#print FAMFUNC "$fam\t$func\n";
	#print FAMFUNC "$fam\t$func\t$qty\n";
    }
    close(FAMFUNC);
}

sub function_of_family {
    my($fig,$fam,$members_in_fam) = @_;
    my($i,$func,@subs,%poss);

    my $set = $members_in_fam->{$fam};
    $set || confess "bad fam: $fam";

    my $ffunc = "";
    for ($i=0; (! $ffunc && ($i < @$set)); $i++)
    {
	my $peg = $set->[$i];
	$func = &function_of_peg($fig,$peg);
	next if (&ignore_function($func));
	#next if ($fig->hypo($func));
	@subs = grep { $fig->usable_subsystem($_) } $fig->peg_to_subsystems($peg);
	if (@subs > 0)
	{
	    $ffunc = $func;
	    
	}
	$poss{$func}++;
    }

    if (! $ffunc)
    {
	my @tmp = sort { $poss{$b} <=> $poss{$a} } grep { $_ } keys(%poss);
	if (@tmp > 0)
	{
	    $ffunc = "$fam (not subsystem-based): " . $tmp[0];
	}
	else
	{
	    $ffunc = "$fam (not subsystem-based): hypothetical protein";;
	}
    }

    # strip off any comments
    if (($ffunc !~ /\#\#/) && ($ffunc =~ /\#/))
    {
        my @parts = split (/\#/, $ffunc);
	$ffunc = $parts[0];
    }

    #
    ## add disambiguating code here
    #
    if ($ffunc !~ /not subsystem-based/ && $ffunc !~ /^FIG\d\d\d\d\d\d/)
    {
	# do not changeanything for now
    }

    return $ffunc;
}

sub ignore_function {
    my $x = (@_ == 1) ? $_[0] : $_[1];

    if (! $x)                             { return 1 }
    if ( ($x =~ /predicted by Psort/) )
    {
        return 1;
    }
    return 0;

}

sub function_of_peg {
    my($fig,$peg) = @_;

    my $func = $fig->function_of($peg,,1);
    $func =~ s/^FIG\d{6} \(not subsystem-based\): //;
    return $func;
}

sub get_current_subsys_based_families {
    my($fig,$tmpdir,$ss_file) = @_;

    my $peg_to_ss_fam = {};
    if ($ss_file)
    {
	system "cp $ss_file $tmpdir/subsys.based.families";
    }
    else
    {
	#&FIG::run("make_subsys_based_families2 > $tmpdir/subsys.based.families");
	&FIG::run("make_subsys_based_families > $tmpdir/subsys.based.families");
    }
    
    # after 8 hours of $fig being open in a process the hash will expire, and if it does it generates an error and program dies ... regenerate $fig hash
    undef $fig;
    $fig = new FIG;

    open(SSFAMS,"cut -f1,2 $tmpdir/subsys.based.families |")
	|| die "could not open $tmpdir/subsys.based.families";

    my $last = <SSFAMS>;
    my $ss_count=0;
    while ($last)
    {
	if ($last =~ /^(\S+)\t(fig\|\d+\.\d+\.peg\.\d+)/)
	{
	    my $fam = $1;
	    my $set = [];
	    print STDERR "SS processing $fam\n" if ($ss_count%100 == 0);
	    $ss_count++;
	    while ($last && ($last =~ /^(\S+)\t(fig\|\d+\.\d+\.peg\.\d+)/) && ($1 eq $fam))
	    {
		push(@$set,$2);
		$last = <SSFAMS>;
	    }
	    &process_ss_fam($fig,$set,$peg_to_ss_fam);
	}
	else
	{
	    print STDERR "BAD SUBSYS-BASED: $last";
	    $last = <SSFAMS>;
	}
    }
    close(SSFAMS);
    return $peg_to_ss_fam;
}

sub existing_families {
    my ($oldFF,$fig) = @_;
    
    my $peg_to_fam = {};
    my $members_in_fam = {};
    my $max_fam = "FIG000000";

    my $ffdata = &FIG::get_figfams_data($oldFF);
    open(FAM2C,"cut -f1,2 $ffdata/families.2c |")
	|| die "could not access $ffdata/families.2c";
    my $last = <FAM2C>;
    while ($last && ($last =~ /^(FIG\d{6})\t(fig\|\d+\.\d+\.peg\.\d+)/))
    {
	my $fam = $1;
	my $set = [];
	while ($last && ($last =~ /^(FIG\d{6})\t(fig\|\d+\.\d+\.peg\.\d+)/) && ($1 eq $fam))
	{
	    my $peg = $2;
	    if (! $fig->is_deleted_fid($peg))
	    {
		push(@$set,$peg);
		$peg_to_fam->{$peg} = $fam;
	    }
	    $last = <FAM2C>;
	}
	if (scalar @$set > 1)
	{
	    $members_in_fam->{$fam} = $set;
	}
	else
	{
	    foreach my $peg (@$set)
	    {
		delete $peg_to_fam->{$peg};
	    }
	}
    
	if ($fam gt $max_fam) { $max_fam = $fam }
    }
    $max_fam++;
    return ($peg_to_fam,$members_in_fam,\$max_fam);
}
	
sub process_ss_fam {
    my($fig,$set,$peg_to_ss_fam) = @_;

    my @clusters = &break_set_into_close($fig,$set);
    foreach my $cluster (@clusters)
    {
	foreach my $peg (@$cluster)
	{
	    $peg_to_ss_fam->{$peg} = $cluster;
	}
    }
}

sub break_set_into_close {
    my($fig,$set) = @_;
    my($peg,$sim,%hits);

    my @clusters = ();
    my %in = map { $_ => 1 } @$set;
    while ($peg = shift @$set)
    {
	my $func1 = $fig->function_of($peg,,1);
	my $ln1 = &length_of($fig,$peg);
#	my @sims = grep { $in{$_->id2} } $fig->sims($peg,1000,1.0e-20,"fig");
	my @sims = grep { $in{$_->id2} } $fig->sims($peg,1000,1.0e-5,"fig");
	my %hits;
	foreach $sim (@sims)
	{
	    my $id2  = $sim->id2;
	    if ((($sim->e1 + 1 - $sim->b1)/$ln1) > 0.7)
	    {
		my $ln2 = &length_of($fig,$id2);
		if ((abs($ln2 - $ln1) / $ln1) <= 0.3)
		{
		    my $func2 = $fig->function_of($id2,,1);
		    if ($func1 eq $func2)
		    {
			$hits{$id2} = 1;
		    }
		}
	    }
	}
	my $clust = [$peg,keys(%hits)];

	if (@$clust > 1)
	{
	    push(@clusters,$clust);
	}
	$set = [grep { ! $hits{$_} } @$set];
    }
    return @clusters;
}
	
sub length_of {
    my($fig,$peg) = @_;

    my $tran = $fig->get_translation($peg);
    return $tran ? length($tran) : 0;
}

sub filter_out_pegs {
    my ($state) = @_;
    my($fig,$tmpdir,$peg_to_fam,$members_in_fam,$peg_to_ss_fam,$next_famP) = @$state;

    foreach my $fam (sort keys(%$members_in_fam)){
	my $set = $members_in_fam->{$fam};
	$set || confess "bad fam: $fam";

	foreach my $peg (@$set){
	    my $func = &function_of_peg($fig,$peg);
	    if ($func =~ /\#.*interrupted|\#.*truncated|\#.*frameshift|\#.*fragment/) {
		delete $peg_to_fam->{$peg};
	    }
	}
    }

}

sub delete_if_necessary {
  my($state) = @_;
  my($fig,$tmpdir,$peg_to_fam,$members_in_fam,$peg_to_ss_fam,$next_famP) = @$state;
  
  foreach my $fam (sort keys(%$members_in_fam)){
    my $fam_func =  &function_of_family($fig,$fam,$members_in_fam);
    next if ($fam_func =~ /not subsystem-based/);
    my $orgs_in_fam = {};
    my $set = $members_in_fam->{$fam};
    $set || confess "bad fam: $fam";

    # group pegs of same genome
    foreach my $peg (@$set){
      push @{$orgs_in_fam->{$fig->genome_of($peg)}}, $peg;
    }

    # delete figfam if the function of pegs in the same group have different functions and in a subsystem
    foreach my $group (keys %$orgs_in_fam){
      my $functionHash = {};
      foreach my $peg (@{$orgs_in_fam->{$group}}){
	my $func = &function_of_peg($fig,$peg);
	$functionHash->{$func} = 1;
      }

      if (scalar (keys %$functionHash) > 1){
	foreach my $peg (@{$orgs_in_fam->{$group}}){
	  my $func = &function_of_peg($fig,$peg);
	  if ($func ne $fam_func){
	    print STDERR "deleting $peg in $fam\n";
	    delete $peg_to_fam->{$peg};
	  }
	}
      }
    }
  }
}

sub split_if_necessary {
  my($state) = @_;
  my($fig,$tmpdir,$peg_to_fam,$members_in_fam,$peg_to_ss_fam,$next_famP) = @$state;
  
  foreach my $fam (sort keys(%$members_in_fam))
    {
      my %ffunc;
      my $set = $members_in_fam->{$fam};
      $set || confess "bad fam: $fam";
      
      foreach my $peg (@$set){
	my $func = &function_of_peg($fig,$peg);
	my @subs = grep { $fig->usable_subsystem($_) } $fig->peg_to_subsystems($peg);
	if (@subs > 0) {
	  $ffunc{$func}++;
	}
      }
      
      if (keys(%ffunc) > 1)
      {
	  &split_fam($state,$fam,[sort { $ffunc{$b} <=> $ffunc{$a} } keys(%ffunc)]);
      }
    }
}

sub split_fam {
    my($state,$fam,$funcs) = @_;
    my($fig,$tmpdir,$peg_to_fam,$members_in_fam,$peg_to_ss_fam,$next_famP) = @$state;

    my $pegs = $members_in_fam->{$fam};
    print STDERR "splitting $fam\n";
#   if ($fam eq 'FIG010670') { print STDERR &Dumper($pegs,$funcs) }
    my($i,$main_ffunc);
    for ($i=0; ($i < @$funcs); $i++)
    {
	my $func = $funcs->[$i];
	my @pegs_with_func = grep { &function_of_peg($fig,$_) eq $func } @$pegs;
#	if ($fam eq 'FIG010670') { print STDERR &Dumper([\@pegs_with_func,$func]) }

	if (@pegs_with_func > 1)
	{
	    my $new_fam;
	    if ($i == 0)
	    {
		$new_fam = $fam;
		$main_ffunc = $func;
	    }
	    else
	    {
		$new_fam = $$next_famP++;
		print STDERR "\tmaking $new_fam\n";
	    }

	    $members_in_fam->{$new_fam} = [@pegs_with_func];
#	    if ($fam eq 'FIG010670') { print STDERR &Dumper([$new_fam,$members_in_fam->{$new_fam}]) }
	    foreach my $peg (@pegs_with_func)
	    {
		$peg_to_fam->{$peg} = $new_fam;
#		if ($fam eq 'FIG010670') { print STDERR "changing connection of $peg from $fam to $new_fam\n"; }
	    }
	}
	else
	{
	    if ($i == 0)
	    {
		delete $members_in_fam->{$fam};
	    }
	    delete $peg_to_fam->{$pegs_with_func[0]};
#	    if ($fam eq 'FIG010670') { print STDERR &Dumper(['deleted peg to fam',$pegs_with_func[0]]) }
	}
    }

    # Now delete any remaining PEGs that were not in subsystems and had different functions; 
    # they cannot be placed
    foreach my $peg (@$pegs)
    {
	if ( ($peg_to_fam->{$peg} && ($peg_to_fam->{$peg} eq $fam)) &&
	     (!$peg_to_ss_fam->{$peg}) && (&function_of_peg($fig,$peg) ne $main_ffunc) )
	{
#	    if ($fam eq 'FIG010670') { print STDERR "disconnecting $peg from $fam\n" }
	    delete $peg_to_fam->{$peg};
	}
    }
}

sub add_new_if_necessary {
    my($state) = @_;
    my($fig,$tmpdir,$peg_to_fam,$members_in_fam,$peg_to_ss_fam,$next_famP) = @$state;

    my($peg,$peg1,$proposed,$i,$fam);
    foreach $peg (keys(%$peg_to_ss_fam))
    {
	if (! $peg_to_fam->{$peg})
	{
	    $proposed = $peg_to_ss_fam->{$peg};
	    for ($i=0; ($i < @$proposed) && (! $peg_to_fam->{$proposed->[$i]}); $i++) {}
	    if ($i == @$proposed)
	    {
		$fam = $$next_famP++;
		print STDERR "adding $fam\n";
		$members_in_fam->{$fam} = $proposed;
		foreach $peg1 (@$proposed)
		{
		    $peg_to_fam->{$peg1} = $fam;
		}
	    }
	    else
	    {
		$fam = $peg_to_fam->{$proposed->[$i]};
		foreach $peg1 (@$proposed)
		{
		    if (! $peg_to_fam->{$peg1})
		    {
			$peg_to_fam->{$peg1} = $fam;
			push(@{$members_in_fam->{$fam}},$peg1);
		    }
		}
	    }

	}
    }
}

sub merge_if_possible {
    my($state) = @_;
    my($fig,$tmpdir,$peg_to_fam,$members_in_fam,$peg_to_ss_fam,$next_famP) = @$state;

    my(%seen,$peg,$peg1,$fam1,$proposed,%fams,$i);

    foreach $peg (keys(%$peg_to_ss_fam))
    {
	if ((! $seen{$peg}) && ($proposed = $peg_to_ss_fam->{$peg}))
	{
	    undef %fams;
	    foreach $peg1 (@$proposed)
	    {
		$seen{$peg1} = 1;
		if ($fam1 = $peg_to_fam->{$peg1})
		{
		    $fams{$fam1} = 1;
		}
	    }

	    my @fams = sort keys(%fams);
	    if (@fams > 1)
	    {
		print STDERR "merging ",join(",",@fams)," into $fams[0]\n";

		my $fam = $fams[0];
		my $members = $members_in_fam->{$fam};
		for ($i=1; ($i < @fams); $i++)
		{
		    my $old = $members_in_fam->{$fams[$i]};
#		    if ($fams[$i] eq 'FIG010670') { print STDERR &Dumper([$fams[$i],$old]) }

		    foreach $peg1 (@$old)
		    {
			push(@$members,$peg1);
			$peg_to_fam->{$peg1} = $fam;
		    }
#		    if ($fams[$i] eq 'FIG010670') { print STDERR &Dumper(['deleting members in fam',$fams[$i]]) }
		    delete $members_in_fam->{$fams[$i]};
		}
	    }
	}
    }
}
    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3