view release on metacpan or  search on metacpan

lib/Bio/BPWrapper/TreeManipulations.pm  view on Meta::CPAN

	    $sites->{$otu}->{$colnames[$i]} = $data[$i] ? 1:0; # force binary
	}
    }
    close BIN;

     foreach my $otu (@otus) {
#    die "otu id not found in rownames\n" unless &_check_id($node->id, \@rownames);
	foreach my $trait_name (@colnames) {
	    $otu->add_tag_value($trait_name, [$sites->{$otu->id}->{$trait_name}]);
	}
    }
#    print Dumper(\@otus); exit;

# Fitch algorithm: post-order traversal
    my @informative;
    for (my $i=0; $i<=$#colnames; $i++) {
	next unless &_is_informative($colnames[$i], $i);
	push @informative, $colnames[$i];
	$rootnode->add_tag_value($colnames[$i], &_fitch_parsimony($rootnode, $i, \@colnames));
	&_penny_parsimony($rootnode, $i, \@colnames) if @{$rootnode->get_tag_values($colnames[$i])} == 2; # only when root state unresolved
	my $ci = &_consistency_index($i, \@colnames);
	print join "\t", ($i+1, $colnames[$i], $ci);
	print "\n";
    }
}

sub _fitch_parsimony {
    my ($node,$index, $refcol)=@_; #warn $node->internal_id, "\t", $node->id() || "inode", "\n";
    my $ref_node_state;
    my @colnames = @{$refcol};
    if ($node->is_Leaf) {
#       warn $node->id, "\t", Dumper($node->get_tag_values($colnames[$index])), "\n";
        return $node->get_tag_values($colnames[$index]);
    } else {
        my @child = $node->each_Descendent;
        my ($ref0, $ref1);
        if ($child[0]->is_Leaf) { # child 0 is an OTU
            $ref0 = $child[0]->get_tag_values($colnames[$index]);
            if ($child[1]->is_Leaf) { # both child 0 & 1 are an OTU
                $ref1 = $child[1]->get_tag_values($colnames[$index]);
#               warn "got sis otu for inode ", $node->internal_id(), "\n";
                $ref_node_state = &_intersect_or_union($ref0, $ref1);
#               warn Dumper($node->get_tag_values($colnames[$index]));
            } else { # child 0 is an OTU, child 1 is an inode
                $ref_node_state = &_intersect_or_union($ref0, &_fitch_parsimony($child[1], $index, \@colnames));
            }
        } else { # child 0 is an inode
            if ($child[1]->is_Leaf) { # child 1 is an inode child 1 is an OTU
                $ref1 = $child[1]->get_tag_values($colnames[$index]);
#               warn "got sis otu for inode ", $node->internal_id(), "\n";
                $ref_node_state = &_intersect_or_union(&_fitch_parsimony($child[0], $index, \@colnames), $ref1);
            } else { # both inodes
                $ref_node_state = &_intersect_or_union(&_fitch_parsimony($child[0], $index, \@colnames), &_fitch_parsimony($child[1], $index, \@colnames));
            }
        }
        $node->add_tag_value($colnames[$index], $ref_node_state);
        return $ref_node_state;
    }
}

sub _intersect_or_union { # from Perl Cookbook
    my ($ref1, $ref2) = @_;
    my (%union, %isect);
    foreach my $e (@$ref1, @$ref2) {
        $union{$e}++ && $isect{$e}++; # Perl Cookbook ideom
    }
 #   warn Dumper(\%union, \%isect);

    if (@$ref1 == @$ref2) { # (0) U (0); (0) U (1); (1) U (1); (0,1) U (0,1)
        return [keys %union];
    } else { # (0) I (0,1); (1) I (0,1)
        return [keys %isect];
    }
}

sub _penny_parsimony {
    my $nd = shift;
    my $index = shift;
    my $refcol = shift;
    my @colnames = @{$refcol};
    my $ref_states = $nd->get_tag_values($colnames[$index]);
    my $ref_new;
    my $pa;
    return if $nd->is_Leaf;
    if ($pa = $nd->ancestor()) {
        my $ref_pa_state = $pa->get_tag_values($colnames[$index]);
        $ref_new = &_intersect_or_union($ref_states, $ref_pa_state); # intersect with with parent
    } else { # is the root
        $ref_new = [1]; # resolve root state using Penny parsimony: one gain only
    }
    $nd->add_tag_value($colnames[$index], $ref_new); # resolve root state using Penny parsimony: one gain only
    &_penny_parsimony($_, $index, \@colnames) for $nd->each_Descendent();
}

sub _is_informative {
    my $col_id = shift;
    my $index = shift;
    my @sts = map {$_->get_tag_values($col_id)->[0]} @otus;
    my %seen;
    $seen{$_}++ for @sts;
    if (keys %seen == 1) { #warn "$index: $col_id is constant\n";
        return 0 }
    foreach (keys %seen) {
        if ($seen{$_} == 1) { #warn "$index: $col_id is singleton\n";
            return 0 }
    }
#   warn "$index: $col_id is informative\n";
    return 1;
}

sub _consistency_index { # ratio of minimum (1 for a binary trait) to actural changes
    my $index = shift;
    my $ci = 0;
    my $refcol = shift;
    my @colnames = @{$refcol};
    foreach my $nd (@nodes) {
        next if $nd eq $rootnode;
        my $pa = $nd->ancestor();
        my $state = $nd->get_tag_values($colnames[$index])->[0]; # assuming fully resolved (only one state)
        my $pa_state = $pa->get_tag_values($colnames[$index])->[0];
        if ($state ne $pa_state) {
#            if ($state) {
#                push @{$gain_loss{'gain'}->{$nd->internal_id}}, $colnames[$index];
#                $fam_events{$colnames[$index]}->{$nd->internal_id}->{'gain'}++;
#            } else {