#! /usr/bin/perl -w # # Global variables # my @ALPHABET = (A, R, N, D, C, Q, E, G, H, I, L, K, M, F, P, S, T, W, Y, V); my $MinSpeed = 1; my $MaxSpeed = 10; my $lapd = 'lapd'; # Application name # # Global options # my $verbose = 0; my $use_ml = 0; my $write_std_dev = 1; my $speed = 5; # Decides the spacing of samples in Octave my $model = 'wag'; my $modelfile = undef; my $octave = 'octave -q'; my $prior = "flat_prior"; my $debug = 0; # If 1, write octave commands to output instead of to an octave process. my $usage = "Usage: $lapd [] The infile should be in either FASTA, STOCKHOLM, or PHYLIP format. Output is a matrix of expected distances and, if possible, estimates of standard deviation. Options: -indels Remove gap columns. A gap can is denoted by '-'. -ml Compute a Maximum Likelihood estimate instead. This option implies -sd. -sd Do not output a matrix with standard deviations after the distance matrix. -id Output percent identity. -jc Use a simplistic Jukes-Cantor model. -jck Use -jc, but use Kimura's correction. -jcss Like -jck, but using Storm-Sonhammer's correction. -wag Default. Use the WAG matrix (see Wheelan and Goldman, 2001). -jtt Use the JTT matrix (see Jones, Taylor, Thornton, 1992). -day Use the Dayhoff matrix (Dayhoff et al, 1978). -arve Use the Arvestad matrix. -mv Use the Muller-Vingron matrix (2000). -f Read matrix and equilibrium distribution from file. -pfam Use a normal distribution as distance prior, estimated from Pfam 7.2. -s \"Speed\". High speed results in low precision. Default is $speed. Valid range is [$MinSpeed, $MaxSpeed]. -v Verbose. Show progress info on STDERR. -octave Point to the Octave binary to run. '$octave' by default. -d Debug option. Output Octave commands to STDOUT. -u, -h This help text. "; # # Reading sequence files # sub read_input { my ($filename) = @_; open(F, "<$filename") or die "Could not open $filename for reading!\n"; my $format = 'fasta'; my $seq_no = 0; my %seqs = (); my @seqnames = (); # The seqnames are put here to preserve order my $seq = undef; my $line=0; my %len_db = (); # Store seqlengths for input checking my $phylip_line = 0; # In the interleaved phylip format, we won't find sequence names on every line! my %line2seqname = (); # Use this to map lineno (on interleaved line) to a name while () { $line++; chomp; if ($format ne 'phylip' and length($_) == 0) { next; } if ($line == 1) { if (m/\d+\s+\d+/) { $format = 'phylip'; next; } elsif (m/^\#\s*STOCKHOLM\s+1.0\s*$/i) { $format = 'stockholm'; next; } } if ($format eq 'phylip') { if (m/^(\S+)\s+(\S+.+)$/) { my $name = $1; push @seqnames, $name; $phylip_line++; $line2seqname{$phylip_line} = $name; my $newblock = $2; $newblock =~ s/\ //g; # Remove blanks from Phylip's format! push @{$len_db{length($newblock)}}, $name; $seqs{$name} = $newblock; } elsif (m/^\s{10,15}(\S+.+)$/) { $phylip_line++; my $newblock = $1; $newblock =~ s/\ //g; # Remove blanks from Phylip's format! $seqs{$line2seqname{$phylip_line}} .= $newblock; } elsif (m/^\s*$/) { # Empty line $phylip_line = 0; } } elsif ($format eq 'stockholm') { if (m/^\#/) { next; } elsif (m/\/\//) { last; } elsif (m/^(\S+)\s+(\S+)$/) { if (exists $seqs{$1}) { $seqs{$1} .= $2; } else { push @seqnames, $1; $seqs{$1} = $2; } } } else { # fasta if (m/^(\#|>\s*)(\S+)/) { if (defined $seq) { $seqs{$curname} = $seq; push @seqnames, $curname; push @{$len_db{length($seq)}}, $curname; } $seq_no++; $seq = ''; $curname = $2; } elsif (m/^([-.a-zA-Z]+)\s*$/) { $seq .= $1; } else { print STDERR "Garbage in input on line $line: '$_'\n"; exit 1; } } } if (defined $seq) { $seqs{$curname} = $seq; push @seqnames, $curname; push @{$len_db{length($seq)}}, $curname; } close(F); if ($format eq 'stockholm') { foreach my $sid (keys %seqs) { push @{$len_db{length($seqs{$sid})}}, $sid; } } if (scalar(keys(%len_db)) > 1) { foreach my $len (keys %len_db) { print "Length = $len:\t", join(', ', @{$len_db{$len}}), "\n"; } die "Sequence lengths are not agreeing. Make sure all sequences are aligned."; } return (\%seqs, \@seqnames); } # # Remove gap columns # sub remove_gaps { my ($seqs) = @_; # First find the actual gap columns my %gap_db = (); for my $name (keys %$seqs) { my $pos = -1; my $s = $seqs->{$name}; while (($pos = index($s, '-', $pos)) > -1) { $gap_db{$pos} = undef; $pos++; } } # Remove those bad columns my @gap_positions = sort {$b <=> $a} keys %gap_db; foreach my $name (keys %$seqs) { foreach my $pos (@gap_positions) { substr($seqs->{$name}, $pos, 1, ''); } } } # # Compute a count matrix # sub count_replacements { my ($s1, $s2) = @_; my %M = (); for (my $i = 1; $i <= length($s1); $i++) { my $c1 = uc(substr($s1, $i, 1)); my $c2 = uc(substr($s2, $i, 1)); $M{$c1}{$c2}++; } return \%M; } # # Compute all count matrices # sub count_for_all_pairs { my ($seqs, $names) = @_; # The following is a matrix with references to count matrices my %counts = (); my $nseqs = scalar(@$names); for (my $i = 1; $i <= $nseqs; $i++) { for (my $j = 1; $j <= $nseqs; $j++) { my $n1 = $names->[$i]; my $n2 = $names->[$j]; $counts{$n1}{$n2} = count_replacements($seqs{$n1}, $seqs{$n2}); } } return \%counts; } # # Output identity based distances # sub output_id_dists { my ($seqs, $names) = @_; my %dists = (); for (my $i = 0; $i < @$names; $i++) { for (my $j = $i+1; $j < @$names; $j++) { $dists{$i}{$j} = 100 * count_id_dist($seqs->{$names->[$i]}, $seqs->{$names->[$j]}); $dists{$j}{$i} = $dists{$i}{$j}; } } print_matrix(\%dists, $names); } # # Output Jukes-Cantor distances # sub output_jc_dists { my ($seqs, $names) = @_; my %dists = (); for (my $i = 0; $i < @$names; $i++) { for (my $j = $i+1; $j < @$names; $j++) { $dists{$i}{$j} = -(1900/20)*log(1 - (20/19) * (1 - count_id_dist($seqs->{$names->[$i]}, $seqs->{$names->[$j]}))); $dists{$j}{$i} = $dists{$i}{$j}; } } print_matrix(\%dists, $names); } # # Output Kimura-corrected Jukes-Cantor distances # sub output_kimura_dists { my ($seqs, $names) = @_; my %dists = (); for (my $i = 0; $i < @$names; $i++) { for (my $j = $i+1; $j < @$names; $j++) { my $diff = 1 - count_id_dist($seqs->{$names->[$i]}, $seqs->{$names->[$j]}); my $adjusted = $diff + 0.2 * $diff * $diff; if ($adjusted > 0.854) { $adjusted = 0.854; } $dists{$i}{$j} = -100 * log(1 - $adjusted); $dists{$j}{$i} = $dists{$i}{$j}; } } print_matrix(\%dists, $names); } # # Output Storm-Sonnhammer corrected Jukes-Cantor distances # sub output_stormsonnhammer_dists { my ($seqs, $names) = @_; my %dists = (); for (my $i = 0; $i < @$names; $i++) { for (my $j = $i+1; $j < @$names; $j++) { my $diff = 1 - count_id_dist($seqs->{$names->[$i]}, $seqs->{$names->[$j]}); my $adjusted; if ($diff > 0.916) { $adjusted = 1000.0; } else { my $diff2 = $diff * $diff; my $diff3 = $diff * $diff2; $adjusted = -100 * log(1 - 0.95844 * $diff - 0.69957 * $diff2 + 2.4955 * $diff3 - 4.6353 * $diff2 * $diff2 + 2.8076 * $diff2 * $diff3); } $dists{$i}{$j} = $adjusted; $dists{$j}{$i} = $dists{$i}{$j}; } } print_matrix(\%dists, $names); } # # Compute the identity based distance # sub count_id_dist { my ($s1, $s2) = @_; my $id = 0; my $len = length($s1); for (my $i=0; $i<$len; $i++) { if (uc(substr($s1, $i, 1)) eq uc(substr($s2, $i, 1))) { $id++; } } return $id / $len; } # # Write a distance matrix to stdout # sub print_matrix { my ($dists, $names) = @_; my @lengthsorted = sort {length($b) <=> length($a)} @$names; my $maxlen = length($lengthsorted[0]); if ($maxlen < 10) { $maxlen = 10; } my $nameformat = '%-' . ($maxlen + 1) . 's'; print " ", scalar(@$names), "\n"; for (my $i = 0; $i < @$names; $i++) { printf($nameformat, $names->[$i]); for (my $j = 0; $j < @$names; $j++) { if ($i == $j) { printf("%12.4f", 0.0); } else { printf("%12.4f", $dists->{$i}{$j}); } } print "\n"; } } # # Write an octave input file # sub write_octave_commands { my ($fh, $seqs, $names) = @_; print $fh 'seqnames = {\'', join('\', \'', @$names), "'}; \nn_seqnames = length(seqnames);\n"; my $nseqs = scalar(@$names); for (my $i = 1; $i <= $nseqs; $i++) { for (my $j = $i+1; $j <= $nseqs; $j++) { my $n1 = $names->[$i-1]; my $n2 = $names->[$j-1]; my $counts = count_replacements($seqs->{$n1}, $seqs->{$n2}); # Print a matrix to work with print $fh "N = [\n"; foreach my $aminoacid1 (@ALPHABET) { foreach my $aminoacid2 (@ALPHABET) { if (exists $counts->{$aminoacid1}{$aminoacid2}) { print $fh $counts->{$aminoacid1}{$aminoacid2}, " "; } else { print $fh "0 "; } } print $fh "\n"; } print $fh "]; \n"; if ($use_ml) { write_ml_commands($fh, $i, $j); } else { write_expectation_commands($fh, $i, $j); } if ($verbose) { print $fh "fprintf(stderr, '.');" } } if ($verbose) { print $fh "fdisp(stderr, '');" } } } # # Sometimes, due to bad alignments or partial sequences, # there is not any shared residues between the sequences. # That has to be detected. # # # sub write_expectation_commands { my ($fh, $i, $j) = @_; print $fh "if (sum(sum(N)) == 0) fdisp(stderr, strcat('Warning: No data between ', seqnames{$i}, ' and ', seqnames{$j}, '.')); d = -1; sd = 0; else [d sd]= expected_distance(N, $prior); endif D($i, $j) = d; D($j, $i) = d; SD($i, $j) = sd; SD($j, $i) = sd; "; } # # # sub write_ml_commands { my ($fh, $i, $j) = @_; print $fh "if (sum(sum(N)) == 0) fdisp(stderr, strcat('Warning: No data between ', seqnames{$i}, ' and ', seqnames{$j}, '.')); d = -1; else d = likelihood_root(N, Q, kimura_distance(N)); endif D($i, $j) = d; D($j, $i) = d;\n "; } # # The Arvestad matrix # sub write_arvestad_Q { my ($fh) = @_; # Rate matrix estimate from Pfam v. 6.6 # Amino acid order is: A, R, N, D, C, Q, E, G, H, I, L, K, M, F, P, S, T, W, Y, V print $fh "Q = [ -0.0111450766711027639 0.0003841380303231597 0.0002175854080946476 0.0003521542195514120 0.0002669638364065408 0.0004365356451808754 0.0007311339616417847 0.0011407998945816499 0.0001271297777259876 0.0003982631605211135 0.0007840683734416032 0.0005266167003455223 0.0002676001252262387 0.0002438571205018998 0.0004925689069618634 0.0021135125224150957 0.0009510534710324932 0.0000536187152705582 0.0001566146976050900 0.0015008621042752290 0.0005873087336090157 -0.0099043170025695064 0.0004713717055640566 0.0003726864413974183 0.0000635652689815325 0.0009409682203956624 0.0006874793030423610 0.0003059821838639237 0.0004092383224601026 0.0001724151049514483 0.0005395675825044375 0.0030981217158649720 0.0001325995386427465 0.0001111832049748879 0.0003052306643452101 0.0006715192145275873 0.0004525885058565650 0.0000831983895503830 0.0002325874455445390 0.0002667054564926617 0.0003979530572472622 0.0005638796671962664 -0.0116013913347076496 0.0019488223361607837 0.0000853739262302032 0.0007971900261696733 0.0007517133188951303 0.0009350412032847155 0.0006128671005269068 0.0002435941539464559 0.0002546035738584158 0.0012127688254852569 0.0000885852481654407 0.0000876065956389980 0.0002630446828045371 0.0018560126013446897 0.0010213419765600033 0.0000140118083078400 0.0002578736510877240 0.0002091075817973467 0.0005237837993513853 0.0003625631022572382 0.0015848542729410798 -0.0095197424649322085 0.0000268311735162764 0.0003539734103315741 0.0027305474006842673 0.0006085359309004409 0.0002646888876366220 0.0000399433337831520 0.0002094883148883043 0.0006476877129334123 0.0000548768079801702 0.0000337120058983304 0.0003768438537459551 0.0009765285436626009 0.0004655955037664709 0.0000200060636948392 0.0001250814855030196 0.0001142008614570688 0.0013929101554434849 0.0002169258986254321 0.0002435532033285662 0.0000941220056158557 -0.0076795875978423606 0.0001005006156276208 0.0000701327778869987 0.0004151748498443514 0.0001169436304578363 0.0004511237100851912 0.0008441968719933680 0.0000934087505809919 0.0002242769192071091 0.0003482568562577715 0.0000763453450698211 0.0010485019118782963 0.0005756568452491751 0.0001236361070450887 0.0001976546961295271 0.0010462664475158745 0.0009454191500801991 0.0013329089543557182 0.0009439836260707480 0.0005154139375353786 0.0000417160073549957 -0.0139608299098068819 0.0027460491642639188 0.0003712699928218229 0.0007551156087912495 0.0001886154334174483 0.0007340335199835685 0.0020720614722589334 0.0003579901348114239 0.0000846257646984294 0.0004753488843336183 0.0009968341186142981 0.0008071176631092741 0.0000543510962053408 0.0002417994381041070 0.0002961759429964076 0.0009243341655148496 0.0005684765079230585 0.0005196156404336309 0.0023209329982487896 0.0000169934837042621 0.0016030079394066519 -0.0107579276605016590 0.0003662106349794556 0.0001917520314020558 0.0001487928294402615 0.0002846174634794862 0.0014695522273540720 0.0000868603305267518 0.0000975920625983819 0.0004276630049951934 0.0007519422606337631 0.0005356453748384806 0.0000334064307736285 0.0001213132605964490 0.0002892190136524380 0.0013338422137138969 0.0002339979063304322 0.0005977556033035967 0.0004783679190314189 0.0000930369260457341 0.0002004380421178836 0.0003386833834019861 -0.0059632541623125448 0.0001068612378319279 0.0000927356214196799 0.0002014286088441646 0.0003471522225063544 0.0000754728753639967 0.0001192001879935420 0.0002198750584771417 0.0009873208897961815 0.0002623500163521789 0.0000437588068069786 0.0000652613602950172 0.0001657152826804325 0.0004541438146163714 0.0009561879309396714 0.0011970446385847673 0.0006357152937858569 0.0000800667163792632 0.0012455319720673028 0.0005418187639007619 0.0003264910877730027 -0.0103512600206307516 0.0002304670607338436 0.0005704545804872165 0.0006250931152755674 0.0001967219792022638 0.0002558892388844414 0.0003720222998335946 0.0006908125787466656 0.0004509387317857185 0.0000997852876886987 0.0011157717815434328 0.0003063031484023090 0.0004940151128905100 0.0001398834246100340 0.0001652094487202673 0.0000333115808777767 0.0001072494521513114 0.0001080295685918424 0.0001459889729700440 0.0000983833712273278 0.0000800263200008674 -0.0131856578498646716 0.0038203662349283935 0.0001974887121429628 0.0007481140593275463 0.0007025751704596634 0.0001433339883284000 0.0001784768199097677 0.0005793322595108167 0.0000855405219838035 0.0002586407968356530 0.0050996920343976846 0.0006134514909510973 0.0002761168977158212 0.0001089152535986633 0.0001101962785925457 0.0001265899879774949 0.0002651780052233342 0.0001761391051674444 0.0001347883812849587 0.0001249399264162423 0.0024096900750003990 -0.0094466285967599265 0.0001816824168615657 0.0011518143140169450 0.0011610383685376794 0.0001775190985247151 0.0001999190894275248 0.0003439737777426158 0.0001248804189304690 0.0003155799496833878 0.0014442157611070218 0.0007268117010250155 0.0027967046523305122 0.0009151730639182973 0.0006009990370173341 0.0000247083745395097 0.0013204599551779668 0.0016042803349567603 0.0004097814051514000 0.0002415048643570266 0.0002197351130172714 0.0003204895633350191 -0.0120917218871214940 0.0002205830889570587 0.0001207110808328812 0.0004397594887497955 0.0008493512702217245 0.0007295889944065780 0.0000373382223934679 0.0001895973336484162 0.0003241443430854602 0.0009089507980163709 0.0002945892685282428 0.0001645179883702517 0.0001253209676962498 0.0001460050827153308 0.0005614620907297761 0.0002333691757772440 0.0002192552359716124 0.0001870513436555280 0.0020485745841802264 0.0050004633009716245 0.0005428739397374917 -0.0143529060915510154 0.0008289667992336709 0.0001098956636173810 0.0004728663190137326 0.0008151688298988788 0.0001108642603318546 0.0003321896200960298 0.0012505208230095160 0.0004237084990701422 0.0001263548071379257 0.0000832274253363663 0.0000393819451917157 0.0001159740275115031 0.0000678937281334163 0.0001341263736654249 0.0001771387962063447 0.0001244622480218373 0.0009841344449232122 0.0025784102334880455 0.0001519678047408064 0.0004240477940901323 -0.0092038134294945093 0.0001062524801689431 0.0003139669611325489 0.0002502828711792517 0.0003412200648477311 0.0019533491544404869 0.0008079137702086748 0.0008643212483194419 0.0003503137661235420 0.0002523687895781103 0.0004445805586778850 0.0000256755696743688 0.0003851376636465385 0.0005935780474720747 0.0003299811853701441 0.0001827389197426735 0.0002027623342046197 0.0003981319120941308 0.0005591085479892998 0.0000567720646749489 0.0001073039076765120 -0.0066821806498501392 0.0009264327111617841 0.0005012267356855777 0.0000445257816289211 0.0001057152567177687 0.0003515056494117976 0.0025368333797801770 0.0005271890827906119 0.0012180523673400961 0.0007880484066499036 0.0002412046493822713 0.0005524657116374846 0.0007139035863536153 0.0010135623238183707 0.0002321139483726297 0.0001727025850097738 0.0003067009302277807 0.0007386644123458261 0.0001670980711420213 0.0002168899024743423 0.0006337132728607953 -0.0131893858199581757 0.0024999379178738844 0.0000519591165400011 0.0002204967118124049 0.0003578494435461841 0.0013606781218838215 0.0004235208867558297 0.0007989499427219233 0.0004478577352694917 0.0001578495943774908 0.0005331907852889951 0.0006061718690695004 0.0003210233163915254 0.0001806017233670178 0.0006682024754233376 0.0006289983599420414 0.0007563128134218651 0.0003433554687452488 0.0002060866507000561 0.0004086735209184456 0.0029798374400965472 -0.0124218575775928257 0.0000511221044367747 0.0001965485696131229 0.0013528761991697894 0.0003133116270765855 0.0003179771218497115 0.0000447664047252082 0.0000785963884164109 0.0001384633261420731 0.0001466438089602022 0.0001544038541871147 0.0002186908405844039 0.0001632226806068736 0.0004029599695157462 0.0009326703381057989 0.0001580834710996163 0.0001907206020261155 0.0011475270588326817 0.0001482734744012995 0.0002529500374620360 0.0002087941682308945 -0.0064804909780848277 0.0011351886976033252 0.0003272471082587312 0.0003477193677211420 0.0003377563167860683 0.0003130410241095721 0.0001867111398830221 0.0000841071335883829 0.0002478831286682606 0.0002130456113690201 0.0001239245930171680 0.0006934667184226061 0.0004629384949811873 0.0008955286357058814 0.0003050010311044384 0.0002171342635187191 0.0024959995897825839 0.0001337598699518758 0.0004078602340921948 0.0003050111990761698 0.0004313246865152644 -0.0086988245674337415 0.0004966115291401862 0.0016737910800153659 0.0001945419228083062 0.0001275052767260266 0.0000856269888433666 0.0002236310408746527 0.0001525125398303491 0.0002551263110076338 0.0001580622040038085 0.0000956237567954554 0.0045849468762958756 0.0020585728809517667 0.0002619217844311924 0.0004105795570076962 0.0005185536949115199 0.0002234007210389635 0.0003324862463942317 0.0010545521135097275 0.0000624562876348475 0.0002494484959738458 -0.0127233397790546300 ]; # Residue distribution to accompany Q: eq = [0.0784206791118385 0.0512922139357810 0.0428774076662180 0.0527243742237138 0.0150300328183314 0.0362097824533222 0.0620295385941313 0.0670711284618227 0.0219525251335432 0.0632208746217339 0.1002315182575809 0.0568202730012320 0.0230874801985448 0.0451334845488915 0.0446912398236554 0.0653346367326373 0.0548125657865228 0.0134205554509911 0.0353211011097199 0.0703185880697882]; "; } sub write_jtt_Q { my ($fh) = @_; print $fh "eq=[ 0.076748 0.051691 0.042645 0.051544 0.019803 0.040752 0.061830 0.073152 0.022944 0.053761 0.091904 0.058676 0.023826 0.040126 0.050901 0.068765 0.058565 0.014261 0.032102 0.066005]; Q = 0.01 * [ -1.24783051471057305 0.02978918457290290 0.02288113515059247 0.04148382800570307 0.01101881019688046 0.02308019485608830 0.06450661211114973 0.13010528043232419 0.00615529094729973 0.01923027641897842 0.02739500325529957 0.02040535862975805 0.01278381817558955 0.00598043556658116 0.09811688057662078 0.25827057827209382 0.27640598656646942 0.00127528771791378 0.00350865335550063 0.19543789990282701 0.04422936947633344 -1.02596500952821201 0.01906761262549373 0.00819433639618826 0.02223438486156236 0.12552386676118199 0.01781611191641278 0.09957778446496321 0.07477538632275223 0.01175183558937570 0.03470033745671279 0.37662461928067720 0.01041644443936926 0.00199347852219372 0.03742602661170071 0.06900880530550654 0.03724206976895588 0.01785402805079294 0.00637936973727387 0.01114914194076530 0.04117906813313802 0.02311229837552811 -1.28056765885934554 0.27041310107421263 0.00668999190524885 0.03482275013374726 0.03563222383282556 0.05887445650848190 0.08913773186645160 0.02510619421366627 0.01095800130211983 0.15333169484646764 0.00710212120866086 0.00398695704438744 0.00758635674561501 0.34367751553138409 0.13500250291246507 0.00113358908259003 0.02232779408045855 0.01049331006189675 0.06176860219970703 0.00821770608907666 0.22372665480579307 -1.04490308529733400 0.00196764467801437 0.01984086926225135 0.47120544275477944 0.09448986847040305 0.02553305874435442 0.00587591779468785 0.00639216742623657 0.01515826641067741 0.00355106060433043 0.00159478281775498 0.00758635674561501 0.04031207438638501 0.02211247892531756 0.00056679454129501 0.01467255039572990 0.02033078824492495 0.04270421880473572 0.05803754925410392 0.01440664065037304 0.00512146024761766 -0.56570991728323761 0.00364424129306657 0.00307174343386427 0.04288386338272138 0.01573018797643263 0.00908096386451759 0.02100283582906301 0.00408107172595161 0.00733885858228289 0.03109826494622206 0.00708059962924068 0.15236597607057387 0.02444010828587730 0.01629534306223165 0.06666441375451196 0.04066157648984991 0.04346679414053458 0.15921805547586029 0.03644032635094357 0.02509515521332655 0.00177088021021293 -1.04414963129574923 0.19843462582763199 0.01889797369408061 0.13610032205696060 0.00480756910474460 0.06574800781271897 0.17023899199683862 0.01017970706574723 0.00159478281775498 0.08294416708539076 0.03621254139793908 0.02967727434713672 0.00255057543582756 0.00765524368472865 0.01311663757737094 0.08007041025887948 0.01489459228645145 0.02457603405063636 0.39281600099227476 0.00098382233900718 0.13078777085116705 -0.94434324711808060 0.08649457190752280 0.00592731720851085 0.00641009213965947 0.00821850097658987 0.10552485462817734 0.00426127272519652 0.00199347852219372 0.00910362809473801 0.02049766494222967 0.01862103488447794 0.00141698635323753 0.00223277940804586 0.02951243454908461 0.13650098510799455 0.07036410838771891 0.03432170272588871 0.06657898321902962 0.01160910360028477 0.01052780817997010 0.07310749372596968 -0.64783144534511261 0.00524339599214421 0.00320504606982974 0.00547900065105991 0.01574127665724193 0.00331432323070840 0.00199347852219372 0.01213817079298401 0.13733435511293876 0.01920294222461788 0.00779342494280644 0.00255174789490955 0.03082409830682171 0.02058953406656901 0.16846297482607153 0.16567636747928993 0.05736035477331783 0.01357674827829914 0.24173467244008273 0.01597306585609422 0.01671743826784054 -1.13569504763143492 0.00854678951954596 0.05113733940989253 0.02623546109540321 0.00781233332952695 0.01594782817754977 0.05816206838304840 0.04987765135942552 0.02676773764643704 0.00113358908259003 0.18276894297289642 0.00721415066755402 0.02745271208875868 0.01129934587248041 0.01991506207551567 0.00563360627237943 0.00334499595262443 0.00364424129306657 0.00737218424127425 0.00436107085248014 0.00364757982062206 -1.27360212913373583 0.20911519151545341 0.01224321517785483 0.11339720196495172 0.03548391769504824 0.00505757116374334 0.02733021992297289 0.14256729833428422 0.00127528771791378 0.01020699157963819 0.63025443559267358 0.02287726007396557 0.01951705196155707 0.00508469670013166 0.00358502217333236 0.00452558275943305 0.02915393034453259 0.00552913818095569 0.00436107085248014 0.01276652937217721 0.12232592499850160 -0.66813928121309374 0.00816214345190322 0.09185410096534712 0.09887653470080858 0.05158722587018206 0.04031207438638501 0.01454768350349839 0.00736832903683518 0.00765524368472865 0.11804973819633845 0.02669013675295983 0.33178988334647019 0.11143960267788557 0.01331579664380592 0.00137735127461006 0.11823538417504885 0.11119711230588668 0.01962481883616063 0.01025881824549955 0.01121766124440408 0.01278433485247313 -0.90073435689846570 0.01538792928543186 0.00159478281775498 0.01062089944386101 0.03211300840949315 0.05993645603441337 0.00141698635323753 0.00255174789490955 0.00918164630415966 0.04117906813313802 0.02259869174496082 0.01271174175032915 0.00768219037142649 0.00609969850184454 0.01741137506687363 0.01105827636191138 0.01017583198912033 0.00752313338003300 0.25586951124140728 0.35430870876854115 0.03789566602669352 -1.18205128769077272 0.01714391529086600 0.00809211386198934 0.01981440944415534 0.13151105887162545 0.00340076724777008 0.00574143276354649 0.21183369687454068 0.01143863003698278 0.00256803315283646 0.00423724725010972 0.00204858409904707 0.01534762848851207 0.00161966279691848 0.00307174343386427 0.00363422571040012 0.00911894955155515 0.04754151670247442 0.22646536024380981 0.00233204098625806 0.01017970706574723 -0.63718230752234473 0.00859787097836368 0.06285950582283764 0.00698288808167923 0.00751002767215893 0.17096710895893974 0.04066157648984991 0.14793961514497733 0.03800689066197956 0.00635587087516457 0.00768219037142649 0.00275470254922012 0.06640617467365757 0.01105827636191138 0.01744428340992056 0.02621697996072106 0.00534174344971623 0.09314301106801855 0.01224321517785483 0.00378779797795246 0.00677782697545865 -0.72767547164022972 0.19472781695118185 0.06866506613651241 0.00085019181194252 0.00318968486863694 0.01508413321397658 0.28825347693196612 0.05187426968729642 0.21313353668051879 0.03021661546094421 0.04387847631972041 0.02146053205916983 0.01843046060318564 0.14609587355808470 0.01664208293158815 0.02136697379886491 0.05387683973542249 0.02740148158853224 0.00686538383503883 0.03668000480836447 0.14414077816668516 -1.44786266571788103 0.27756980124674929 0.00495945223633137 0.02009501467241270 0.02492161139700479 0.36222328450445485 0.03287082435630664 0.09830413620254544 0.01946154894094712 0.00826410764766035 0.02065070066071058 0.01965915797673135 0.02398588968864077 0.01048679198428842 0.13087271451804755 0.02282916937941631 0.06005005539614512 0.05350264643857847 0.00478434845326493 0.05967933973217140 0.32591287258145168 -1.33538948074266073 0.00170038362388504 0.00669833822413757 0.07345317043327726 0.00686317802218967 0.06471443545147870 0.00338979780008777 0.00204858409904707 0.02262791379716523 0.00728848258613315 0.00614348686772855 0.03997648281440128 0.00182378991031103 0.00480756910474460 0.04748467230918593 0.00583010246564516 0.00568169696692869 0.02113087233525345 0.00303454269824600 0.02391394243260128 0.00698288808167923 -0.31278499728186271 0.02264676256732225 0.01639579697171367 0.00838832869378738 0.01027213261134583 0.02966073075076802 0.02355871713904125 0.04112377377050029 0.00971797678151086 0.00430044080740998 0.00581476113664019 0.13062895232602753 0.01709357903909193 0.02191600260423966 0.00466408197251613 0.00426127272519652 0.21370089757916697 0.00505757116374334 0.04304509637868230 0.01222005414293865 0.01006060310798650 -0.60597828279249011 0.01049331006189675 0.22724745006805797 0.00873131271964395 0.00677959560017555 0.01587652676761475 0.01219939700368908 0.00809831398459239 0.02764569090477845 0.03416172167776110 0.00250771112667767 0.51334154551772937 0.16437001953179742 0.00816214345190322 0.07646617167991525 0.02471913367520214 0.01163241367660968 0.02596370892682425 0.06517362209567280 0.00354246588309384 0.00510349578981910 -1.24172244008155808 ]; "; }sub write_day_Q { my ($fh) = @_; print $fh "eq=[0.087127 0.040904 0.040432 0.046872 0.033474 0.038255 0.049530 0.088612 0.033618 0.036886 0.085357 0.080482 0.014753 0.039772 0.050680 0.069577 0.058542 0.010494 0.029916 0.064718]; Q = 0.01 *[ -1.331753078285893 0.010980908021398 0.039396714951244 0.055924671401760 0.011981710512689 0.033852201936162 0.097508444443872 0.211452337525720 0.007687912270517 0.023838758197173 0.034796191267431 0.020805627470855 0.010561401397724 0.007118010851865 0.125975277824057 0.282942194097256 0.215948491636422 0.000000000000000 0.007138771418566 0.133843453061184 0.023389731399871 -0.864729027242428 0.012864233453467 0.000000000000000 0.007654981716440 0.093569007598831 0.000492466891131 0.007929462657215 0.080221693257568 0.023472008071062 0.012730313878328 0.371300428710638 0.013201751747155 0.005536230662562 0.051901814463511 0.106535691664982 0.015133856556730 0.020972288686868 0.002379590472855 0.015443475353214 0.084896062118051 0.013014409506842 -1.787482595563020 0.421765230154942 0.000000000000000 0.039177267409266 0.072885099887338 0.122466145483646 0.178827524553330 0.028239759710497 0.028855378124211 0.254468828297377 0.000146686130524 0.005536230662562 0.021163846674442 0.342436151780298 0.133294351980433 0.002399814128348 0.028257636865156 0.009652172095758 0.103954361777205 0.000000000000000 0.363816602355876 -1.413318779610123 0.000000000000000 0.050968483813997 0.567814325473657 0.110131425794646 0.028746106750629 0.008802003026648 0.000000000000000 0.056815367324257 0.000000000000000 0.000000000000000 0.006550714446851 0.065720069533592 0.038416712797854 0.000000000000000 0.000000000000000 0.011582606514910 0.031186308533162 0.009354106833043 0.000000000000000 0.000000000000000 -0.266079185167489 0.000000000000000 0.000000000000000 0.009691565469929 0.009359197546716 0.016137005548855 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.009574121114628 0.111378223104299 0.009313142496449 0.000000000000000 0.028555085674263 0.031530428846144 0.077099484984761 0.100048273083847 0.041406751428348 0.062449216398632 0.000000000000000 -1.237758906610106 0.352606294049556 0.024669439378001 0.202559775475360 0.006601502269986 0.061954194207865 0.122433115501568 0.016722218879729 0.000000000000000 0.077096870028323 0.038740251514539 0.030849784519489 0.000000000000000 0.000000000000000 0.022521734890103 0.171524696932389 0.000406700297089 0.059497079722287 0.537342884385246 0.000000000000000 0.272339062767328 -1.359867886826981 0.071365163914931 0.014373053375314 0.022371757692731 0.009335563510774 0.066417964618498 0.004400583915718 0.000000000000000 0.025698956676108 0.054651426243724 0.019790427804955 0.000000000000000 0.006543873800352 0.023808691169537 0.207908723554410 0.003660302673799 0.055879014063499 0.058254866043500 0.003661078212211 0.010650130946208 0.039889818181584 -0.650440008835634 0.003342570552399 0.000000000000000 0.005940813143220 0.021605843912041 0.002493664218907 0.005931675709888 0.017132637784072 0.161878908114323 0.017462142180843 0.000000000000000 0.000000000000000 0.034747819544730 0.019924586007298 0.097608071301314 0.215073903050159 0.040079347837928 0.009319108176536 0.230499262621510 0.021176076318619 0.008810514063572 -0.875478506647675 0.002567250882772 0.037342254043096 0.020805627470855 0.000000000000000 0.018981362271640 0.047366704461845 0.024212657196587 0.012805570932618 0.002817173107191 0.037775998756577 0.028313038147558 0.056308612629319 0.026028819013684 0.030954561747406 0.011184934280352 0.014644312848842 0.006846512751134 0.030040480358971 0.000000000000000 0.002339799386679 -1.277529629931566 0.218112711115359 0.036809956294589 0.049286539856044 0.077507229275864 0.006046813335555 0.016602964934802 0.111757709957394 0.000000000000000 0.011005605936955 0.572052066208617 0.035517740273878 0.006100504456332 0.013668248044309 0.000000000000000 0.000000000000000 0.027766412824043 0.005417135802437 0.006167359844500 0.014707310430554 0.094254782410360 -0.530222348062144 0.014403895941361 0.077303590786116 0.062084872430157 0.016124835561479 0.011760433495485 0.019208356398927 0.004799628256696 0.008328566654993 0.112608674450515 0.022523445051728 0.188708937849208 0.127838319943833 0.033088763912708 0.000000000000000 0.058195358384639 0.040874751963845 0.023788387971644 0.008690683436237 0.016870505801076 0.015276376653994 -0.744004184233594 0.035644729717317 0.000000000000000 0.016628736672775 0.066411859739209 0.079161711219820 0.000000000000000 0.003866834518390 0.006434781397172 0.062372617066323 0.036603026737993 0.000402007295421 0.000000000000000 0.000000000000000 0.043361247423848 0.014774006733920 0.014977873908072 0.000000000000000 0.123228042373078 0.447258360925269 0.194452595208373 -1.251820819940516 0.036380944353977 0.008566318892036 0.042890992748239 0.060535426226922 0.000000000000000 0.000000000000000 0.166017360047045 0.015593154266581 0.005693804159243 0.005628102135892 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.013215771095358 0.016044338651514 0.071883024717629 0.133243951926503 0.000000000000000 0.013495124008202 -0.543000287953162 0.005542912224258 0.031822349458371 0.007566928278365 0.007929820598020 0.207619268756619 0.007721737676607 0.216571587035844 0.041890130600147 0.016884306407676 0.006058506068524 0.006323680548364 0.058195358384639 0.025115811447664 0.029955747816144 0.031420163192548 0.004401001513324 0.027158002940434 0.026407142559162 0.002493664218907 0.004349895520584 -0.743002119006685 0.169488600376107 0.045401569670191 0.000000000000000 0.000000000000000 0.030886950706427 0.354311116390641 0.062631845751677 0.198993611233324 0.044273698193060 0.053584872015081 0.021300261892417 0.038904884399323 0.206166029087577 0.011698996933395 0.008802003026648 0.014427689062105 0.076820778353925 0.009094540092484 0.018190472176989 0.123455772267576 -1.600038928746133 0.320139273315451 0.007825480853309 0.010113259509635 0.019304344191517 0.321392235161193 0.010574207724309 0.092059670651376 0.030758569270968 0.005325204672306 0.020159176433894 0.016743874298443 0.026431542190715 0.007353655215277 0.070416024213187 0.028006690532322 0.108829436001394 0.015255357574490 0.005140785615236 0.039304286681106 0.380484613089220 -1.291754247243531 0.000000000000000 0.012492849982490 0.101026067935605 0.000000000000000 0.081746759714851 0.009246167794680 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.009024940491476 0.000000000000000 0.039039629226874 0.000000000000000 0.000000000000000 0.030053823596764 0.000000000000000 0.051884265421257 0.000000000000000 -0.239139963601423 0.018144377355521 0.000000000000000 0.020790872355441 0.003253602376710 0.038190693064981 0.000000000000000 0.031951228033837 0.000000000000000 0.010834271604875 0.000000000000000 0.042450646015463 0.013569754666083 0.023763252572880 0.010402813735427 0.000000000000000 0.276020643033436 0.000000000000000 0.023520866990970 0.024446999053180 0.006364724427358 -0.543577755842724 0.018017387912082 0.180187560413822 0.009760807130131 0.006030109431313 0.008388700710264 0.016308439308938 0.013312663682760 0.018221274971835 0.047576775943287 0.014707310430554 0.326040862112102 0.148520328580497 0.008002164411867 0.037845021675176 0.004745340567910 0.024187253342219 0.020753706168503 0.091385210746410 0.000000000000000 0.008328566654993 -0.984302096282584 ]; "; } sub write_mullervingron_Q { my ($fh) = @_; print $fh " Q=[ -1.1411028e-02 3.1690849e-04 2.9013185e-04 3.7559234e-04 1.8174787e-04 4.6178663e-04 7.5261309e-04 1.1149029e-03 1.0736606e-04 2.7734892e-04 5.4006580e-04 4.5046348e-04 2.1266883e-04 1.4604906e-04 6.6618241e-04 2.3011367e-03 1.5203925e-03 3.5828436e-05 9.0782088e-05 1.5690613e-03 4.8515338e-04 -9.4310795e-03 3.7157420e-04 2.0862332e-04 1.0698086e-04 1.1457039e-03 4.3513405e-04 4.3792289e-04 5.4580859e-04 1.3386759e-04 4.1678950e-04 3.0053505e-03 1.3873245e-04 8.9463876e-05 3.1842073e-04 6.9178341e-04 3.8597097e-04 9.7296830e-05 1.4808794e-04 2.6841445e-04 4.7133402e-04 3.9815154e-04 -1.1726112e-02 2.1483740e-03 7.0898738e-05 6.2388333e-04 6.4327205e-04 8.0278519e-04 6.9453676e-04 2.7580234e-04 2.0203855e-04 1.3322474e-03 1.0330597e-04 1.0495604e-04 1.6829203e-04 2.2168092e-03 1.1232738e-03 1.1705931e-05 1.7992664e-04 1.5451808e-04 5.4026550e-04 1.9030448e-04 1.8482646e-03 -9.4852140e-03 1.6121436e-05 3.4126882e-04 3.1283569e-03 6.4080303e-04 2.3959964e-04 6.5748579e-05 1.7748815e-04 4.3565237e-04 4.4435446e-05 3.4406955e-05 2.5650008e-04 7.9593887e-04 4.3488289e-04 1.6205275e-05 1.2191013e-04 1.5706075e-04 9.6915446e-04 3.6338430e-04 2.3324679e-04 6.4306296e-05 -7.1248064e-03 9.6669423e-05 1.5314963e-04 3.6723966e-04 1.6753784e-04 2.8723233e-04 7.0770531e-04 1.0405299e-04 1.3671763e-04 3.5325069e-04 1.0478952e-04 1.2909862e-03 5.0895544e-04 1.0361573e-04 3.4721745e-04 7.6559465e-04 8.6606635e-04 1.4039880e-03 7.1184667e-04 4.5705561e-04 3.3255947e-05 -1.2203904e-02 2.1800277e-03 2.9299738e-04 8.8586834e-04 1.5295037e-04 6.6006051e-04 1.8222396e-03 2.1344208e-04 7.5142567e-05 5.5952649e-04 8.3535948e-04 6.2439829e-04 5.8321649e-05 1.5335622e-04 2.1800027e-04 9.1689095e-04 3.4572356e-04 4.7530871e-04 2.6435549e-03 3.1482952e-05 1.4257760e-03 -1.0221110e-02 4.6196480e-04 1.8186272e-04 1.2312026e-04 2.7838809e-04 1.3028962e-03 8.2632147e-05 4.7231590e-05 3.6310249e-04 6.5812121e-04 4.5850011e-04 2.7102651e-05 9.4673080e-05 3.0277800e-04 1.3061551e-03 3.3375765e-04 5.6785003e-04 5.2136252e-04 8.1226747e-05 1.8220058e-04 4.4657482e-04 -6.0392495e-03 9.5333435e-05 9.0083396e-05 1.7968531e-04 2.7063990e-04 6.6044973e-05 7.7277925e-05 1.8681855e-04 1.0999271e-03 2.0786522e-04 6.0238783e-05 5.5960607e-05 2.1024679e-04 3.7545458e-04 1.2582925e-03 1.4462256e-03 5.8480694e-04 1.1268229e-04 1.6686434e-03 5.2742624e-04 2.8513421e-04 -1.1074327e-02 1.7795292e-04 5.3479989e-04 5.4264283e-04 1.5322229e-04 2.8807895e-04 4.4767125e-04 7.0764859e-04 4.6757253e-04 4.1052918e-05 1.2719774e-03 1.8304165e-04 3.3708652e-04 1.0690461e-04 2.1617044e-04 5.8498739e-05 6.4008045e-05 9.7381506e-05 1.2822283e-04 9.2782545e-05 6.3004980e-05 -1.1933615e-02 2.9746570e-03 1.8559743e-04 7.6681895e-04 4.7882056e-04 1.1253651e-04 2.1606170e-04 7.6813720e-04 3.5385501e-05 1.3969211e-04 5.0918473e-03 4.1849541e-04 2.1334500e-04 9.7025736e-05 9.7735423e-05 1.0220511e-04 2.7738163e-04 1.8085194e-04 1.1794951e-04 1.1968124e-04 1.8095751e-03 -7.7950168e-03 1.6772097e-04 9.0331644e-04 9.1848270e-04 2.4874759e-04 3.2551321e-04 2.9282907e-04 9.8179482e-05 2.0430110e-04 1.2016803e-03 5.8392973e-04 2.5594471e-03 1.0300280e-03 3.9854915e-04 2.4145555e-05 1.2580055e-03 1.3948555e-03 2.9672954e-04 2.0185458e-04 1.9032080e-04 2.8364564e-04 -1.0705235e-02 1.9196038e-04 8.2046736e-05 3.3396747e-04 7.1072265e-04 7.3293322e-04 3.7994193e-05 1.3008230e-04 2.6401692e-04 7.4695104e-04 3.2119228e-04 2.2454052e-04 1.1461209e-04 8.9963267e-05 3.9900399e-04 2.4417605e-04 2.0081696e-04 1.5532179e-04 2.0123174e-03 3.9794234e-03 5.1634628e-04 -1.2629915e-02 4.9590302e-04 1.2031547e-04 3.8847933e-04 9.6921938e-04 9.6522863e-05 1.8190775e-04 1.3729018e-03 2.6727787e-04 1.0681172e-04 1.1984499e-04 4.6871750e-05 1.2570651e-04 7.4283604e-05 7.6922727e-05 1.2374145e-04 1.5082430e-04 6.9087735e-04 2.1614810e-03 1.1796694e-04 2.5601803e-04 -7.5146356e-03 9.0231324e-05 4.0287852e-04 1.9927068e-04 2.1534349e-04 1.7084738e-03 5.7980949e-04 1.0673233e-03 3.3232223e-04 1.7082816e-04 2.9675883e-04 3.1823736e-05 4.7739850e-04 4.9254580e-04 2.5952155e-04 2.0333378e-04 1.4993039e-04 4.9197306e-04 4.2337327e-04 5.5478483e-05 7.9307785e-05 -6.6735580e-03 1.2657325e-03 5.6199322e-04 2.8477195e-05 8.0421627e-05 2.0501465e-04 2.5249187e-03 4.8812317e-04 1.4299581e-03 6.0693586e-04 2.7068680e-04 4.8149286e-04 5.8778867e-04 1.0231000e-03 2.1940381e-04 1.8320430e-04 4.3966992e-04 5.9575276e-04 1.1796059e-04 2.3032997e-04 8.6644039e-04 -1.3135752e-02 2.5542055e-03 4.3499655e-05 1.5660170e-04 3.1567974e-04 2.0815748e-03 3.4009930e-04 9.1199691e-04 4.1213833e-04 1.2847834e-04 4.4749228e-04 5.0887725e-04 2.3650120e-04 1.7884214e-04 7.8496685e-04 5.0838853e-04 7.7179204e-04 3.7794714e-04 1.4427915e-04 4.7286045e-04 3.1818541e-03 -1.2822830e-02 2.9422169e-05 1.4161443e-04 1.1637050e-03 2.1277266e-04 3.8894994e-04 4.2421934e-05 6.7936489e-05 1.2129952e-04 1.8478401e-04 1.3530304e-04 3.1243606e-04 6.9499070e-05 1.7246686e-04 7.5284909e-04 1.7522655e-04 1.6175407e-04 7.0542811e-04 1.0286795e-04 2.4621928e-04 1.2911224e-04 -4.9206227e-03 7.4845330e-04 1.9084255e-04 2.0803073e-04 2.2481166e-04 2.5528818e-04 2.0482983e-04 1.6215334e-04 1.8817460e-04 1.8956899e-04 1.1165455e-04 8.6663198e-04 2.6387448e-04 6.2024109e-04 2.3385307e-04 1.2059541e-04 2.1845517e-03 1.1522980e-04 3.4453968e-04 2.4610391e-04 2.9246012e-04 -7.1715591e-03 3.3896594e-04 1.8374874e-03 1.9671115e-04 1.0886470e-04 1.2507276e-04 1.6720854e-04 1.3154049e-04 2.8935373e-04 2.0693484e-04 5.8813907e-05 4.4486864e-03 1.7673446e-03 2.3608946e-04 4.5802956e-04 3.5822488e-04 1.4415477e-04 3.3561074e-04 9.8963626e-04 3.6869577e-05 1.6522412e-04 -1.2061858e-02 ]; eq = [ 7.7076463e-02 5.0081937e-02 4.6237740e-02 5.3792986e-02 1.4453339e-02 4.0892361e-02 6.3357935e-02 6.5567235e-02 2.1880269e-02 5.9196969e-02 9.7646129e-02 5.9207941e-02 2.2069588e-02 4.1350852e-02 4.7687160e-02 7.0729517e-02 5.6775916e-02 1.2701980e-02 3.2374605e-02 6.6919080e-02]; "; } sub write_wag_Q { my ($fh) = @_; print $fh " Q = 0.01 * [ -1.064444707752500108 0.024253680012000001 0.019929652411200000 0.042156214809800002 0.019829882912000000 0.033371078203800003 0.091898529864999995 0.117944490095999999 0.007743598260200000 0.009370174110000001 0.034303854235000000 0.056214349178999999 0.017425584439200000 0.008089684358599998 0.065832507504999999 0.234330242140999984 0.129414648096999985 0.001627520024700000 0.008491734536999999 0.142217282555999996 0.047781437430899999 -0.928350780929100061 0.024835293932400002 0.008402971410399999 0.010198206189800000 0.111488147550000000 0.025496972347300001 0.048674413647000002 0.052213352795000001 0.009062124214000001 0.042903719238999993 0.331941090612000000 0.013323503537400002 0.003947378880899999 0.031095523055900000 0.085103118000999983 0.033826234045100004 0.016744036728000001 0.013458271348600000 0.017854985964400001 0.044167061559199992 0.027937434311999996 -1.383913476725120217 0.309721806842000014 0.005121509796800000 0.056694964283999998 0.054993273962199996 0.093704896007999985 0.096657307876999998 0.026861601975999998 0.011338897351999999 0.186830763486000001 0.003865844696700000 0.003695692210990000 0.008927511311100001 0.276280123716999948 0.123859441762000019 0.001034586454530000 0.038307781200000002 0.013912977917600001 0.064017844844200006 0.006477251487999999 0.212232770148000044 -0.942973292519679984 0.000584927870220000 0.022653267702299997 0.358464938023999979 0.072061426051200000 0.022737624558799997 0.001911353642000000 0.007310928382299999 0.029764733852999999 0.002023483135800000 0.001795938059760000 0.019402822190400001 0.074506504503999996 0.022871586798199998 0.001866815085300000 0.011489194956199999 0.010799881226000001 0.088970318415999991 0.023225614651999998 0.010368697886400000 0.001728175599990000 -0.464363904347719902 0.003629393712990000 0.001239673632800000 0.025531162513199999 0.006082709623600000 0.008245762910000001 0.033128997982999994 0.004592219169540000 0.007615453301400001 0.015296664837999999 0.005006666192400000 0.097857567113999983 0.031298538896800003 0.010315697313000000 0.019183274008599999 0.071047316584000000 0.078709936684199996 0.133477005999999981 0.060339961416000003 0.035184447913299999 0.001907956249620000 -1.314688531726940335 0.317551411783000004 0.027477423093600000 0.104910689642999988 0.005521101322000000 0.074957777200999998 0.241595194139999997 0.030136742202000005 0.003840146193520000 0.042713996173200000 0.071524881772999996 0.052344503685600001 0.003103570908300000 0.008032288081999999 0.021359497263600001 0.137118971514999993 0.019310611604000001 0.037025401501200005 0.352205574615999994 0.000412260145600000 0.200883241107000010 -1.178538388149710237 0.047263462140600003 0.013926451782500000 0.006174326070000001 0.013298858966999999 0.160308574698000000 0.006145768834800001 0.003118129931410000 0.031226680100500004 0.049005878908099994 0.050199114115500001 0.002252213346300000 0.006924431282600000 0.041738437483600004 0.122727478487999989 0.025708889379999999 0.043997465063999999 0.049377325838400002 0.005921200257200000 0.012122182861200001 0.032961024531300002 -0.474138393494500099 0.006093410533000000 0.001475794546600000 0.005284930673300000 0.023171279758800001 0.003395420070000000 0.001918943198900000 0.011146518267000001 0.093280508577999993 0.013778681079100002 0.004847803739700000 0.003654548216800000 0.013274988413200000 0.027457059416600000 0.093974759800000002 0.154649002326000012 0.053090505487599998 0.004807101581600000 0.157714501490999998 0.033095024472499998 0.020763831438000001 -0.945524792749799947 0.006697516540000001 0.043058119557999999 0.055232250355200005 0.007881840680699999 0.026109518334899998 0.031860178693800000 0.051454994525099995 0.028877737998900001 0.003777291377100000 0.136632497248000001 0.008391061424800001 0.016748205046500000 0.008221840588000000 0.021664752698400000 0.002249687608700000 0.003284932553000000 0.004183954967700000 0.007396413565500000 0.002535025635180000 0.003376161347000000 -1.174983186929160039 0.273366152729999978 0.020086845595200003 0.083031965142000008 0.040717445092999995 0.004573051667280000 0.022206797975999996 0.088966278632000004 0.003056759189700000 0.014821160613999998 0.554495756280000029 0.034470540828499999 0.021883589211999996 0.005141350603200001 0.004837692591970000 0.007419736538599999 0.031934678940899998 0.008956340090700000 0.005103643371660000 0.012202505960599999 0.153684232020000006 -0.691758700294730078 0.015975776072999999 0.094666495854000013 0.081290001922999996 0.019030310556400001 0.023965531328100000 0.019928090099400001 0.009571068743100000 0.014060931055599999 0.127636184504000022 0.078507833793499987 0.235312640239999982 0.117737663694000017 0.027373376460500001 0.001429431734420000 0.143052276689999980 0.150049162926999990 0.031099375904400001 0.021754411321599998 0.015694841712000002 0.022203558994999997 -1.071523983755319875 0.018220904545200003 0.003414136268400000 0.025485287337599999 0.067232846626999992 0.084623394646000008 0.001978133179500000 0.004700780988799999 0.021653926690400001 0.077401682138399985 0.030039999463999999 0.007748339957400001 0.005918657305400000 0.007539348359599999 0.056754463806000001 0.018295752803600001 0.014494138380000000 0.009873690013299999 0.206342056359999998 0.418460210180000014 0.057951832293600000 -1.259723418632500147 0.045758173096999998 0.007840546159900000 0.034335238399499993 0.092502574724000003 0.007418894945400000 0.015112772425400000 0.145935047820000025 0.018234653182599998 0.004516408091999999 0.003758918791740000 0.002665740341040000 0.007684890555999999 0.003669901134480000 0.004710544986710000 0.004156845625800001 0.016597916712299998 0.051348273019999995 0.182346690529999994 0.005510372709600000 0.023220499701000003 -0.679999371059869961 0.007388177916400000 0.037951976664900001 0.010488266168100000 0.022005248076000003 0.227669563575999989 0.046074483275200001 0.124618565544999990 0.029878490307999998 0.007625599241400001 0.024186209678400003 0.002112350551200000 0.034280980153199997 0.039616780709499998 0.020277640926000002 0.017009022197399998 0.004843149220800000 0.035849495395999999 0.034543479225599999 0.003341378088300000 0.006204599663599999 -0.577018522993099969 0.112151837711999991 0.048528525376800004 0.002005466389500000 0.007620849813199999 0.022324102797199999 0.292004459040999975 0.053830082679999995 0.155350266162000011 0.061138656376000002 0.027178817748000000 0.037788440246999996 0.040927982907099994 0.111708930275999999 0.018083290889700000 0.015481979040000001 0.029719604450999998 0.059989719918000006 0.009632481043500001 0.020981165598900002 0.073828693968000000 -1.326554610767000097 0.267114820854000012 0.007534500037800001 0.027760548480600000 0.016500171048400000 0.183747304968999964 0.024378648435999999 0.079353827364000012 0.021384268456599999 0.009904592475199999 0.031510065376799998 0.047768830858499998 0.018801003749400001 0.011563505309099999 0.070671182560000004 0.028157755998000002 0.086032427628000008 0.029568433524000001 0.006606558905700000 0.036399237530400003 0.304350756557999980 -1.100482689685899684 0.001594878417600000 0.010270012781600000 0.098419398788000004 0.009800474210699999 0.051179890239999998 0.002811180652980000 0.007402571491699999 0.013845044145999999 0.007923610109700000 0.009089527207300000 0.028054441319399999 0.006414902009700000 0.010298201078000000 0.057355623581000000 0.008529242643000000 0.010057659406200001 0.058786971516000000 0.006379604955500000 0.036409443981800002 0.006764111972800001 -0.444675698936180031 0.087670143938000003 0.025903054476400001 0.020854367506500000 0.016776769075999998 0.042451088400000003 0.018580216566100000 0.010500218797400001 0.008363355651000000 0.011397136246699999 0.008625219487200001 0.094633174671999998 0.020363959220000002 0.034364459161999997 0.008266179350399999 0.008355678279900001 0.248050243531999975 0.009886934702599999 0.054710100674700002 0.017763725579600001 0.035754572001000004 -0.692010371093099930 0.022312972188000001 0.173776433678999975 0.011074304227999999 0.007671138392400001 0.008689965308500001 0.019349118692000001 0.011065478696100001 0.034181074255900001 0.015588649794600000 0.002891639805400000 0.379067125800000015 0.155205511059999995 0.018945643412400001 0.040145332815000004 0.024976584354799998 0.014410205269699999 0.016179526528099997 0.084699660521000003 0.005256161897100000 0.011101848966000000 -1.034275403476000221 ]; eq = [0.0866279 0.043972 0.0390894 0.0570451 0.0193078 0.0367281 0.0580589 0.0832518 0.0244313 0.048466 0.086209 0.0620286 0.0195027 0.0384319 0.0457631 0.0695179 0.0610127 0.0143859 0.0352742 0.0708956 ]; "; } # # Basic Octave initialization # sub write_octave_prologue { my ($fh, $m) = @_; print $fh "PS1=\"> \";split_long_rows=0;more off;\n"; print $fh "StepSize = $speed;\n"; if ($m eq 'jtt') { write_jtt_Q($fh); } elsif ($m eq 'wag') { write_wag_Q($fh); } elsif ($m eq 'day') { write_day_Q($fh); } elsif ($m eq 'mv') { write_mullervingron_Q($fh); } elsif ($m eq 'file') { open(MF, "<$modelfile") or die "Could not open '$modelfile' for reading!\n"; while () { print $fh $_; } close(MF); } else { write_arvestad_Q($fh); } while () { print $fh $_; } } sub write_octave_ending { my ($fh) = @_; # First write the number of sequences print $fh " n = length(seqnames); disp(strcat(' ', int2str(n)));"; # Then write the distance matrix print $fh "write_matrix(seqnames, D);"; # And here comes the standard deviations if ($write_std_dev) { print $fh "\ndisp('');\n"; print $fh "write_matrix(seqnames, SD);\n"; } } ### Main program ### if (@ARGV < 1) { print STDERR $usage; exit 2; } my $remove_indel_columns = 0; while (substr($ARGV[0], 0, 1) eq '-') { my $option = shift @ARGV; if ($option eq '-indels') { $remove_indel_columns = 1; } elsif ($option eq '-s') { # Speed, which affects precision, $speed = shift @ARGV; if ($speed =~ m/^\d+$/) { if ($speed < 1 || $speed > 10) { print STDERR "$lapd error: Invalid range for '-s' option. Speed must be in [$MinSpeed, $MaxSpeed].\nComputation aborted.\n"; exit 5; } } else { print STDERR "$lapd error: The argument to '-s' must be an integer, found '$speed'.\n"; exit 6; } } elsif ($option eq '-sd') { $write_std_dev = 0; } elsif ($option eq '-ml') { $write_std_dev = 0; # Don't have standard deviations in ML estimation $use_ml = 1; } elsif ($option eq '-id') { $model = 'id'; } elsif ($option eq '-jc') { $model = 'jc'; } elsif ($option eq '-jck') { $model = 'jck'; } elsif ($option eq '-jcss') { $model = 'jcss'; } elsif ($option eq '-jtt') { $model = 'jtt'; } elsif ($option eq '-day') { $model = 'day'; } elsif ($option eq '-wag') { $model = 'wag'; } elsif ($option eq '-arve') { $model = 'arve'; } elsif ($option eq '-mv') { $model = 'mv'; } elsif ($option eq '-f') { $modelfile = shift @ARGV; if (!defined $modelfile) { print STDERR "$lapd error: Missing argument to option '-f'!\n\n$usage"; exit 7; } $model = 'file'; } elsif ($option eq '-octave') { $octave = shift @ARGV; } elsif ($option eq '-pfam') { $prior = 'norm_prior'; } elsif ($option eq '-v') { $verbose = 1; } elsif ($option eq '-d') { $debug = 1; } elsif ($option eq '-u' || $option eq '-h') { print $usage; exit 0; } else { print STDERR "$lapd: Unknown option: '$option'\n\n$usage"; exit 4; } } if (@ARGV == 0) { print STDERR "$lapd: No infile given on command line!\n\n $usage"; exit 3; } my $infile = shift @ARGV; # my $outfile = undef; # if (@ARGV == 1) { # $outfile = shift @ARGV; # } my ($seqs, $names) = read_input($infile); if (@$names == 0) { print STDERR "No sequences read from '$infile'!\n"; exit 8; } if ($remove_indel_columns) { remove_gaps($seqs); } if ($model eq 'id') { output_id_dists($seqs, $names); } elsif ($model eq 'jc') { output_jc_dists($seqs, $names); } elsif ($model eq 'jck') { output_kimura_dists($seqs, $names); } elsif ($model eq 'jcss') { output_stormsonnhammer_dists($seqs, $names); } else { my $fh; if ($debug) { $fh = STDOUT; } else { open($fh, "|$octave") or die "Could not call '$octave'!\n" } write_octave_prologue($fh, $model); write_octave_commands($fh, $seqs, $names); write_octave_ending($fh); } exit 0; ########################################################### __DATA__ # # How to sample from the divergence space # global MaxDistance = 400; global DSamples = 1:StepSize:MaxDistance; # These are the distance samples we use global DSamples2 = DSamples .^2; global DSamples3 = DSamples .^3; global DSamples4 = DSamples .^4; global DS_prev = shift(DSamples, 1); DS_prev(1) = 0; global DS2_prev = shift(DSamples2, 1); DS2_prev(1) = 0; global DS3_prev = shift(DSamples3, 1); DS3_prev(1) = 0; global DS4_prev = shift(DSamples4, 1); DS4_prev(1) = 0; global SampleGap = DSamples - shift(DSamples, 1); SampleGap(1) = DSamples(1); SampleGap = log(SampleGap ./ 2); # # This describes the prior distribution for, for example, pairwise distances in Pfam # global NormMean = 115; global NormVar = 2500; # # Compute a list (actually a cell array) of P(t) matrices # global prePt = {}; for d = DSamples global prePt; prePt(d) = log(diag(eq) * expm(Q*d)); endfor # # Some precomputations in case we want to use the Norm prior # global norm_prior = normal_pdf(DSamples, NormMean, NormVar); global flat_prior = ones(size(DSamples)); # # Compute the posterior probability of observing a set of replacements # # The integration code demands that the samples are uniformly distributed. # Numerical integration using simple linear interpolation. # function posterior_vec = posterior(identical, N, prior) global DSamples; global SampleGap; global prePt; l=[]; # middle datapoints for n = 1 : columns(DSamples) d = DSamples(n); p = loglikelihood(N, prePt{d}) + prior(n); # log-prob! l(1,n) = p; endfor # Integrate for total posterior probability if (identical) # Special treatment at first point ptot = log(1 + exp(p)) - log(2); else ptot = p - log(2); # Half interval endif for n = 2 : columns(DSamples) term = logprob_add(l(n), l(n-1)) + SampleGap(n); ptot = logprob_add(ptot, term); endfor # setup return value posterior_vec = exp(l .- ptot); endfunction function p = loglikelihood(N, Pt) p = sum(sum(N .* Pt)); endfunction # # Addition of log-probabilities # function logratio = logprob_add(p, q) if (p < q) logratio = q + log(1 + exp(p - q)); else logratio = p + log(1 + exp(q - p)); endif endfunction # # Compute expected distance with uniform prior # Returns mean and standard deviation. # function [d sigma] = expected_distance(N, prior) global DSamples; global DSamples2; global DSamples3; global DSamples4; global DS_prev; global DS2_prev; global DS3_prev; global DS4_prev; # Make sure there is data to compare! if (sum(sum(N)) == 0) d = -1; sigma = 0; return; endif identical = 0; if (sum(sum(diag(N))) == sum(sum(N))) identical = 1; endif pv = posterior(identical, N, prior); # clg; # plot(DSamples, pv, '+-'); # hold on; # Compute slant $k_{i+1}=(y_{i+1}-y_i) / (x_{i+1} - x_i)$ pv_prev = shift(pv,1); pv_prev(1) = identical; # First point gets messed up, so we fix it k = (pv - pv_prev) ./ (DSamples - shift(DSamples,1)); # plot(DSamples, k, '-@'); # Integrate for expected d: d = 0; d += sum((DSamples3 - DS3_prev) .* k / 3); d += sum((DSamples2 - DS2_prev) .* (pv_prev - DS_prev .* k) / 2); moment2 = 0; moment2 += sum((DSamples4 - DS4_prev) .* k / 4); moment2 += sum((DSamples3 - DS3_prev) .* (pv_prev - DS_prev .* k) / 3); sigma = sqrt(moment2 - d ^ 2); endfunction # Kimura distance, not corrected # Result is returned in PAMs! function D = kimura_distance(N) tot = sum(sum(N)); diagonal = triu(tril(N)); p = sum(sum(N-diagonal)) / tot; adjusted = p + 0.2*p*p;; if (adjusted > 0.854) # Infinite distance! adjusted = 0.854; endif D = -100 * log(1 - adjusted); endfunction # # Derivative of likelihood # function l = likelihood_deriv (N, Q, t) Pt = expm(Q .* t); direction = Pt * Q; l = sum(sum(N .* direction ./ Pt)); endfunction # # Computing the distance with highest likelihood # using Newton's method # # This is almost straight from Octave's manual. # # function distance = likelihood_root (N, Q, t_start) t = t_start; # Initial guess # The derivative is not always trustworthy around t == 0. if (sum(sum(N)) == sum(diag(N))) distance = 0; return; endif if (t == 0) t = 1; endif # Go ahead and do Newton-Raphson. delta = tol = 0.001; maxit = 50; Ld = likelihood_deriv(N, Q, t); for i = 1:maxit if (abs (Ld) < tol) distance = t; return; else L_new = likelihood_deriv(N, Q, t + delta); deriv = (L_new - Ld) / delta; t = t - Ld / deriv; if (t < 1) distance = 1; return; endif if (t > 500) distance = 500; return; endif if (abs(Ld) < abs(L_new)) distance = t; return; endif Ld = L_new; endif endfor distance = t; endfunction function write_matrix(names, M) n = length(names); maxlen = 0; for i = 1 : n if (maxlen < length(names{i})) maxlen = length(names{i}); endif endfor if (maxlen < 10) maxlen = 10; endif nameformat = strcat('%-', int2str(maxlen + 1), 's'); for i = 1 : n printf(nameformat, names{i}); for j = 1 : n if (i != j) printf("%12.4f", M(i,j)); else printf("%12.4f", 0.0); endif endfor printf("\n"); endfor endfunction