# this heuristic seems to be faulty
        my $min_mm = scalar( uniq(@colbp) ) - 1;
[download]

To me this seems perfectly right. If you want to speed up things a bit more (this is only a micro optimization) you can use a hash in the first place:

        my %chars;
        foreach my $site ( @{$tfbs} ) {
            my $bp = substr($site,$pos,1);
            $chars{$bp} = 1;

        }
        
       # this heuristic seems to be faulty
        my $min_mm = keys %chars - 1;
        push @mincol, $min_mm;
[download]

Instead I think that your sample output is wrong:

# is: 
00100000302011000000100
# should be:
00100000301011000000100
[download]

If not, could you please explain how to get the 2 there?

Hi Moritz,
It is based on phylogenetic tree, the substitution comparison from lower nodes must be compared to higher level.

E.g. For column 11, the bases for row 3 is compared to row 1. Same result will also be given if we compare row 3 to row 2. Both row 1 and row 2 has higher level than row 3 row 4 in phylogenetic tree.


---
neversaint and everlastingly indebted.......

I still don't quite understand. We have

human:  C
chimp:  T
mouse:  T
rat:    C
[download]

So which comparisons exactly lead to to a number of 2? With which rules?

As a poor perl programmer I'm not familiar with the phylogenetic tree, so if that's relevant here please explain in which relation these four examples appear in that tree.

Have you changed your example?

A
C
T
G
[download]

How does that become 3? And the column moritz quoted above used to bee:

C
T
T
T
[download]

If you do want to work on the code yourself bioperl may be a good place to look. I myself have not worked with the tools there, but I get the sense others have done similar things to what you are trying to do. So there maybe something there that is helpful for you.

There are also numerous programs that grappled with Fitch's algorithm (and other algorithms for that matter), that might avoid the trouble of reinventing the wheel from scratch. For some helpful links to programs you might want to check out the Workshop on Molecular Evolution.

Sorry I can't help more than that.
Bio.

I think that the idea is that given the tree:

the values in column 11 should be treated like this: (T,C,(C,T)) and so the value of 2 results from T ne C ne (C,T); whereas for instance column 3 leads to (G,A,(A,A)) and G ne A eq (A,A).

This is just a wild guess but, if I understand the logic right, (T,C,(A,G))=>3 and (T,C,(C,A))=>2. If it is so, the algorithm should not be so simple when more complex trees are involved.

Just a simple question: am I wrong or, from a cladistic point of view, any node of the tree should have exactly zero or two descendant, one of which is a terminal node?

For if this is the case the problem becomes quite a trivial one, but your example doesn't match my hypothesis.

Update: forget the second condition, it's bullshit.

But I think that the first one holds, for you can have ((rat,mouse),(human,chimp)) or ((rat,(mouse,(human,chimp)) or (mouse,(rat,(human,chimp)) but it make no sense (rat,mouse,(human,chimp)), particularly if you are checking the variation score independently on each base. For the total score to make sense you should assume one more common ancestor or the score for base X will not be consistent with that of base Y.

Dear psini,

Thank you so much for your thoughtful replies.

You are right, that the phylogenetic tree is always a binary tree.
I guess the problem I gave you is the tree "representation". Newick's format gives a general-ordered tree.
And it seems that we always can convert this general-ordered tree into a binary tree.



---
neversaint and everlastingly indebted.......

Dear dwm042,

You are right. It is a small parsimony problem.
BTW, from the start I have already included the "genetic tree input" in my sample code.


---
neversaint and everlastingly indebted.......

Assuming (as in my previous post) that the phylogenetic tree is a binary tree, the following code should work.

use strict;
use warnings;
use Data::Dumper;

my $sites={'mouse'=>'CGGCGGAAAACTGTCCTCCGTGC',
           'rat'=>  'CGACGGAACATTCTCCTCCGCGC', 
           'human'=>'CGACGGAATATTCCCCTCCGTGC',
           'chimp'=>'CGACGGAAGACTCTCCTCCGTGC'};
my $tree = [['mouse','rat'],['human','chimp']]; 

# Sanity checks
my $len=-1;
treeCheck($tree,$sites,\$len);
              
my @val;
$val[$_]=parsimony(extractLocus($tree,$sites,$_))->[0] for (0..$len-1)
+;
print join('',@val)."\n";


sub treeCheck {
  my ($subTree,$sites,$len)=@_;
  if (ref($subTree) eq 'ARRAY') {
    die "Not a binary tree" if @$subTree != 2;
    treeCheck($subTree->[0],$sites,$len);
    treeCheck($subTree->[1],$sites,$len);
  } else {
    die "Invalid char" if $sites->{$subTree} !~ /^[ATCG]+$/;
    if ($$len<0) {
      $$len=length($sites->{$subTree});
    } else {
      die "Invalid length" if length($sites->{$subTree}) != $$len;
    }
  }
}

sub extractLocus {
  my ($subTree,$sites,$index)=@_;
  if (ref($subTree) eq 'ARRAY') {
    return [extractLocus($subTree->[0],$sites,$index),extractLocus($su
+bTree->[1],$sites,$index)];
  } else {
    my $base=substr($sites->{$subTree},$index,1);
    $base =~ tr/ATCG/1248/;    
    return $base;
  }
}
  
sub parsimony {
  my $locusTree=shift;
  if (ref($locusTree) eq 'ARRAY') {
    my $left=parsimony($locusTree->[0]);
    my $right=parsimony($locusTree->[1]);
    my $count=$left->[0]+$right->[0];
    $count++ if ($left->[1] & $right->[1])==0;
    my $bases=$left->[1] | $right->[1];
    return [$count,$bases]
  } else {
    return[0,$locusTree];
  }
}
[download]

I have significantly changed the input data structures, so $tree is "really" a tree (an AoA..oA) containing the names of the species and $sites is an hash using the species as key. It works with your original data but changing the structure of the tree gives (as to be expected) different results:

With $tree = [['mouse','rat'],['human','chimp']]; 
gives 00100000302011000000100

With $tree = ['mouse',['rat',['human','chimp']]]; 
gives 00100000301011000000100
[download]

I have no understanding of the problem domain so the following code is not a solution to your actual problem, but is an alternative way of finding changes between adjacent strings in a list and counting those changes by column.

use strict;
use warnings;

my $sites = [
    'CGGCGGAAAACTGTCCTCCGTGC',    # mouse
    'CGACGGAACATTCTCCTCCGCGC',    # rat
    'CGACGGAATATTCCCCTCCGTGC',    # human
    'CGACGGAAGACTCTCCTCCGTGC',    # chimp
    ];

my @val = my_parsimony ($sites);
print @val;


sub my_parsimony {
    my ($sites) = @_;
    my $lastSite = $sites->[0];
    my @mincol = (0) x length $lastSite;

    for my $nextSite (@$sites) {
        my $changes = $lastSite ^ $nextSite;
        
        $lastSite = $nextSite;
        next unless $changes =~ s/[^\0]/+/g;

        my $pos = 0;

        ++$mincol[$pos++] while -1 != ($pos = index $changes, '+', $po
+s);
    }

    return @mincol;
}
[download]

Prints:

00100000302012000000200
[download]