@madM , I had a look at the code. I noticed that Freq was missing from it, so there was no way to reproduce the numbers you described aboe. However, after a bit of googling, I found this page which refered to this language called Darwin and how you can compute mutation matrices and mutation rates and it also included the exact code that computes the 1-PAM matrix which madM also wrote above.
FWIW this language called Darwin is closed source. I found that one can obtain a copy if he signs a license agreement and sends it by fax or regular mail!.
Anyway, I had a look over the page and I noticed that one can obtain the MutationMatrix from the CountMatrix. Afterwards, one can find the Freq frequencies through the CountMatrix. So .. what I did is I used the CountMatrix featured in the page mentioned above, and I wrote some PDL code to compute the various matrices. The results weren't anywhere near the ones in the articles. Right off the bat the Mutation matrix was different, the frequencies were surprisingly the same as in the article (and without a copy of this Darwin programming language I can't check more than this). Maybe someone wants to have a second look(although I'll mention that I'd much rather trust PDL). I'll just post the Perl PDL code here:
#!/usr/bin/env perl
use strict;
use warnings;
use Data::Dumper;
use PDL;
use PDL::NiceSlice;
use Storable qw/dclone/;
my $CountMatrix = pdl([
[ 4372, 121, 113, 162, 96, 140, 324, 387, 40, 136, 212, 21
+9, 72, 50, 213, 596, 300, 6, 39, 428],
[ 121, 2748, 100, 70, 20, 188, 170, 88, 69, 48, 85, 67
+9, 38, 26, 46, 129, 110, 10, 30, 59],
[ 113, 100, 1727, 322, 16, 92, 166, 190, 89, 37, 42, 20
+9, 17, 18, 38, 207, 160, 4, 40, 45],
[ 162, 70, 322, 3092, 8, 127, 732, 154, 68, 18, 36, 17
+6, 11, 14, 67, 204, 120, 4, 24, 28],
[ 96, 20, 16, 8, 680, 6, 6, 20, 6, 30, 30, 1
+2, 10, 22, 6, 61, 54, 7, 16, 83],
[ 140, 188, 92, 127, 6, 1367, 346, 58, 70, 35, 100, 31
+0, 52, 18, 48, 128, 103, 6, 30, 66],
[ 324, 170, 166, 732, 6, 346, 3378, 138, 63, 40, 68, 41
+1, 29, 13, 93, 212, 192, 6, 28, 104],
[ 387, 88, 190, 154, 20, 58, 138, 5788, 30, 19, 36, 13
+6, 14, 17, 46, 230, 78, 8, 12, 40],
[ 40, 69, 89, 68, 6, 70, 63, 30, 1059, 14, 39, 9
+2, 16, 54, 34, 48, 40, 13, 96, 27],
[ 136, 48, 37, 18, 30, 35, 40, 19, 14, 3217, 858, 7
+2, 198, 144, 24, 55, 156, 12, 37, 1300],
[ 212, 85, 42, 36, 30, 100, 68, 36, 39, 858, 5061, 12
+4, 410, 331, 62, 104, 133, 20, 92, 593],
[ 219, 679, 209, 176, 12, 310, 411, 136, 92, 72, 124, 305
+4, 48, 23, 84, 216, 172, 4, 32, 84],
[ 72, 38, 17, 11, 10, 52, 29, 14, 16, 198, 410, 4
+8, 923, 82, 12, 38, 64, 5, 24, 134],
[ 50, 26, 18, 14, 22, 18, 13, 17, 54, 144, 331, 2
+3, 82, 2191, 14, 38, 39, 59, 367, 101],
[ 213, 46, 38, 67, 6, 48, 93, 46, 34, 24, 62, 8
+4, 12, 14, 2757, 147, 66, 2, 10, 54],
[ 596, 129, 207, 204, 61, 128, 212, 230, 48, 55, 104, 21
+6, 38, 38, 147, 2326, 490, 13, 28, 100],
[ 300, 110, 160, 120, 54, 103, 192, 78, 40, 156, 133, 17
+2, 64, 39, 66, 490, 2808, 4, 31, 299],
[ 6, 10, 4, 4, 7, 6, 6, 8, 13, 12, 20,
+4, 5, 59, 2, 13, 4, 582, 48, 14],
[ 39, 30, 40, 24, 16, 30, 28, 12, 96, 37, 92, 3
+2, 24, 367, 10, 28, 31, 48, 1612, 42],
[ 428, 59, 45, 28, 83, 66, 104, 40, 27, 1300, 593, 8
+4, 134, 101, 54, 100, 299, 14, 42, 3784],
]);
# $r1 is a row vector of 1 - M[i,i]
# $r2 is a column vector of values from $F (actually from $Freq)
# $r1 x $r2 is the product which and it's a scalar (actually a 1x1 mat
+rix, so that's
# why we're using ->at(0,0) to pull out that number from that 1x1 matr
+ix).
sub RateMutation {
my ($M,$F) = @_;
my $Mdim = shape($M)->at(1,0);
my $r1 = pdl(ones($Mdim) - $M->diagonal(0,1));
my $r2 = $F->(*1);
return pdl($r1 x $r2)->at(0,0);
}
# $Tot holds the total of all elements in the matrix $M
# $Mdim is the dimension of the matrix minus one (in our case 20 since
+ we have a 20 x 20 square matrix)
# That's because indexing starts at 0 for piddles (the primite PDL dat
+a structures)
sub ComputeFrequencies {
my ($M) = @_;
my $Tot = sum($M);
my $Mdim = shape($M)->at(1,0) - 1;
my $F = pdl(map { $M->slice(",$_")->sum /$Tot } (0..$Mdim));
return $F;
}
sub ComputeMutationMatrix {
my ($C) = @_;
my $N = zeroes(20,20);
for(0..19) {
$N->set($_,$_,1/$C->slice(":,$_")->sum );
};
my $M = $C x $N;
return $M;
};
my $MutationMatrix = ComputeMutationMatrix($CountMatrix);
my $Freq = ComputeFrequencies($CountMatrix);
my $DBL_EPSILON; # e is equal to the DBL_EPSILON constant
$DBL_EPSILON=1;$DBL_EPSILON/=2 while $DBL_EPSILON/2+1>1;
print "dbl_epsilon=$DBL_EPSILON\n";
my $PAM1 = $MutationMatrix->copy();
my $lnPAM1 = log($PAM1);
my $k = 1;
my $r = RateMutation($PAM1,$Freq);
while(abs($r -0.01) > $DBL_EPSILON) {
$k *= $r * 100;
$lnPAM1 /= $r * 100;
$PAM1 = $lnPAM1->exp;
$r = RateMutation($PAM1,$Freq);
print "r=$r\n";
};
# print the 10 x 10 upper left sub-matrix of $PAM1
# not sure if this is the same as
# the Darwin code ==> PrintMatrix(transpose(transpose(PAM1[1..10])[1
+..10]),'% .5f ');
print $PAM1->slice('0:9,0:9');
And, here are the results I got from this code above ^ :
Freq = [0.084746479 0.051042173 0.038350263 0.057409245 0.012554643 0.
+034739087 0.068834076 0.079076299 0.020769539 0.068105505 0.089075666
+ 0.06501172 0.023198108 0.038234114 0.040367031 0.056701793 0.0572191
+84 0.0087322873 0.027854624 0.077978164]
Mutation Matrix
[
[ 0.54472963 0.02503103 0.031112335 0.029795843 0.08074011
+8 0.042553191 0.049700874 0.051675791 0.020335536 0.0210852
+71 0.025130394 0.035569271 0.032771962 0.01380834 0.055715
+407 0.11098696 0.055360768 0.0072551391 0.014783927 0.05795
+5315]
[ 0.015076003 0.56847331 0.02753304 0.012874747 0.01682085
+8 0.057142857 0.026077619 0.011750567 0.0350788 0.00744186
+05 0.010075865 0.11028098 0.017296313 0.0071803369 0.012032
+435 0.024022346 0.020298948 0.012091898 0.011372252 0.007989
+1672]
[ 0.014079242 0.020686802 0.47549559 0.059223837 0.01345668
+6 0.027963526 0.025464028 0.025370543 0.045246568 0.00573643
+41 0.0049786629 0.033945103 0.0077378243 0.0049710025 0.0099398
+378 0.038547486 0.029525743 0.0048367594 0.015163002 0.006093
+4326]
[ 0.020184401 0.014480761 0.088656388 0.56869597 0.006728343
+1 0.038601824 0.11228716 0.020563493 0.034570412 0.00279069
+77 0.0042674253 0.02858535 0.0050068275 0.0038663353 0.017525
+504 0.037988827 0.022144307 0.0048367594 0.0090978014 0.003791
+4692]
[ 0.011961126 0.0041373604 0.0044052863 0.0014713997 0.5719091
+7 0.0018237082 0.00092038656 0.0026705835 0.0030503305 0.00465116
+28 0.0035561878 0.0019490011 0.0045516614 0.0060756697 0.0015694
+481 0.011359404 0.0099649382 0.0084643289 0.0060652009 0.01123
+8998]
[ 0.017443309 0.038891187 0.025330396 0.02335847 0.005046257
+4 0.41550152 0.053075625 0.0077446922 0.035587189 0.00542635
+66 0.011853959 0.050349196 0.023668639 0.0049710025 0.012555
+585 0.023836127 0.019007197 0.0072551391 0.011372252 0.008937
+0345]
[ 0.040368801 0.035167563 0.045704846 0.13463307 0.005046257
+4 0.10516717 0.51817763 0.018427026 0.03202847 0.00620155
+04 0.0080606923 0.066753289 0.013199818 0.0035901685 0.024326
+445 0.039478585 0.035430891 0.0072551391 0.010614102 0.014
+0826]
[ 0.048218291 0.018204386 0.052312775 0.028324444 0.01682085
+8 0.017629179 0.021168891 0.77286687 0.015251652 0.00294573
+64 0.0042674253 0.02208868 0.0063723259 0.0046948357 0.012032
+435 0.04283054 0.0143938 0.0096735187 0.0045489007 0.005416
+3846]
[ 0.0049838026 0.014273893 0.024504405 0.012506897 0.005046257
+4 0.021276596 0.0096640589 0.0040058753 0.53838332 0.00217054
+26 0.0046230441 0.014942342 0.0072826582 0.014913007 0.0088935
+391 0.0089385475 0.0073814357 0.015719468 0.036391205 0.003656
+0596]
[ 0.016944929 0.0099296649 0.010187225 0.0033106493 0.02523128
+7 0.010638298 0.0061359104 0.0025370543 0.0071174377 0.498759
+69 0.10170697 0.011694007 0.090122895 0.03976802 0.0062777
+923 0.010242086 0.028787599 0.014510278 0.014025777 0.176
+0325]
[ 0.026414154 0.017583782 0.011563877 0.0066212985 0.02523128
+7 0.030395137 0.010431048 0.0048070503 0.019827148 0.133023
+26 0.59992888 0.020139678 0.18661812 0.091411212 0.01621
+763 0.019366853 0.024543274 0.024183797 0.034874905 0.08029
+7901]
[ 0.027286319 0.14046338 0.057544053 0.032370793 0.01009251
+5 0.094224924 0.06304648 0.018159968 0.046771734 0.0111627
+91 0.014698909 0.49602079 0.021847975 0.0063518365 0.021972
+273 0.040223464 0.031740173 0.0048367594 0.012130402 0.01137
+4408]
[ 0.0089708448 0.0078609847 0.0046806167 0.0020231745 0.008410428
+9 0.015805471 0.0044485351 0.0018694085 0.0081342145 0.0306976
+74 0.048601233 0.0077960045 0.42011834 0.022645678 0.0031388
+962 0.0070763501 0.011810297 0.0060459492 0.0090978014 0.01814
+4888]
[ 0.0062297533 0.0053785685 0.0049559471 0.0025749494 0.01850294
+4 0.0054711246 0.0019941709 0.002269996 0.027452974 0.0223255
+81 0.039236605 0.0037355855 0.037323623 0.60508147 0.0036620
+455 0.0070763501 0.0071968998 0.071342201 0.13912055 0.01367
+6371]
[ 0.026538749 0.0095159288 0.010462555 0.012322972 0.005046257
+4 0.014589666 0.014265992 0.0061423421 0.017285206 0.00372093
+02 0.0073494547 0.013643008 0.0054619936 0.0038663353 0.72116
+139 0.027374302 0.012179369 0.0024183797 0.0037907506 0.007312
+1192]
[ 0.074258659 0.026685974 0.056993392 0.037520692 0.05130361
+6 0.038905775 0.032520325 0.030711711 0.024402644 0.00852713
+18 0.012328118 0.03508202 0.017296313 0.010494339 0.038451
+478 0.43314711 0.090422587 0.015719468 0.010614102 0.01354
+0961]
[ 0.03737852 0.022755482 0.044052863 0.022070995 0.04541631
+6 0.031306991 0.02945237 0.010415276 0.020335536 0.0241860
+47 0.015765766 0.027935683 0.029130633 0.010770505 0.017263
+929 0.091247672 0.51817679 0.0048367594 0.011751327 0.04048
+7475]
[ 0.0007475704 0.0020686802 0.0011013216 0.00073569983 0.005887300
+3 0.0018237082 0.00092038656 0.0010682334 0.0066090493 0.00186046
+51 0.0023707918 0.00064966705 0.0022758307 0.016293841 0.00052314
+936 0.0024208566 0.00073814357 0.70374849 0.018195603 0.001895
+7346]
[ 0.0048592076 0.0062060405 0.011013216 0.004414199 0.01345668
+6 0.009118541 0.0042951373 0.0016023501 0.048805287 0.00573643
+41 0.010905642 0.0051973364 0.010923987 0.10135322 0.0026157
+468 0.0052141527 0.0057206127 0.058041112 0.61106899 0.005687
+2038]
[ 0.053326688 0.012205213 0.012389868 0.0051498988 0.0698065
+6 0.02006079 0.015953367 0.005341167 0.013726487 0.201550
+39 0.070293978 0.013643008 0.060992262 0.027892847 0.014125
+033 0.018621974 0.055176232 0.016928658 0.015921152 0.5123
+8998]
]
dbl_epsilon=2.22044604925031e-16
r=0.0134242898847
r=0.0100183936697395
r=0.0100000985667805
r=0.0100000005281864
r=0.0100000000028304
r=0.0100000000000152
r=0.0100000000000001
[
[0.98990983 0.94029293 0.94371323 0.9430323 0.95885789 0.94865982 0.
+95112207 0.95174101 0.93703741 0.93760395]
[0.93236752 0.99061516 0.94178966 0.92991394 0.93407374 0.95334023 0.
+94093615 0.9284966 0.9456056 0.92144301]
[ 0.9313034 0.93730535 0.98766594 0.9539097 0.93060051 0.94203362 0.
+94056219 0.94050444 0.9496323 0.91744771]
[0.93692071 0.9317407 0.96035631 0.99062164 0.91989381 0.94711763 0
+.9641522 0.93721181 0.94537517 0.90647761]
[0.92877194 0.91245623 0.91341257 0.89684268 0.99071482 0.9000624 0
+.8898454 0.90581207 0.90782484 0.91424115]
[0.93464055 0.94723572 0.94047958 0.93920794 0.91548638 0.98544458 0.
+95216579 0.9220568 0.94583278 0.91659698]
[0.94782559 0.9456455 0.94979208 0.96707799 0.91548638 0.96309835 0.
+98908433 0.93549698 0.94417056 0.91864259]
[0.95064132 0.93530715 0.95193569 0.94223533 0.93407374 0.93480595 0
+.9376659 0.99570789 0.93254784 0.90729619]
[0.91529606 0.9315169 0.9399592 0.92946403 0.91548638 0.93774535 0
+.9254713 0.9119644 0.98971618 0.90268235]
[0.93418835 0.9258903 0.92628621 0.90906687 0.94041803 0.92695645 0.
+91847941 0.90503673 0.92075759 0.98845386]
]
-
Are you posting in the right place? Check out Where do I post X? to know for sure.
-
Posts may use any of the Perl Monks Approved HTML tags. Currently these include the following:
<code> <a> <b> <big>
<blockquote> <br /> <dd>
<dl> <dt> <em> <font>
<h1> <h2> <h3> <h4>
<h5> <h6> <hr /> <i>
<li> <nbsp> <ol> <p>
<small> <strike> <strong>
<sub> <sup> <table>
<td> <th> <tr> <tt>
<u> <ul>
-
Snippets of code should be wrapped in
<code> tags not
<pre> tags. In fact, <pre>
tags should generally be avoided. If they must
be used, extreme care should be
taken to ensure that their contents do not
have long lines (<70 chars), in order to prevent
horizontal scrolling (and possible janitor
intervention).
-
Want more info? How to link
or How to display code and escape characters
are good places to start.
