PAM1 := MutationMatrix:
lnPAM1 := ln(MutationMatrix):
k := 1:
RateMutation := sum(Freq[i] * (1 - PAM1[i,i]), i= 1..20):
iter := 0:
# The following loop goes until the difference reach the computer's
# precision.
while abs(RateMutation - 0.01) > DBL_EPSILON do
iter := iter + 1;
k:= k * RateMutation * 100:
lnPAM1 := lnPAM1 / (RateMutation * 100):
PAM1 := exp(lnPAM1):
RateMutation := sum(Freq[i] * (1 - PAM1[i,i]), i= 1..20);
printf('After %d iteration(s), k=%f, RateMutation=%f
', iter, k, RateMutation);
od:
PrintMatrix(transpose(transpose(PAM1[1..10])[1..10]),'% .5f ');
####
the mutation matrix i used was:
Mutationmatrix
0.5447 0.0250 0.0311 0.0298 0.0807 0.0426 0.0497 0.0517 0.0203 0.0211 0.0251 0.0356 0.0328 0.0138 0.0557 0.1110 0.0554 0.0073 0.0148 0.0580
0.0151 0.5685 0.0275 0.0129 0.0168 0.0571 0.0261 0.0118 0.0351 0.0074 0.0101 0.1103 0.0173 0.0072 0.0120 0.0240 0.0203 0.0121 0.0114 0.0080
0.0141 0.0207 0.4755 0.0592 0.0135 0.0280 0.0255 0.0254 0.0452 0.0057 0.0050 0.0339 0.0077 0.0050 0.0099 0.0385 0.0295 0.0048 0.0152 0.0061
0.0202 0.0145 0.0887 0.5687 0.0067 0.0386 0.1123 0.0206 0.0346 0.0028 0.0043 0.0286 0.0050 0.0039 0.0175 0.0380 0.0221 0.0048 0.0091 0.0038
0.0120 0.0041 0.0044 0.0015 0.5719 0.0018 0.0009 0.0027 0.0031 0.0047 0.0036 0.0019 0.0046 0.0061 0.0016 0.0114 0.0100 0.0085 0.0061 0.0112
0.0174 0.0389 0.0253 0.0234 0.0050 0.4155 0.0531 0.0077 0.0356 0.0054 0.0119 0.0503 0.0237 0.0050 0.0126 0.0238 0.0190 0.0073 0.0114 0.0089
0.0404 0.0352 0.0457 0.1346 0.0050 0.1052 0.5182 0.0184 0.0320 0.0062 0.0081 0.0668 0.0132 0.0036 0.0243 0.0395 0.0354 0.0073 0.0106 0.0141
0.0482 0.0182 0.0523 0.0283 0.0168 0.0176 0.0212 0.7729 0.0153 0.0029 0.0043 0.0221 0.0064 0.0047 0.0120 0.0428 0.0144 0.0097 0.0045 0.0054
0.0050 0.0143 0.0245 0.0125 0.0050 0.0213 0.0097 0.0040 0.5384 0.0022 0.0046 0.0149 0.0073 0.0149 0.0089 0.0089 0.0074 0.0157 0.0364 0.0037
0.0169 0.0099 0.0102 0.0033 0.0252 0.0106 0.0061 0.0025 0.0071 0.4988 0.1017 0.0117 0.0901 0.0398 0.0063 0.0102 0.0288 0.0145 0.0140 0.1760
0.0264 0.0176 0.0116 0.0066 0.0252 0.0304 0.0104 0.0048 0.0198 0.1330 0.5999 0.0201 0.1866 0.0914 0.0162 0.0194 0.0245 0.0242 0.0349 0.0803
0.0273 0.1405 0.0575 0.0324 0.0101 0.0942 0.0630 0.0182 0.0468 0.0112 0.0147 0.4960 0.0218 0.0064 0.0220 0.0402 0.0317 0.0048 0.0121 0.0114
0.0090 0.0079 0.0047 0.0020 0.0084 0.0158 0.0044 0.0019 0.0081 0.0307 0.0486 0.0078 0.4201 0.0226 0.0031 0.0071 0.0118 0.0060 0.0091 0.0181
0.0062 0.0054 0.0050 0.0026 0.0185 0.0055 0.0020 0.0023 0.0275 0.0223 0.0392 0.0037 0.0373 0.6051 0.0037 0.0071 0.0072 0.0713 0.1391 0.0137
0.0265 0.0095 0.0105 0.0123 0.0050 0.0146 0.0143 0.0061 0.0173 0.0037 0.0073 0.0136 0.0055 0.0039 0.7212 0.0274 0.0122 0.0024 0.0038 0.0073
0.0743 0.0267 0.0570 0.0375 0.0513 0.0389 0.0325 0.0307 0.0244 0.0085 0.0123 0.0351 0.0173 0.0105 0.0385 0.4331 0.0904 0.0157 0.0106 0.0135
0.0374 0.0228 0.0441 0.0221 0.0454 0.0313 0.0295 0.0104 0.0203 0.0242 0.0158 0.0279 0.0291 0.0108 0.0173 0.0912 0.5182 0.0048 0.0118 0.0405
0.0007 0.0021 0.0011 0.0007 0.0059 0.0018 0.0009 0.0011 0.0066 0.0019 0.0024 0.0006 0.0023 0.0163 0.0005 0.0024 0.0007 0.7037 0.0182 0.0019
0.0049 0.0062 0.0110 0.0044 0.0135 0.0091 0.0043 0.0016 0.0488 0.0057 0.0109 0.0052 0.0109 0.1014 0.0026 0.0052 0.0057 0.0580 0.6111 0.0057
0.0533 0.0122 0.0124 0.0051 0.0698 0.0201 0.0160 0.0053 0.0137 0.2016 0.0703 0.0136 0.0610 0.0279 0.0141 0.0186 0.0552 0.0169 0.0159 0.5124
my @PAM1= @mutMatrix;
my $e; # e is equal to the DBL_EPSILON constant
$e=1;$e/=2 while $e/2+1>1;
my $k = 1;
my $mutRate = sum map{ $freq[ $_ ] * ( 1 - $mutMatrix[ $_ ][ $_ ] ) } 0..19;
my $iter = 0;
my @lnPAM=(); # in case the position in the matrix is 0
LINE: for (my $i=0; $i<20; $i++){
for (my $j=0; $j<20; $j++){
if ($mutMatrix[$i][$j] > 0){
$lnPAM[$i][$j] = log $mutMatrix[$i][$j];
}else{
next LINE;
}
}
}
# The following loop goes until the difference reach the computer's
# precision.
while( abs( $mutRate - 0.01) > $e ) {
++$iter;
$k *= $mutRate * 100;
map $_ /= ( $mutRate * 100), @$_ for @lnPAM;
@PAM1 = map[ map exp( $_ ), @$_ ], @lnPAM;
$mutRate = sum map{ $freq[ $_ ] * ( 1 - $PAM1[ $_ ][ $_ ] ) } 0..19;
printf "After %d iteration(s), k=%f, RateMutation=%f\n", $iter, $k, $mutRate;
}
print "\n";
print "PAM1 Matrix:\n\n";
for (my $i=0; $i<20; $i++){ #prints Pam1 matrix print PAM1 matrix
for (my $j=0; $j<20; $j++){
printf("%.4f ", $PAM1[$i][$j]);
}
print "\n";
}
####
0.99002 0.00040 0.00048 0.00048 0.00193 0.00090 0.00115 0.00112 0.00031 0.00018
0.00024 0.99072 0.00052 0.00013 0.00035 0.00140 0.00045 0.00019 0.00075 0.00011
0.00022 0.00039 0.98805 0.00162 0.00028 0.00065 0.00046 0.00057 0.00121 0.00011
0.00033 0.00015 0.00242 0.99064 0.00006 0.00071 0.00311 0.00036 0.00071 0.00002
0.00029 0.00009 0.00009 0.00001 0.99123 0.00001 -0.00001 0.00004 0.00005 0.00006
0.00037 0.00095 0.00059 0.00043 0.00003 0.98585 0.00157 0.00011 0.00095 0.00004
0.00094 0.00060 0.00083 0.00373 -0.00005 0.00312 0.98900 0.00029 0.00057 0.00004
0.00104 0.00030 0.00117 0.00050 0.00027 0.00026 0.00034 0.99589 0.00024 0.00002
0.00008 0.00031 0.00065 0.00026 0.00008 0.00057 0.00017 0.00006 0.99023 0.00001
0.00015 0.00014 0.00019 0.00003 0.00032 0.00008 0.00004 0.00002 0.00005 0.98765
####
0.9876 0.9712 0.9670 0.9662 0.9586 0.9726 0.9654 0.9654 0.9677 0.9605 0.9644 0.9670 0.9694 0.9646 0.9607 0.9717 0.9710 0.9700 0.9714 0.9615
0.9607 0.9891 0.9504 0.9522 0.9552 0.9458 0.9567 0.9588 0.9543 0.9515 0.9513 0.9513 0.9462 0.9505 0.9555 0.9614 0.9589 0.9572 0.9592 0.9606
0.9663 0.9601 0.9908 0.9770 0.9504 0.9667 0.9638 0.9548 0.9665 0.9470 0.9532 0.9720 0.9621 0.9655 0.9577 0.9690 0.9657 0.9516 0.9567 0.9551
0.9682 0.9647 0.9798 0.9895 0.9601 0.9697 0.9711 0.9587 0.9716 0.9557 0.9605 0.9726 0.9659 0.9751 0.9662 0.9702 0.9716 0.9592 0.9517 0.9613
0.9591 0.9662 0.9515 0.9585 0.9941 0.9426 0.9711 0.9630 0.9565 0.9682 0.9668 0.9533 0.9549 0.9481 0.9547 0.9566 0.9609 0.9612 0.9671 0.9827
0.9716 0.9552 0.9665 0.9666 0.9412 0.9932 0.9567 0.9504 0.9664 0.9477 0.9520 0.9676 0.9592 0.9620 0.9638 0.9683 0.9636 0.9572 0.9478 0.9502
0.9579 0.9596 0.9571 0.9614 0.9630 0.9501 0.9848 0.9564 0.9631 0.9546 0.9538 0.9627 0.9525 0.9685 0.9615 0.9595 0.9582 0.9528 0.9412 0.9730
0.9676 0.9715 0.9577 0.9587 0.9647 0.9534 0.9661 0.9883 0.9649 0.9766 0.9746 0.9635 0.9659 0.9656 0.9645 0.9634 0.9673 0.9808 0.9618 0.9615
0.9691 0.9662 0.9687 0.9709 0.9574 0.9687 0.9721 0.9641 0.9863 0.9570 0.9629 0.9726 0.9692 0.9758 0.9785 0.9709 0.9685 0.9644 0.9346 0.9646
0.9685 0.9701 0.9557 0.9617 0.9759 0.9566 0.9702 0.9826 0.9636 0.9940 0.9859 0.9631 0.9647 0.9637 0.9651 0.9659 0.9710 0.9775 0.9584 0.9732
0.9579 0.9552 0.9474 0.9519 0.9598 0.9464 0.9547 0.9658 0.9549 0.9711 0.9870 0.9502 0.9552 0.9527 0.9536 0.9551 0.9593 0.9673 0.9412 0.9575
0.9694 0.9642 0.9752 0.9729 0.9552 0.9709 0.9727 0.9637 0.9736 0.9575 0.9590 0.9896 0.9674 0.9733 0.9678 0.9724 0.9717 0.9598 0.9451 0.9653
0.9672 0.9545 0.9607 0.9617 0.9524 0.9581 0.9579 0.9616 0.9657 0.9546 0.9596 0.9630 0.9927 0.9622 0.9630 0.9675 0.9637 0.9578 0.9643 0.9566
0.9638 0.9601 0.9654 0.9722 0.9468 0.9622 0.9753 0.9626 0.9735 0.9548 0.9584 0.9701 0.9634 0.9876 0.9708 0.9692 0.9661 0.9609 0.9634 0.9602
0.9625 0.9679 0.9603 0.9660 0.9561 0.9666 0.9710 0.9642 0.9789 0.9589 0.9620 0.9673 0.9670 0.9736 0.9918 0.9660 0.9660 0.9624 0.9478 0.9654
0.9770 0.9772 0.9749 0.9733 0.9613 0.9745 0.9722 0.9664 0.9747 0.9630 0.9668 0.9753 0.9748 0.9752 0.9693 0.9886 0.9784 0.9683 0.9648 0.9691
0.9718 0.9702 0.9673 0.9704 0.9613 0.9654 0.9666 0.9659 0.9679 0.9637 0.9667 0.9701 0.9666 0.9678 0.9650 0.9740 0.9873 0.9697 0.9674 0.9615
0.9726 0.9704 0.9550 0.9598 0.9633 0.9608 0.9630 0.9813 0.9656 0.9720 0.9766 0.9600 0.9625 0.9643 0.9632 0.9657 0.9715 0.9900 0.9592 0.9661
0.9561 0.9545 0.9424 0.9348 0.9514 0.9338 0.9337 0.9446 0.9185 0.9355 0.9327 0.9279 0.9513 0.9491 0.9311 0.9445 0.9514 0.9416 0.9950 0.9621
0.9567 0.9664 0.9512 0.9545 0.9774 0.9464 0.9758 0.9546 0.9584 0.9603 0.9593 0.9582 0.9539 0.9563 0.9588 0.9592 0.9559 0.9587 0.9726 0.9888