in reply to translating code and creating matrix
@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:
And, here are the results I got from this code above ^ :#!/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');
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] ]
|
|---|
| Replies are listed 'Best First'. | |
|---|---|
|
Re^2: translating code and creating matrix
by BrowserUk (Patriarch) on Aug 17, 2013 at 17:10 UTC |
In Section
Seekers of Perl Wisdom