Beefy Boxes and Bandwidth Generously Provided by pair Networks
more useful options
 
PerlMonks  

Re: translating code and creating matrix

by spx2 (Deacon)
on Aug 17, 2013 at 16:29 UTC ( #1049856=note: print w/replies, xml ) Need Help??


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:

#!/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] ]

Replies are listed 'Best First'.
Re^2: translating code and creating matrix
by BrowserUk (Patriarch) on Aug 17, 2013 at 17:10 UTC

    You seem to have got similar results to those I got with pure perl.


    With the rise and rise of 'Social' network sites: 'Computers are making people easier to use everyday'
    Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
    "Science is about questioning the status quo. Questioning authority".
    In the absence of evidence, opinion is indistinguishable from prejudice.

Log In?
Username:
Password:

What's my password?
Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: note [id://1049856]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others exploiting the Monastery: (4)
As of 2022-12-10 07:51 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    No recent polls found

    Notices?