forester/archive/RIO/others/phylip_mod/src/protdist.c

   1 /*Modified by Christian Zmasek. Use at your own risk.*/
   2
   3 #include "phylip.h"
   4 #include "seq.h"
   5
   6 /* version 3.6. (c) Copyright 1993-2004 by the University of Washington.
   7    Written by Joseph Felsenstein, Akiko Fuseki, Sean Lamont, and Andrew Keeffe.
   8    Permission is granted to copy and use this program provided no fee is
   9    charged for it and provided that this copyright notice is not removed. */
  10
  11 #define nmlngth         26    /*changed from to 10 to 26 by CZ 2006-07-28 */ /* number of characters in species name     */
  12 #define protepsilon .00001
  13 typedef long *steparray;
  14 typedef enum {
  15   universal, ciliate, mito, vertmito, flymito, yeastmito
  16 } codetype;
  17 typedef enum {
  18   chemical, hall, george
  19 } cattype;
  20
  21 typedef double matrix[20][20];
  22
  23 #ifndef OLDC
  24 /* function prototypes */
  25 void   protdist_uppercase(Char *);
  26 void   protdist_inputnumbers(void);
  27 void   getoptions(void);
  28 void   transition(void);
  29 void   doinit(void);
  30 void   printcategories(void);
  31 void   inputoptions(void);
  32 void   protdist_inputdata(void);
  33 void   doinput(void);
  34 void   code(void);
  35 void   protdist_cats(void);
  36 void   maketrans(void);
  37 void   givens(matrix, long, long, long, double, double, boolean);
  38 void   coeffs(double, double, double *, double *, double);
  39 void   tridiag(matrix, long, double);
  40 void   shiftqr(matrix, long, double);
  41 void   qreigen(matrix, long);
  42 void   pmbeigen(void);
  43 void   pameigen(void);
  44 void   jtteigen(void);
  45 void   predict(long, long, long);
  46 void   makedists(void);
  47 void   reallocchars(void);
  48 /* function prototypes */
  49 #endif
  50
  51 long chars, datasets, ith, ctgry, categs;
  52 /* spp = number of species
  53    chars = number of positions in actual sequences */
  54 double freqa, freqc, freqg, freqt, cvi, invarfrac, ttratio, xi, xv,
  55   ease, fracchange;
  56 boolean weights, justwts, progress, mulsets, gama, invar, basesequal,
  57   usepmb, usejtt, usepam, kimura, similarity, firstset;
  58 codetype whichcode;
  59 cattype whichcat;
  60 steptr oldweight;
  61 double rate[maxcategs];
  62 aas **gnode;
  63 aas trans[4][4][4];
  64 double pie[20];
  65 long cat[(long)ser - (long)ala + 1], numaa[(long)ser - (long)ala + 1];
  66 double eig[20];
  67 matrix prob, eigvecs;
  68 double **d;
  69 char infilename[100], outfilename[100], catfilename[100], weightfilename[100];
  70
  71 /* Local variables for makedists, propagated globally for c version: */
  72   double tt, p, dp, d2p, q, elambdat;
  73
  74
  75 /* this jtt matrix decomposition due to Elisabeth  Tillier */
  76 static double jtteigs[] =
  77 {0.0,        -0.007031123, -0.006484345, -0.006086499, -0.005514432,
  78 -0.00772664, -0.008643413, -0.010620756, -0.009965552, -0.011671808,
  79 -0.012222418,-0.004589201, -0.013103714, -0.014048038, -0.003170582,
  80 -0.00347935, -0.015311677, -0.016021194, -0.017991454, -0.018911888};
  81
  82 static double jttprobs[20][20] =
  83 {{0.076999996, 0.051000003, 0.043000004, 0.051999998, 0.019999996, 0.041,
  84   0.061999994, 0.073999997, 0.022999999, 0.052000004, 0.090999997, 0.058999988,
  85   0.024000007, 0.04, 0.050999992, 0.069, 0.059000006, 0.014000008, 0.032000004,
  86   0.066000005},
  87  {0.015604455, -0.068062363, 0.020106264, 0.070723273, 0.011702977, 0.009674053,
  88   0.074000798, -0.169750458, 0.005560808, -0.008208636, -0.012305869,
  89  -0.063730179, -0.005674643, -0.02116828, 0.104586169, 0.016480839, 0.016765139,
  90   0.005936994, 0.006046367, -0.0082877},
  91  {-0.049778281, -0.007118197, 0.003801272, 0.070749616, 0.047506147,
  92    0.006447017, 0.090522425, -0.053620432, -0.008508175, 0.037170603,
  93    0.051805545, 0.015413608, 0.019939916, -0.008431976, -0.143511376,
  94   -0.052486072, -0.032116542, -0.000860626, -0.02535993, 0.03843545},
  95  {-0.028906423, 0.092952047, -0.009615343, -0.067870117, 0.031970392,
  96    0.048338335, -0.054396304, -0.135916654, 0.017780083, 0.000129242,
  97    0.031267424, 0.116333586, 0.007499746, -0.032153596, 0.033517051,
  98   -0.013719269, -0.00347293, -0.003291821, -0.02158326, -0.008862168},
  99  {0.037181176, -0.023106564, -0.004482225, -0.029899635, 0.118139633,
 100  -0.032298569, -0.04683198, 0.05566988, -0.012622847, 0.002023096,
 101  -0.043921088, -0.04792557, -0.003452711, -0.037744513, 0.020822974,
 102   0.036580187, 0.02331425, -0.004807711, -0.017504496, 0.01086673},
 103  {0.044754061, -0.002503471, 0.019452517, -0.015611487, -0.02152807,
 104  -0.013131425, -0.03465365, -0.047928912, 0.020608851, 0.067843095,
 105  -0.122130014, 0.002521499, 0.013021646, -0.082891087, -0.061590119,
 106   0.016270856, 0.051468938, 0.002079063, 0.081019713, 0.082927944},
 107  {0.058917882, 0.007320741, 0.025278141, 0.000357541, -0.002831285,
 108  -0.032453034, -0.010177288, -0.069447924, -0.034467324, 0.011422358,
 109  -0.128478324, 0.04309667, -0.015319944, 0.113302422, -0.035052393,
 110   0.046885372, 0.06185183, 0.00175743, -0.06224497, 0.020282093},
 111  {-0.014562092, 0.022522921, -0.007094389, 0.03480089, -0.000326144,
 112  -0.124039037, 0.020577906, -0.005056454, -0.081841576, -0.004381786,
 113   0.030826152, 0.091261631, 0.008878828, -0.02829487, 0.042718836,
 114  -0.011180886, -0.012719227, -0.000753926, 0.048062375, -0.009399129},
 115  {0.033789571, -0.013512235, 0.088010984, 0.017580292, -0.006608005,
 116  -0.037836971, -0.061344686, -0.034268357, 0.018190209, -0.068484614,
 117   0.120024744, -0.00319321, -0.001349477, -0.03000546, -0.073063759,
 118   0.081912399, 0.0635245, 0.000197, -0.002481798, -0.09108114},
 119  {-0.113947615, 0.019230545, 0.088819683, 0.064832765, 0.001801467,
 120  -0.063829682, -0.072001633, 0.018429333, 0.057465965, 0.043901014,
 121  -0.048050874, -0.001705918, 0.022637173, 0.017404665, 0.043877902,
 122  -0.017089594, -0.058489485, 0.000127498, -0.029357194, 0.025943972},
 123  {0.01512923, 0.023603725, 0.006681954, 0.012360216, -0.000181447,
 124  -0.023011838, -0.008960024, -0.008533239, 0.012569835, 0.03216118,
 125   0.061986403, -0.001919083, -0.1400832, -0.010669741, -0.003919454,
 126  -0.003707024, -0.026806029, -0.000611603, -0.001402648, 0.065312824},
 127  {-0.036405351, 0.020816769, 0.011408213, 0.019787053, 0.038897829,
 128    0.017641789, 0.020858533, -0.006067252, 0.028617353, -0.064259496,
 129   -0.081676567, 0.024421823, -0.028751676, 0.07095096, -0.024199434,
 130   -0.007513119, -0.028108766, -0.01198095, 0.111761119, -0.076198809},
 131  {0.060831772, 0.144097327, -0.069151377, 0.023754576, -0.003322955,
 132  -0.071618574, 0.03353154, -0.02795295, 0.039519769, -0.023453968,
 133  -0.000630308, -0.098024591, 0.017672997, 0.003813378, -0.009266499,
 134  -0.011192111, 0.016013873, -0.002072968, -0.010022044, -0.012526904},
 135  {-0.050776604, 0.092833081, 0.044069596, 0.050523021, -0.002628417,
 136    0.076542572, -0.06388631, -0.00854892, -0.084725311, 0.017401063,
 137   -0.006262541, -0.094457679, -0.002818678, -0.0044122, -0.002883973,
 138    0.028729685, -0.004961596, -0.001498627, 0.017994575, -0.000232779},
 139  {-0.01894566, -0.007760205, -0.015160993, -0.027254587, 0.009800903,
 140   -0.013443561, -0.032896517, -0.022734138, -0.001983861, 0.00256111,
 141    0.024823166, -0.021256768, 0.001980052, 0.028136263, -0.012364384,
 142   -0.013782446, -0.013061091, 0.111173981, 0.021702122, 0.00046654},
 143  {-0.009444193, -0.042106824, -0.02535015, -0.055125574, 0.006369612,
 144   -0.02945416, -0.069922064, -0.067221068, -0.003004999, 0.053624311,
 145    0.128862984, -0.057245803, 0.025550508, 0.087741073, -0.001119043,
 146   -0.012036202, -0.000913488, -0.034864475, 0.050124813, 0.055534723},
 147  {0.145782464, -0.024348311, -0.031216873, 0.106174443, 0.00202862,
 148   0.02653866, -0.113657267, -0.00755018, 0.000307232, -0.051241158,
 149   0.001310685, 0.035275877, 0.013308898, 0.002957626, -0.002925034,
 150  -0.065362319, -0.071844582, 0.000475894, -0.000112419, 0.034097762},
 151  {0.079840455, 0.018769331, 0.078685899, -0.084329807, -0.00277264,
 152  -0.010099754, 0.059700608, -0.019209715, -0.010442992, -0.042100476,
 153  -0.006020556, -0.023061786, 0.017246106, -0.001572858, -0.006703785,
 154   0.056301316, -0.156787357, -0.000303638, 0.001498195, 0.051363455},
 155  {0.049628261, 0.016475144, 0.094141653, -0.04444633, 0.005206131,
 156  -0.001827555, 0.02195624, 0.013066683, -0.010415582, -0.022338403,
 157   0.007837197, -0.023397671, -0.002507095, 0.005177694, 0.017109561,
 158  -0.202340113, 0.069681441, 0.000120736, 0.002201146, 0.004670849},
 159  {0.089153689, 0.000233354, 0.010826822, -0.004273519, 0.001440618,
 160   0.000436077, 0.001182351, -0.002255508, -0.000700465, 0.150589876,
 161  -0.003911914, -0.00050154, -0.004564983, 0.00012701, -0.001486973,
 162  -0.018902754, -0.054748555, 0.000217377, -0.000319302, -0.162541651}};
 163
 164 /* PMB matrix decomposition courtesy of Elisabeth Tillier */
 165 static double pmbeigs[] =
 166 {0.0000001586972220,-1.8416770496147100, -1.6025046986139100,-1.5801012515121300,
 167 -1.4987794099715900,-1.3520794233801900,-1.3003469390479700,-1.2439503327631300,
 168 -1.1962574080244200,-1.1383730501367500,-1.1153278910708000,-0.4934843510654760,
 169 -0.5419014550215590,-0.9657997830826700,-0.6276075673757390,-0.6675927795018510,
 170 -0.6932641383465870,-0.8897872681859630,-0.8382698977371710,-0.8074694642446040};
 171 static double pmbprobs[20][20] =
 172 {{0.0771762457248147,0.0531913844998640,0.0393445076407294,0.0466756566755510,
 173 0.0286348361997465,0.0312327748383639,0.0505410248721427,0.0767106611472993,
 174 0.0258916271688597,0.0673140562194124,0.0965705469252199,0.0515979465932174,
 175 0.0250628079438675,0.0503492018628350,0.0399908189418273,0.0641898881894471,
 176 0.0517539616710987,0.0143507440546115,0.0357994592438322,0.0736218495862984},
 177 {0.0368263046116572,-0.0006728917107827,0.0008590805287740,-0.0002764255356960,
 178 0.0020152937187455,0.0055743720652960,0.0003213317669367,0.0000449190281568,
 179 -0.0004226254397134,0.1805040629634510,-0.0272246813586204,0.0005904606533477,
 180 -0.0183743200073889,-0.0009194625608688,0.0008173657533167,-0.0262629806302238,
 181 0.0265738757209787,0.0002176606241904,0.0021315644838566,-0.1823229927207580},
 182 {-0.0194800075560895,0.0012068088610652,-0.0008803318319596,-0.0016044273960017,
 183 -0.0002938633803197,-0.0535796754602196,0.0155163896648621,-0.0015006360762140,
 184 0.0021601372013703,0.0268513218744797,-0.1085292493742730,0.0149753083138452,
 185 0.1346457366717310,-0.0009371698759829,0.0013501708044116,0.0346352293103622,
 186 -0.0276963770242276,0.0003643142783940,0.0002074817333067,-0.0174108903914110},
 187 {0.0557839400850153,0.0023271577185437,0.0183481103396687,0.0023339480096311,
 188 0.0002013267015151,-0.0227406863569852,0.0098644845475047,0.0064721276774396,
 189 0.0001389408104210,-0.0473713878768274,-0.0086984445005797,0.0026913674934634,
 190 0.0283724052562196,0.0001063665179457,0.0027442574779383,-0.1875312134708470,
 191 0.1279864877057640,0.0005103347834563,0.0003155113168637,0.0081451082759554},
 192 {0.0037510125027265,0.0107095920636885,0.0147305410328404,-0.0112351252180332,
 193 -0.0001500408626446,-0.1523450933729730,0.0611532413339872,-0.0005496748939503,
 194 0.0048714378736644,-0.0003826320053999,0.0552010244407311,0.0482555671001955,
 195 -0.0461664995115847,-0.0021165008617978,-0.0004574454232187,0.0233755883688949,
 196 -0.0035484915422384,0.0009090698422851,0.0013840637687758,-0.0073895139302231},
 197 {-0.0111512564930024,0.1025460064723080,0.0396772456883791,-0.0298408501361294,
 198 -0.0001656742634733,-0.0079876311843289,0.0712644184507945,-0.0010780604625230,
 199 -0.0035880882043592,0.0021070399334252,0.0016716329894279,-0.1810123023850110,
 200 0.0015141703608724,-0.0032700852781804,0.0035503782441679,0.0118634302028026,
 201 0.0044561606458028,-0.0001576678495964,0.0023470722225751,-0.0027457045397157},
 202 {0.1474525743949170,-0.0054432538500293,0.0853848892349828,-0.0137787746207348,
 203 -0.0008274830358513,0.0042248844582553,0.0019556229305563,-0.0164191435175148,
 204 -0.0024501858854849,0.0120908948084233,-0.0381456105972653,0.0101271614855119,
 205 -0.0061945941321859,0.0178841099895867,-0.0014577779202600,-0.0752120602555032,
 206 -0.1426985695849920,0.0002862275078983,-0.0081191734261838,0.0313401149422531},
 207 {0.0542034611735289,-0.0078763926211829,0.0060433542506096,0.0033396210615510,
 208 0.0013965072374079,0.0067798903832256,-0.0135291136622509,-0.0089982442731848,
 209 -0.0056744537593887,-0.0766524225176246,0.1881210263933930,-0.0065875518675173,
 210 0.0416627569300375,-0.0953804133524747,-0.0012559228448735,0.0101622644292547,
 211 -0.0304742453119050,0.0011702318499737,0.0454733434783982,-0.1119239362388150},
 212 {0.1069409037912470,0.0805064400880297,-0.1127352030714600,0.1001181253523260,
 213 -0.0021480427488769,-0.0332884841459003,-0.0679837575848452,-0.0043812841356657,
 214 0.0153418716846395,-0.0079441315103188,-0.0121766182046363,-0.0381127991037620,
 215 -0.0036338726532673,0.0195324059593791,-0.0020165963699984,-0.0061222685010268,
 216 -0.0253761448771437,-0.0005246410999057,-0.0112205170502433,0.0052248485517237},
 217 {-0.0325247648326262,0.0238753651653669,0.0203684886605797,0.0295666232678825,
 218 -0.0003946714764213,-0.0157242718469554,-0.0511737848084862,0.0084725632040180,
 219 -0.0167068828528921,0.0686962159427527,-0.0659702890616198,-0.0014289912494271,
 220 -0.0167000964093416,-0.1276689083678200,0.0036575057830967,-0.0205958145531018,
 221 0.0000368919612829,0.0014413626622426,0.1064360941926030,0.0863372661517408},
 222 {-0.0463777468104402,0.0394712148670596,0.1118686750747160,0.0440711686389031,
 223 -0.0026076286506751,-0.0268454015202516,-0.1464943067133240,-0.0137514051835380,
 224 -0.0094395514284145,-0.0144124844774228,0.0249103379323744,-0.0071832157138676,
 225 0.0035592787728526,0.0415627419826693,0.0027040097365669,0.0337523666612066,
 226 0.0316121324137152,-0.0011350177559026,-0.0349998884574440,-0.0302651879823361},
 227 {0.0142360925194728,0.0413145623127025,0.0324976427846929,0.0580930922002398,
 228 -0.0586974207121084,0.0202001168873069,0.0492204086749069,0.1126593173463060,
 229 0.0116620013776662,-0.0780333711712066,-0.1109786767320410,0.0407775100936731,
 230 -0.0205013161312652,-0.0653458585025237,0.0347351829703865,0.0304448983224773,
 231 0.0068813748197884,-0.0189002309261882,-0.0334507528405279,-0.0668143558699485},
 232 {-0.0131548829657936,0.0044244322828034,-0.0050639951827271,-0.0038668197633889,
 233 -0.1536822386530220,0.0026336969165336,0.0021585651200470,-0.0459233839062969,
 234 0.0046854727140565,0.0393815434593599,0.0619554007991097,0.0027456299925622,
 235 0.0117574347936383,0.0373018612990383,0.0024818527553328,-0.0133956606027299,
 236 -0.0020457128424105,0.0154178819990401,0.0246524142683911,0.0275363065682921},
 237 {-0.1542307272455030,0.0364861558267547,-0.0090880407008181,0.0531673937889863,
 238 0.0157585615170580,0.0029986538457297,0.0180194047699875,0.0652152443589317,
 239 0.0266842840376180,0.0388457366405908,0.0856237634510719,0.0126955778952183,
 240 0.0099593861698250,-0.0013941794862563,0.0294065511237513,-0.1151906949298290,
 241 -0.0852991447389655,0.0028699120202636,-0.0332087026659522,0.0006811857297899},
 242 {0.0281300736924501,-0.0584072081898638,-0.0178386569847853,-0.0536470338171487,
 243 -0.0186881656029960,-0.0240008730656106,-0.0541064820498883,0.2217137098936020,
 244 -0.0260500001542033,0.0234505236798375,0.0311127151218573,-0.0494139126682672,
 245 0.0057093465049849,0.0124937286655911,-0.0298322975915689,0.0006520211333102,
 246 -0.0061018680727128,-0.0007081999479528,-0.0060523759094034,0.0215845995364623},
 247 {0.0295321046399105,-0.0088296411830544,-0.0065057049917325,-0.0053478115612781,
 248 -0.0100646496794634,-0.0015473619084872,0.0008539960632865,-0.0376381933046211,
 249 -0.0328135588935604,0.0672161874239480,0.0667626853916552,-0.0026511651464901,
 250 0.0140451514222062,-0.0544836996133137,0.0427485157912094,0.0097455780205802,
 251 0.0177309072915667,-0.0828759701187452,-0.0729504795471370,0.0670731961252313},
 252 {0.0082646581043963,-0.0319918630534466,-0.0188454445200422,-0.0374976353856606,
 253 0.0037131290686848,-0.0132507796987883,-0.0306958830735725,-0.0044119395527308,
 254 -0.0140786756619672,-0.0180512599925078,-0.0208243802903953,-0.0232202769398931,
 255 -0.0063135878270273,0.0110442171178168,0.1824538048228460,-0.0006644614422758,
 256 -0.0069909097436659,0.0255407650654681,0.0099119399501151,-0.0140911517070698},
 257 {0.0261344441524861,-0.0714454044548650,0.0159436926233439,0.0028462736216688,
 258 -0.0044572637889080,-0.0089474834434532,-0.0177570282144517,-0.0153693244094452,
 259 0.1160919467206400,0.0304911481385036,0.0047047513411774,-0.0456535116423972,
 260 0.0004491494948617,-0.0767108879444462,-0.0012688533741441,0.0192445965934123,
 261 0.0202321954782039,0.0281039933233607,-0.0590403018490048,0.0364080426546883},
 262 {0.0115826306265004,0.1340228176509380,-0.0236200652949049,-0.1284484655137340,
 263 -0.0004742338006503,0.0127617346949511,-0.0428560878860394,0.0060030732454125,
 264 0.0089182609926781,0.0085353834972860,0.0048464809638033,0.0709740071429510,
 265 0.0029940462557054,-0.0483434904493132,-0.0071713680727884,-0.0036840391887209,
 266 0.0031454003250096,0.0246243550241551,-0.0449551277644180,0.0111449232769393},
 267 {0.0140356721886765,-0.0196518236826680,0.0030517022326582,0.0582672093364850,
 268 -0.0000973895685457,0.0021704767224292,0.0341806268602705,-0.0152035987563018,
 269 -0.0903198657739177,0.0259623214586925,0.0155832497882743,-0.0040543568451651,
 270 0.0036477631918247,-0.0532892744763217,-0.0142569373662724,0.0104500681408622,
 271 0.0103483945857315,0.0679534422398752,-0.0768068882938636,0.0280289727046158}}
 272 ;
 273
 274 static double pameigs[] = {0.0, -0.002350753691875762, -0.002701991863800379,
 275          -0.002931612442853115, -0.004262492032364507, -0.005395980482561625,
 276          -0.007141172690079523, -0.007392844756151318, -0.007781761342200766,
 277          -0.00810032066366362, -0.00875299712761124, -0.01048227332164386,
 278          -0.01109594097332267, -0.01298616073142234, -0.01342036228188581,
 279          -0.01552599145527578, -0.01658762802054814, -0.0174893445623765,
 280          -0.01933280832903272, -0.02206353522613025};
 281
 282 static double pamprobs[20][20] =
 283  {{0.087683339901135, 0.04051291829598762, 0.04087846315185977,
 284    0.04771603459744777, 0.03247095396561266, 0.03784612688594957,
 285    0.0504933695604875, 0.0898249006830755, 0.03285885059543713,
 286    0.0357514442352119, 0.0852464099207521, 0.07910313444070642,
 287    0.01488243946396588, 0.04100101908956829, 0.05158026947089499,
 288    0.06975497205982451, 0.05832757042475474, 0.00931264523877807,
 289    0.03171540880870517, 0.06303972920984541},
 290   {0.01943453646811026, -0.004492574160484092, 0.007694891061220776,
 291    0.01278399096887701, 0.0106157418450234, 0.007542140341575122,
 292    0.01326994069032819, 0.02615565199894889, 0.003123125764490066,
 293    0.002204507682495444, -0.004782898215768979, 0.01204241965177619,
 294    0.0007847400096924341, -0.03043626073172116, 0.01221202591902536,
 295    0.01100527004684405, 0.01116495631339549, -0.0925364931988571,
 296    -0.02622065387931562, 0.00843494142432107},
 297   {0.01855357100209072, 0.01493642835763868, 0.0127983090766285,
 298    0.0200533250704364, -0.1681898360107787, 0.01551657969909255,
 299    0.02128060163107209, 0.03100633591848964, 0.00845480845269879,
 300    0.000927149370785571, 0.00937207565817036, 0.03490557769673472,
 301    0.00300443019551563, -0.02590837220264415, 0.01329376859943192,
 302    0.006854110889741407, 0.01102593860528263, 0.003360844186685888,
 303    -0.03459712356647764, 0.003351477369404443},
 304   {0.02690642688200102, 0.02131745801890152, 0.0143626616005213,
 305    0.02405101425725929, 0.05041008641436849, 0.01430925051050233,
 306    0.02362114036816964, 0.04688381789373886, 0.005250115453626377,
 307    -0.02040112168595516, -0.0942720776915669, 0.03773004996758644,
 308    -0.00822831940782616, -0.1164872809439224, 0.02286281877257392,
 309    0.02849551240669926, 0.01468856796295663, 0.02377110964207936,
 310    -0.094380545436577, -0.02089068498518036},
 311   {0.00930172577225213, 0.01493463068441099, 0.020186920775608,
 312    0.02892154953912524, -0.01224593358361567, 0.01404228329986624,
 313    0.02671186617119041, 0.04537535161795231, 0.02229995804098249,
 314    -0.04635704133961575, -0.1966910360247138, 0.02796648065439046,
 315    -0.02263484732621436, 0.0440490503242072, 0.01148782948302166,
 316    0.01989170531824069, 0.001306805142981245, -0.005676690969116321,
 317    0.07680476281625202, -0.07967537039721849},
 318   {0.06602274245435476, -0.0966661981471856, -0.005241648783844579,
 319    0.00859135188171146, -0.007762129660943368, -0.02888965572526196,
 320    0.003592291525888222, 0.1668410669287673, -0.04082039290551406,
 321    0.005233775047553415, -0.01758244726137135, -0.1493955762326898,
 322    -0.00855819137835548, 0.004211419253492328, 0.01929306335052688,
 323    0.03008056746359405, 0.0190444422412472, 0.005577189741419315,
 324    0.0000874156155112068, 0.02634091459108298},
 325   {0.01933897472880726, 0.05874583569377844, -0.02293534606228405,
 326    -0.07206314017962175, -0.004580681581546643, -0.0628814337610561,
 327    -0.0850783812795136, 0.07988417636610614, -0.0852798990133397,
 328    0.01649047166155952, -0.05416647263757423, 0.1089834536254064,
 329    0.005093403979413865, 0.02520300254161142, 0.0005951431406455604,
 330    0.02441251821224675, 0.02796099482240553, -0.002574933994926502,
 331    -0.007172237553012804, 0.03002455129086954},
 332   {0.04041118479094272, -0.002476225672095412, -0.01494505811263243,
 333    -0.03759443758599911, -0.00892246902492875, -0.003634714029239211,
 334    -0.03085671837973749, -0.126176309029931, 0.005814031139083794,
 335    0.01313561962646063, -0.04760487162503322, -0.0490563712725484,
 336    -0.005082243450421558, -0.01213634309383557, 0.1806666927079249,
 337    0.02111663336185495, 0.02963486860587087, -0.0000175020101657785,
 338    0.01197155383597686, 0.0357526792184636},
 339   {-0.01184769557720525, 0.01582776076338872, -0.006570708266564639,
 340    -0.01471915653734024, 0.00894343616503608, 0.00562664968033149,
 341    -0.01465878888356943, 0.05365282692645818, 0.00893509735776116,
 342    -0.05879312944436473, 0.0806048683392995, -0.007722897986905326,
 343    -0.001819943882718859, 0.0942535573077267, 0.07483883782251654,
 344    0.004354639673913651, -0.02828804845740341, -0.001318222184691827,
 345    -0.07613149604246563, -0.1251675867732172},
 346   {0.00834167031558193, -0.01509357596974962, 0.007098172811092488,
 347    0.03127677418040319, 0.001992448468465455, 0.00915441566808454,
 348    0.03430175973499201, -0.0730648147535803, -0.001402707145575659,
 349    0.04780949194330815, -0.1115035603461273, -0.01292297197609604,
 350    -0.005056270550868528, 0.1112053349612027, -0.03801929822379964,
 351    -0.001191241001736563, 0.01872874622910247, 0.0005314214903865993,
 352    -0.0882576318311789, 0.07607183599610171},
 353   {-0.01539460099727769, 0.04988596184297883, -0.01187240760647617,
 354    -0.06987843637091853, -0.002490472846497859, 0.01009857892494956,
 355    -0.07473588067847209, 0.0906009925879084, 0.1243612446505172,
 356    0.02152806401345371, -0.03504879644860233, -0.06680752427613573,
 357    -0.005574485153629651, 0.001518282948127752, -0.01999168507510701,
 358    -0.01478606199529457, -0.02203749419458996, -0.00132680708294333,
 359    -0.01137505997867614, 0.05332658773667142},
 360   {-0.06104378736432388, 0.0869446603393548, -0.03298331234537257,
 361    0.03128515657456024, 0.003906358569208259, 0.03578694104193928,
 362    0.06241936133189683, 0.06182827284921748, -0.05566564263245907,
 363    0.02640868588189002, -0.01349751243059039, -0.05507866642582638,
 364    -0.006671347738489326, -0.001470096466016046, 0.05185743641479938,
 365    -0.07494697511168257, -0.1175185439057584, -0.001188074094105709,
 366    0.00937934805737347, 0.05024773745437657},
 367   {-0.07252555582124737, -0.116554459356382, 0.003605361887406413,
 368    -0.00836518656029184, 0.004615715410745561, 0.005105376617651312,
 369    -0.00944938657024391, 0.05602449420950007, 0.02722719610561933,
 370    0.01959357494748446, -0.0258655103753962, 0.1440733975689835,
 371    0.01446782819722976, 0.003718896062070054, 0.05825843045655135,
 372    -0.06230154142733073, -0.07833704962300169, 0.003160836143568724,
 373    -0.001169873777936648, 0.03471745590503304},
 374   {-0.03204352258752698, 0.01019272923862322, 0.04509668708733181,
 375    0.05756522429120813, -0.0004601149081726732, -0.0984718150777423,
 376    -0.01107826100664925, -0.005680277810520585, 0.01962359392320817,
 377    0.01550006899131986, 0.05143956925922197, 0.02462476682588468,
 378    -0.0888843861002653, -0.00171553583659411, 0.01606331750661664,
 379    0.001176847743518958, -0.02070972978912828, -0.000341523293579971,
 380    -0.002654732745607882, 0.02075709428885848},
 381   {0.03595199666430258, -0.02800219615234468, -0.04341570015493925,
 382    -0.0748275906176658, 0.0001051403676377422, 0.1137431321746627,
 383    0.005852087565974318, 0.003443037513847801, -0.02481931657706633,
 384    -0.003651181839831423, 0.03195794176786321, 0.04135411406392523,
 385    -0.07562030263210619, 0.001769332364699, -0.01984381173403915,
 386    -0.005029750745010152, 0.02649253902476472, 0.000518085571702734,
 387    0.001062936684474851, 0.01295950668914449},
 388   {-0.16164552322896, -0.0006050035060464324, 0.0258380054414968,
 389    0.003188424740960557, -0.0002058911341821877, 0.03157555987384681,
 390    -0.01678913462596107, 0.03096216145389774, -0.0133791110666919,
 391    0.1125249625204277, -0.00769017706442472, -0.02653938062180483,
 392    -0.002555329863523985, -0.00861833362947954, 0.01775148884754278,
 393    0.02529310679774722, 0.0826243417011238, -0.0001036728183032624,
 394    0.001963562313294209, -0.0935900561309786},
 395   {0.1652394174588469, -0.002814245280784351, -0.0328982001821263,
 396    -0.02000104712964131, 0.0002208121995725443, -0.02733462178511839,
 397    0.02648078162927627, -0.01788316626401427, 0.01630747623755998,
 398    0.1053849023838147, -0.005447706553811218, 0.01810876922536839,
 399    -0.001808914710282444, -0.007687912115607397, -0.01332593672114388,
 400    -0.02110750894891371, -0.07456116592983384, 0.000219072589592394,
 401    0.001270886972191055, -0.1083616930749109},
 402   {0.02453279389716254, -0.005820072356487439, 0.100260287284095,
 403    0.01277522280305745, -0.003184943445296999, 0.05814689527984152,
 404    -0.0934012278200201, -0.03017986487349484, -0.03136625380994165,
 405    0.00988668352785117, -0.00358900410973142, -0.02017443675004764,
 406    0.000915384582922184, -0.001460963415183106, -0.01370112443251124,
 407    0.1130040979284457, -0.1196161771323699, -0.0005800211204222045,
 408    -0.0006153403201024954, 0.00416806428223025},
 409   {-0.0778089244252535, -0.007055161182430869, -0.0349307504860869,
 410    -0.0811915584276571, -0.004689825871599125, -0.03726108871471753,
 411    0.1072225647141469, -0.00917015113070944, 0.01381628985996913,
 412    -0.00123227881492089, 0.001815954515275675, 0.005708744099349901,
 413    -0.0001448985044877925, -0.001306578795561384, -0.006992743514185243,
 414    0.1744720240732789, -0.05353628497814023, -0.0007613684227234787,
 415    -0.0003550282315997644, 0.01340106423804634},
 416   {-0.0159527329868513, -0.007622151568160798, -0.1389875105184963,
 417    0.1165051999914764, -0.002217810389087748, 0.01550003226513692,
 418    -0.07427664222230566, -0.003371438498619264, 0.01385754771325365,
 419    0.004759020167383304, 0.001624078805220564, 0.02011638303109029,
 420    -0.001717827082842178, -0.0007424036708598594, -0.003978884451898934,
 421    0.0866418927301209, -0.01280817739158123, -0.00023039242454603,
 422    0.002309205802479111, 0.0005926106991001195}};
 423
 424
 425 void protdist_uppercase(Char *ch)
 426 {
 427  (*ch) = (isupper(*ch) ? (*ch) : toupper(*ch));
 428 }  /* protdist_uppercase */
 429
 430
 431 void protdist_inputnumbers()
 432 {
 433   /* input the numbers of species and of characters */
 434   long i;
 435
 436   fscanf(infile, "%ld%ld", &spp, &chars);
 437
 438   if (printdata)
 439     fprintf(outfile, "%2ld species, %3ld  positions\n\n", spp, chars);
 440   gnode = (aas **)Malloc(spp * sizeof(aas *));
 441   if (firstset) {
 442     for (i = 0; i < spp; i++)
 443       gnode[i] = (aas *)Malloc(chars * sizeof(aas ));
 444   }
 445   weight = (steparray)Malloc(chars*sizeof(long));
 446   oldweight = (steparray)Malloc(chars*sizeof(long));
 447   category = (steparray)Malloc(chars*sizeof(long));
 448   d      = (double **)Malloc(spp*sizeof(double *));
 449   nayme  = (naym *)Malloc(spp*sizeof(naym));
 450
 451   for (i = 0; i < spp; ++i)
 452     d[i] = (double *)Malloc(spp*sizeof(double));
 453 }  /* protdist_inputnumbers */
 454
 455
 456 void getoptions()
 457 {
 458   /* interactively set options */
 459   long loopcount, loopcount2;
 460   Char ch, ch2;
 461   Char in[100];
 462   boolean done;
 463
 464   if (printdata)
 465     fprintf(outfile, "\nProtein distance algorithm, version %s\n\n",VERSION);
 466   putchar('\n');
 467   weights = false;
 468   printdata = false;
 469   progress = true;
 470   interleaved = true;
 471   similarity = false;
 472   ttratio = 2.0;
 473   whichcode = universal;
 474   whichcat = george;
 475   basesequal = true;
 476   freqa = 0.25;
 477   freqc = 0.25;
 478   freqg = 0.25;
 479   freqt = 0.25;
 480   usejtt = true;
 481   usepmb = false;
 482   usepam = false;
 483   kimura = false;
 484   gama = false;
 485   invar = false;
 486   invarfrac = 0.0;
 487   ease = 0.457;
 488   loopcount = 0;
 489   do {
 490     cleerhome();
 491     printf("\nProtein distance algorithm, version %s\n\n",VERSION);
 492     printf("Settings for this run:\n");
 493     printf("  P  Use JTT, PMB, PAM, Kimura, categories model?  %s\n",
 494            usejtt ? "Jones-Taylor-Thornton matrix" :
 495            usepmb ? "Henikoff/Tillier PMB matrix" :
 496            usepam ? "Dayhoff PAM matrix" :
 497            kimura ? "Kimura formula" :
 498            similarity ? "Similarity table" : "Categories model");
 499     if (!kimura && !similarity) {
 500       printf("  G  Gamma distribution of rates among positions?");
 501       if (gama)
 502         printf("  Yes\n");
 503       else {
 504         if (invar)
 505           printf("  Gamma+Invariant\n");
 506         else
 507           printf("  No\n");
 508       }
 509     }
 510     printf("  C           One category of substitution rates?");
 511     if (!ctgry || categs == 1)
 512       printf("  Yes\n");
 513     else
 514       printf("  %ld categories\n", categs);
 515     printf("  W                    Use weights for positions?");
 516     if (weights)
 517       printf("  Yes\n");
 518     else
 519       printf("  No\n");
 520     if (!(usejtt || usepmb || usepam || kimura || similarity)) {
 521       printf("  U                       Use which genetic code?  %s\n",
 522              (whichcode == universal) ? "Universal"                  :
 523              (whichcode == ciliate)   ? "Ciliate"                    :
 524              (whichcode == mito)      ? "Universal mitochondrial"    :
 525              (whichcode == vertmito)  ? "Vertebrate mitochondrial"   :
 526              (whichcode == flymito)   ? "Fly mitochondrial"        :
 527              (whichcode == yeastmito) ? "Yeast mitochondrial"        : "");
 528       printf("  A          Which categorization of amino acids?  %s\n",
 529              (whichcat == chemical) ? "Chemical"              :
 530              (whichcat == george)   ? "George/Hunt/Barker"    : "Hall");
 531
 532       printf("  E              Prob change category (1.0=easy):%8.4f\n",ease);
 533       printf("  T                Transition/transversion ratio:%7.3f\n",ttratio);
 534       printf("  F                             Base Frequencies:");
 535       if (basesequal)
 536         printf("  Equal\n");
 537       else
 538         printf("%7.3f%6.3f%6.3f%6.3f\n", freqa, freqc, freqg, freqt);
 539     }
 540     printf("  M                   Analyze multiple data sets?");
 541     if (mulsets)
 542       printf("  Yes, %2ld %s\n", datasets,
 543                (justwts ? "sets of weights" : "data sets"));
 544     else
 545       printf("  No\n");
 546     printf("  I                  Input sequences interleaved?  %s\n",
 547            (interleaved ? "Yes" : "No, sequential"));
 548     printf("  0                 Terminal type (IBM PC, ANSI)?  %s\n",
 549            ibmpc ? "IBM PC" :
 550            ansi  ? "ANSI"   : "(none)");
 551     printf("  1            Print out the data at start of run  %s\n",
 552            (printdata ? "Yes" : "No"));
 553     printf("  2          Print indications of progress of run  %s\n",
 554            progress ? "Yes" : "No");
 555     printf("\nAre these settings correct? (type Y or the letter for one to change)\n");
 556     in[0] = '\0';
 557     getstryng(in);
 558     ch=in[0];
 559     if (ch == '\n')
 560       ch = ' ';
 561     protdist_uppercase(&ch);
 562     done = (ch == 'Y');
 563     if (!done) {
 564       if (((strchr("CPGMWI120",ch) != NULL) && (usejtt || usepmb || usepam)) ||
 565           ((strchr("CPMWI120",ch) != NULL) && (kimura || similarity)) ||
 566           ((strchr("CUAPGETFMWI120",ch) != NULL) &&
 567             (! (usejtt || usepmb || usepam || kimura || similarity)))) {
 568         switch (ch) {
 569
 570         case 'U':
 571           printf("Which genetic code?\n");
 572           printf(" type         for\n\n");
 573           printf("   U           Universal\n");
 574           printf("   M           Mitochondrial\n");
 575           printf("   V           Vertebrate mitochondrial\n");
 576           printf("   F           Fly mitochondrial\n");
 577           printf("   Y           Yeast mitochondrial\n\n");
 578           loopcount2 = 0;
 579           do {
 580             printf("type U, M, V, F, or Y\n");
 581             scanf("%c%*[^\n]", &ch);
 582             getchar();
 583             if (ch == '\n')
 584               ch = ' ';
 585             protdist_uppercase(&ch);
 586             countup(&loopcount2, 10);
 587           } while (ch != 'U' && ch != 'M' && ch != 'V' && ch != 'F' && ch != 'Y');
 588           switch (ch) {
 589
 590           case 'U':
 591             whichcode = universal;
 592             break;
 593
 594           case 'M':
 595             whichcode = mito;
 596             break;
 597
 598           case 'V':
 599             whichcode = vertmito;
 600             break;
 601
 602           case 'F':
 603             whichcode = flymito;
 604             break;
 605
 606           case 'Y':
 607             whichcode = yeastmito;
 608             break;
 609           }
 610           break;
 611
 612         case 'A':
 613           printf(
 614             "Which of these categorizations of amino acids do you want to use:\n\n");
 615           printf(
 616             " all have groups: (Glu Gln Asp Asn), (Lys Arg His), (Phe Tyr Trp)\n");
 617           printf(" plus:\n");
 618           printf("George/Hunt/Barker:");
 619           printf(" (Cys), (Met   Val  Leu  Ileu), (Gly  Ala  Ser  Thr    Pro)\n");
 620           printf("Chemical:          ");
 621           printf(" (Cys   Met), (Val  Leu  Ileu    Gly  Ala  Ser  Thr), (Pro)\n");
 622           printf("Hall:              ");
 623           printf(" (Cys), (Met   Val  Leu  Ileu), (Gly  Ala  Ser  Thr), (Pro)\n\n");
 624           printf("Which do you want to use (type C, H, or G)\n");
 625           loopcount2 = 0;
 626           do {
 627             scanf("%c%*[^\n]", &ch);
 628             getchar();
 629             if (ch == '\n')
 630               ch = ' ';
 631             protdist_uppercase(&ch);
 632             countup(&loopcount2, 10);
 633           } while (ch != 'C' && ch != 'H' && ch != 'G');
 634           switch (ch) {
 635
 636           case 'C':
 637             whichcat = chemical;
 638             break;
 639
 640           case 'H':
 641             whichcat = hall;
 642             break;
 643
 644           case 'G':
 645             whichcat = george;
 646             break;
 647           }
 648           break;
 649
 650         case 'C':
 651           ctgry = !ctgry;
 652           if (ctgry) {
 653             initcatn(&categs);
 654             initcategs(categs, rate);
 655           }
 656           break;
 657
 658         case 'W':
 659           weights = !weights;
 660           break;
 661
 662         case 'P':
 663           if (usejtt) {
 664             usejtt = false;
 665             usepmb = true;
 666           } else {
 667             if (usepmb) {
 668               usepmb = false;
 669               usepam = true;
 670             } else {
 671               if (usepam) {
 672                 usepam = false;
 673                 kimura = true;
 674               } else {
 675                 if (kimura) {
 676                   kimura = false;
 677                   similarity = true;
 678                 } else {
 679                   if (similarity)
 680                     similarity = false;
 681                   else
 682                     usejtt = true;
 683                 }
 684               }
 685             }
 686           }
 687           break;
 688
 689         case 'G':
 690           if (!(gama || invar))
 691             gama = true;
 692           else {
 693             if (gama) {
 694               gama = false;
 695               invar = true;
 696             } else {
 697               if (invar)
 698                 invar = false;
 699             }
 700           }
 701           break;
 702
 703
 704         case 'E':
 705           printf("Ease of changing category of amino acid?\n");
 706           loopcount2 = 0;
 707           do {
 708             printf(" (1.0 if no difficulty of changing,\n");
 709             printf(" less if less easy. Can't be negative\n");
 710             scanf("%lf%*[^\n]", &ease);
 711             getchar();
 712             countup(&loopcount2, 10);
 713           } while (ease > 1.0 || ease < 0.0);
 714           break;
 715
 716         case 'T':
 717           loopcount2 = 0;
 718           do {
 719             printf("Transition/transversion ratio?\n");
 720             scanf("%lf%*[^\n]", &ttratio);
 721             getchar();
 722             countup(&loopcount2, 10);
 723           } while (ttratio < 0.0);
 724           break;
 725
 726         case 'F':
 727           loopcount2 = 0;
 728           do {
 729             basesequal = false;
 730             printf("Frequencies of bases A,C,G,T ?\n");
 731             scanf("%lf%lf%lf%lf%*[^\n]", &freqa, &freqc, &freqg, &freqt);
 732             getchar();
 733             if (fabs(freqa + freqc + freqg + freqt - 1.0) >= 1.0e-3)
 734               printf("FREQUENCIES MUST SUM TO 1\n");
 735             countup(&loopcount2, 10);
 736           } while (fabs(freqa + freqc + freqg + freqt - 1.0) >= 1.0e-3);
 737           break;
 738
 739         case 'M':
 740           mulsets = !mulsets;
 741           if (mulsets) {
 742             printf("Multiple data sets or multiple weights?");
 743             loopcount2 = 0;
 744             do {
 745               printf(" (type D or W)\n");
 746               scanf("%c%*[^\n]", &ch2);
 747               getchar();
 748               if (ch2 == '\n')
 749                   ch2 = ' ';
 750               uppercase(&ch2);
 751               countup(&loopcount2, 10);
 752             } while ((ch2 != 'W') && (ch2 != 'D'));
 753             justwts = (ch2 == 'W');
 754             if (justwts)
 755               justweights(&datasets);
 756             else
 757               initdatasets(&datasets);
 758           }
 759           break;
 760
 761         case 'I':
 762           interleaved = !interleaved;
 763           break;
 764
 765         case '0':
 766           if (ibmpc) {
 767             ibmpc = false;
 768             ansi = true;
 769             } else if (ansi)
 770               ansi = false;
 771             else
 772               ibmpc = true;
 773           break;
 774
 775         case '1':
 776           printdata = !printdata;
 777           break;
 778
 779         case '2':
 780           progress = !progress;
 781           break;
 782         }
 783       } else {
 784         if (strchr("CUAPGETFMWI120",ch) == NULL)
 785           printf("Not a possible option!\n");
 786         else
 787           printf("That option not allowed with these settings\n");
 788         printf("\nPress Enter or Return key to continue\n");
 789         getchar();
 790       }
 791     }
 792     countup(&loopcount, 100);
 793   } while (!done);
 794   if (gama || invar) {
 795     loopcount = 0;
 796     do {
 797       printf(
 798 "\nCoefficient of variation of substitution rate among positions (must be positive)\n");
 799       printf(
 800   " In gamma distribution parameters, this is 1/(square root of alpha)\n");
 801       scanf("%lf%*[^\n]", &cvi);
 802       getchar();
 803       countup(&loopcount, 10);
 804     } while (cvi <= 0.0);
 805     cvi = 1.0 / (cvi * cvi);
 806   }
 807   if (invar) {
 808     loopcount = 0;
 809     do {
 810       printf("Fraction of invariant positions?\n");
 811       scanf("%lf%*[^\n]", &invarfrac);
 812       getchar();
 813       countup (&loopcount, 10);
 814     } while ((invarfrac <= 0.0) || (invarfrac >= 1.0));
 815   }
 816 }  /* getoptions */
 817
 818
 819 void transition()
 820 {
 821   /* calculations related to transition-transversion ratio */
 822   double aa, bb, freqr, freqy, freqgr, freqty;
 823
 824   freqr = freqa + freqg;
 825   freqy = freqc + freqt;
 826   freqgr = freqg / freqr;
 827   freqty = freqt / freqy;
 828   aa = ttratio * freqr * freqy - freqa * freqg - freqc * freqt;
 829   bb = freqa * freqgr + freqc * freqty;
 830   xi = aa / (aa + bb);
 831   xv = 1.0 - xi;
 832   if (xi <= 0.0 && xi >= -epsilon)
 833     xi = 0.0;
 834   if (xi < 0.0){
 835     printf("THIS TRANSITION-TRANSVERSION RATIO IS IMPOSSIBLE WITH");
 836     printf(" THESE BASE FREQUENCIES\n");
 837     exxit(-1);}
 838 }  /* transition */
 839
 840
 841 void doinit()
 842 {
 843   /* initializes variables */
 844   protdist_inputnumbers();
 845   getoptions();
 846   transition();
 847 }  /* doinit*/
 848
 849
 850 void printcategories()
 851 { /* print out list of categories of positions */
 852   long i, j;
 853
 854   fprintf(outfile, "Rate categories\n\n");
 855   for (i = 1; i <= nmlngth + 3; i++)
 856     putc(' ', outfile);
 857   for (i = 1; i <= chars; i++) {
 858     fprintf(outfile, "%ld", category[i - 1]);
 859     if (i % 60 == 0) {
 860       putc('\n', outfile);
 861       for (j = 1; j <= nmlngth + 3; j++)
 862         putc(' ', outfile);
 863     } else if (i % 10 == 0)
 864       putc(' ', outfile);
 865   }
 866   fprintf(outfile, "\n\n");
 867 }  /* printcategories */
 868
 869 void reallocchars(void)
 870 {
 871   int i;
 872
 873   free(weight);
 874   free(oldweight);
 875   free(category);
 876   for (i = 0; i < spp; i++) {
 877     free(gnode[i]);
 878     gnode[i] = (aas *)Malloc(chars * sizeof(aas ));
 879   }
 880   weight = (steparray)Malloc(chars*sizeof(long));
 881   oldweight = (steparray)Malloc(chars*sizeof(long));
 882   category = (steparray)Malloc(chars*sizeof(long));
 883 }
 884
 885 void inputoptions()
 886 { /* input the information on the options */
 887   long i;
 888
 889   if (!firstset && !justwts) {
 890     samenumsp(&chars, ith);
 891     reallocchars();
 892   } if (firstset || !justwts) {
 893     for (i = 0; i < chars; i++) {
 894       category[i] = 1;
 895       oldweight[i] = 1;
 896       weight[i] = 1;
 897     }
 898   }
 899   /*  if (!justwts && weights) {*/
 900   if (justwts || weights)
 901     inputweights(chars, oldweight, &weights);
 902   if (printdata)
 903     putc('\n', outfile);
 904   if (usejtt && printdata)
 905     fprintf(outfile, "  Jones-Taylor-Thornton model distance\n");
 906   if (usepmb && printdata)
 907     fprintf(outfile, "  Henikoff/Tillier PMB model distance\n");
 908   if (usepam && printdata)
 909     fprintf(outfile, "  Dayhoff PAM model distance\n");
 910   if (kimura && printdata)
 911     fprintf(outfile, "  Kimura protein distance\n");
 912   if (!(usejtt || usepmb || usepam || kimura || similarity) && printdata)
 913     fprintf(outfile, "  Categories model distance\n");
 914   if (similarity)
 915     fprintf(outfile, "  \n  Table of similarity between sequences\n");
 916   if ((ctgry && categs > 1) && (firstset || !justwts)) {
 917     inputcategs(0, chars, category, categs, "ProtDist");
 918     if (printdata)
 919       printcategs(outfile, chars, category, "Position categories");
 920   } else if (printdata && (categs > 1)) {
 921     fprintf(outfile, "\nPosition category   Rate of change\n\n");
 922     for (i = 1; i <= categs; i++)
 923       fprintf(outfile, "%15ld%13.3f\n", i, rate[i - 1]);
 924     putc('\n', outfile);
 925     printcategories();
 926   }
 927   if (weights && printdata)
 928     printweights(outfile, 0, chars, oldweight, "Positions");
 929 }  /* inputoptions */
 930
 931
 932 void protdist_inputdata()
 933 {
 934   /* input the names and sequences for each species */
 935   long i, j, k, l, aasread=0, aasnew=0;
 936   Char charstate;
 937   boolean allread, done;
 938   aas aa=0;   /* temporary amino acid for input */
 939
 940   if (progress)
 941     putchar('\n');
 942   j = nmlngth + (chars + (chars - 1) / 10) / 2 - 5;
 943   if (j < nmlngth - 1)
 944     j = nmlngth - 1;
 945   if (j > 37)
 946     j = 37;
 947   if (printdata) {
 948     fprintf(outfile, "\nName");
 949     for (i = 1; i <= j; i++)
 950       putc(' ', outfile);
 951     fprintf(outfile, "Sequences\n");
 952     fprintf(outfile, "----");
 953     for (i = 1; i <= j; i++)
 954       putc(' ', outfile);
 955     fprintf(outfile, "---------\n\n");
 956   }
 957   aasread = 0;
 958   allread = false;
 959   while (!(allread)) {
 960     /* eat white space -- if the separator line has spaces on it*/
 961     do {
 962       charstate = gettc(infile);
 963     } while (charstate == ' ' || charstate == '\t');
 964     ungetc(charstate, infile);
 965     if (eoln(infile))
 966       scan_eoln(infile);
 967     i = 1;
 968     while (i <= spp) {
 969       if ((interleaved && aasread == 0) || !interleaved)
 970         initname(i-1);
 971       if (interleaved)
 972         j = aasread;
 973       else
 974         j = 0;
 975       done = false;
 976       while (((!done) && (!(eoln(infile) || eoff(infile))))) {
 977         if (interleaved)
 978           done = true;
 979         while (((j < chars) & (!(eoln(infile) | eoff(infile))))) {
 980           charstate = gettc(infile);
 981           if (charstate == '\n' || charstate == '\t')
 982             charstate = ' ';
 983           if (charstate == ' ' || (charstate >= '0' && charstate <= '9'))
 984             continue;
 985           protdist_uppercase(&charstate);
 986           if ((!isalpha(charstate) && charstate != '.' && charstate != '?' &&
 987                charstate != '-' && charstate != '*') || charstate == 'J' ||
 988               charstate == 'O' || charstate == 'U' || charstate == '.') {
 989         printf("ERROR -- bad amino acid: %c at position %ld of species %3ld\n",
 990                    charstate, j, i);
 991             if (charstate == '.') {
 992           printf("         Periods (.) may not be used as gap characters.\n");
 993           printf("         The correct gap character is (-)\n");
 994             }
 995             exxit(-1);
 996           }
 997           j++;
 998
 999           switch (charstate) {
1000
1001           case 'A':
1002             aa = ala;
1003             break;
1004
1005           case 'B':
1006             aa = asx;
1007             break;
1008
1009           case 'C':
1010             aa = cys;
1011             break;
1012
1013           case 'D':
1014             aa = asp;
1015             break;
1016
1017           case 'E':
1018             aa = glu;
1019             break;
1020
1021           case 'F':
1022             aa = phe;
1023             break;
1024
1025           case 'G':
1026             aa = gly;
1027             break;
1028
1029           case 'H':
1030             aa = his;
1031             break;
1032
1033           case 'I':
1034             aa = ileu;
1035             break;
1036
1037           case 'K':
1038             aa = lys;
1039             break;
1040
1041           case 'L':
1042             aa = leu;
1043             break;
1044
1045           case 'M':
1046             aa = met;
1047             break;
1048
1049           case 'N':
1050             aa = asn;
1051             break;
1052
1053           case 'P':
1054             aa = pro;
1055             break;
1056
1057           case 'Q':
1058             aa = gln;
1059             break;
1060
1061           case 'R':
1062             aa = arg;
1063             break;
1064
1065           case 'S':
1066             aa = ser;
1067             break;
1068
1069           case 'T':
1070             aa = thr;
1071             break;
1072
1073           case 'V':
1074             aa = val;
1075             break;
1076
1077           case 'W':
1078             aa = trp;
1079             break;
1080
1081           case 'X':
1082             aa = unk;
1083             break;
1084
1085           case 'Y':
1086             aa = tyr;
1087             break;
1088
1089           case 'Z':
1090             aa = glx;
1091             break;
1092
1093           case '*':
1094             aa = stop;
1095             break;
1096
1097           case '?':
1098             aa = quest;
1099             break;
1100
1101           case '-':
1102             aa = del;
1103             break;
1104           }
1105           gnode[i - 1][j - 1] = aa;
1106         }
1107         if (interleaved)
1108           continue;
1109         if (j < chars)
1110           scan_eoln(infile);
1111         else if (j == chars)
1112           done = true;
1113       }
1114       if (interleaved && i == 1)
1115         aasnew = j;
1116       scan_eoln(infile);
1117       if ((interleaved && j != aasnew) || ((!interleaved) && j != chars)){
1118         printf("ERROR: SEQUENCES OUT OF ALIGNMENT\n");
1119         exxit(-1);}
1120       i++;
1121     }
1122     if (interleaved) {
1123       aasread = aasnew;
1124       allread = (aasread == chars);
1125     } else
1126       allread = (i > spp);
1127   }
1128   if ( printdata) {
1129     for (i = 1; i <= ((chars - 1) / 60 + 1); i++) {
1130       for (j = 1; j <= spp; j++) {
1131         for (k = 0; k < nmlngth; k++)
1132           putc(nayme[j - 1][k], outfile);
1133         fprintf(outfile, "   ");
1134         l = i * 60;
1135         if (l > chars)
1136           l = chars;
1137         for (k = (i - 1) * 60 + 1; k <= l; k++) {
1138           if (j > 1 && gnode[j - 1][k - 1] == gnode[0][k - 1])
1139             charstate = '.';
1140           else {
1141             switch (gnode[j - 1][k - 1]) {
1142
1143             case ala:
1144               charstate = 'A';
1145               break;
1146
1147             case asx:
1148               charstate = 'B';
1149               break;
1150
1151             case cys:
1152               charstate = 'C';
1153               break;
1154
1155             case asp:
1156               charstate = 'D';
1157               break;
1158
1159             case glu:
1160               charstate = 'E';
1161               break;
1162
1163             case phe:
1164               charstate = 'F';
1165               break;
1166
1167             case gly:
1168               charstate = 'G';
1169               break;
1170
1171             case his:
1172               charstate = 'H';
1173               break;
1174
1175             case ileu:
1176               charstate = 'I';
1177               break;
1178
1179             case lys:
1180               charstate = 'K';
1181               break;
1182
1183             case leu:
1184               charstate = 'L';
1185               break;
1186
1187             case met:
1188               charstate = 'M';
1189               break;
1190
1191             case asn:
1192               charstate = 'N';
1193               break;
1194
1195             case pro:
1196               charstate = 'P';
1197               break;
1198
1199             case gln:
1200               charstate = 'Q';
1201               break;
1202
1203             case arg:
1204               charstate = 'R';
1205               break;
1206
1207             case ser:
1208               charstate = 'S';
1209               break;
1210
1211             case thr:
1212               charstate = 'T';
1213               break;
1214
1215             case val:
1216               charstate = 'V';
1217               break;
1218
1219             case trp:
1220               charstate = 'W';
1221               break;
1222
1223             case tyr:
1224               charstate = 'Y';
1225               break;
1226
1227             case glx:
1228               charstate = 'Z';
1229               break;
1230
1231             case del:
1232               charstate = '-';
1233               break;
1234
1235             case stop:
1236               charstate = '*';
1237               break;
1238
1239             case unk:
1240               charstate = 'X';
1241               break;
1242
1243             case quest:
1244               charstate = '?';
1245               break;
1246
1247             default:        /*cases ser1 and ser2 cannot occur*/
1248               break;
1249             }
1250           }
1251           putc(charstate, outfile);
1252           if (k % 10 == 0 && k % 60 != 0)
1253             putc(' ', outfile);
1254         }
1255         putc('\n', outfile);
1256       }
1257       putc('\n', outfile);
1258     }
1259     putc('\n', outfile);
1260   }
1261   if (printdata)
1262     putc('\n', outfile);
1263 }  /* protdist_inputdata */
1264
1265
1266 void doinput()
1267 { /* reads the input data */
1268   long i;
1269   double sumrates, weightsum;
1270
1271   inputoptions();
1272   if(!justwts || firstset)
1273     protdist_inputdata();
1274   if (!ctgry) {
1275     categs = 1;
1276     rate[0] = 1.0;
1277   }
1278   weightsum = 0;
1279   for (i = 0; i < chars; i++)
1280     weightsum += oldweight[i];
1281   sumrates = 0.0;
1282   for (i = 0; i < chars; i++)
1283     sumrates += oldweight[i] * rate[category[i] - 1];
1284   for (i = 0; i < categs; i++)
1285     rate[i] *= weightsum / sumrates;
1286 }  /* doinput */
1287
1288
1289 void code()
1290 {
1291   /* make up table of the code 1 = u, 2 = c, 3 = a, 4 = g */
1292   long n;
1293   aas b;
1294
1295   trans[0][0][0] = phe;
1296   trans[0][0][1] = phe;
1297   trans[0][0][2] = leu;
1298   trans[0][0][3] = leu;
1299   trans[0][1][0] = ser;
1300   trans[0][1][1] = ser;
1301   trans[0][1][2] = ser;
1302   trans[0][1][3] = ser;
1303   trans[0][2][0] = tyr;
1304   trans[0][2][1] = tyr;
1305   trans[0][2][2] = stop;
1306   trans[0][2][3] = stop;
1307   trans[0][3][0] = cys;
1308   trans[0][3][1] = cys;
1309   trans[0][3][2] = stop;
1310   trans[0][3][3] = trp;
1311   trans[1][0][0] = leu;
1312   trans[1][0][1] = leu;
1313   trans[1][0][2] = leu;
1314   trans[1][0][3] = leu;
1315   trans[1][1][0] = pro;
1316   trans[1][1][1] = pro;
1317   trans[1][1][2] = pro;
1318   trans[1][1][3] = pro;
1319   trans[1][2][0] = his;
1320   trans[1][2][1] = his;
1321   trans[1][2][2] = gln;
1322   trans[1][2][3] = gln;
1323   trans[1][3][0] = arg;
1324   trans[1][3][1] = arg;
1325   trans[1][3][2] = arg;
1326   trans[1][3][3] = arg;
1327   trans[2][0][0] = ileu;
1328   trans[2][0][1] = ileu;
1329   trans[2][0][2] = ileu;
1330   trans[2][0][3] = met;
1331   trans[2][1][0] = thr;
1332   trans[2][1][1] = thr;
1333   trans[2][1][2] = thr;
1334   trans[2][1][3] = thr;
1335   trans[2][2][0] = asn;
1336   trans[2][2][1] = asn;
1337   trans[2][2][2] = lys;
1338   trans[2][2][3] = lys;
1339   trans[2][3][0] = ser;
1340   trans[2][3][1] = ser;
1341   trans[2][3][2] = arg;
1342   trans[2][3][3] = arg;
1343   trans[3][0][0] = val;
1344   trans[3][0][1] = val;
1345   trans[3][0][2] = val;
1346   trans[3][0][3] = val;
1347   trans[3][1][0] = ala;
1348   trans[3][1][1] = ala;
1349   trans[3][1][2] = ala;
1350   trans[3][1][3] = ala;
1351   trans[3][2][0] = asp;
1352   trans[3][2][1] = asp;
1353   trans[3][2][2] = glu;
1354   trans[3][2][3] = glu;
1355   trans[3][3][0] = gly;
1356   trans[3][3][1] = gly;
1357   trans[3][3][2] = gly;
1358   trans[3][3][3] = gly;
1359   if (whichcode == mito)
1360     trans[0][3][2] = trp;
1361   if (whichcode == vertmito) {
1362     trans[0][3][2] = trp;
1363     trans[2][3][2] = stop;
1364     trans[2][3][3] = stop;
1365     trans[2][0][2] = met;
1366   }
1367   if (whichcode == flymito) {
1368     trans[0][3][2] = trp;
1369     trans[2][0][2] = met;
1370     trans[2][3][2] = ser;
1371   }
1372   if (whichcode == yeastmito) {
1373     trans[0][3][2] = trp;
1374     trans[1][0][2] = thr;
1375     trans[2][0][2] = met;
1376   }
1377   n = 0;
1378   for (b = ala; (long)b <= (long)val; b = (aas)((long)b + 1)) {
1379     if (b != ser2) {
1380       n++;
1381       numaa[(long)b - (long)ala] = n;
1382     }
1383   }
1384   numaa[(long)ser - (long)ala] = (long)ser1 - (long)(ala) + 1;
1385 }  /* code */
1386
1387
1388 void protdist_cats()
1389 {
1390   /* define categories of amino acids */
1391   aas b;
1392
1393   /* fundamental subgroups */
1394   cat[0] = 1;                        /* for alanine */
1395   cat[(long)cys - (long)ala] = 1;
1396   cat[(long)met - (long)ala] = 2;
1397   cat[(long)val - (long)ala] = 3;
1398   cat[(long)leu - (long)ala] = 3;
1399   cat[(long)ileu - (long)ala] = 3;
1400   cat[(long)gly - (long)ala] = 4;
1401   cat[0] = 4;
1402   cat[(long)ser - (long)ala] = 4;
1403   cat[(long)thr - (long)ala] = 4;
1404   cat[(long)pro - (long)ala] = 5;
1405   cat[(long)phe - (long)ala] = 6;
1406   cat[(long)tyr - (long)ala] = 6;
1407   cat[(long)trp - (long)ala] = 6;
1408   cat[(long)glu - (long)ala] = 7;
1409   cat[(long)gln - (long)ala] = 7;
1410   cat[(long)asp - (long)ala] = 7;
1411   cat[(long)asn - (long)ala] = 7;
1412   cat[(long)lys - (long)ala] = 8;
1413   cat[(long)arg - (long)ala] = 8;
1414   cat[(long)his - (long)ala] = 8;
1415   if (whichcat == george) {
1416   /* George, Hunt and Barker: sulfhydryl, small hydrophobic, small hydrophilic,
1417                               aromatic, acid/acid-amide/hydrophilic, basic */
1418     for (b = ala; (long)b <= (long)val; b = (aas)((long)b + 1)) {
1419       if (cat[(long)b - (long)ala] == 3)
1420         cat[(long)b - (long)ala] = 2;
1421       if (cat[(long)b - (long)ala] == 5)
1422         cat[(long)b - (long)ala] = 4;
1423     }
1424   }
1425   if (whichcat == chemical) {
1426     /* Conn and Stumpf:  monoamino, aliphatic, heterocyclic,
1427                          aromatic, dicarboxylic, basic */
1428     for (b = ala; (long)b <= (long)val; b = (aas)((long)b + 1)) {
1429       if (cat[(long)b - (long)ala] == 2)
1430         cat[(long)b - (long)ala] = 1;
1431       if (cat[(long)b - (long)ala] == 4)
1432         cat[(long)b - (long)ala] = 3;
1433     }
1434   }
1435   /* Ben Hall's personal opinion */
1436   if (whichcat != hall)
1437     return;
1438   for (b = ala; (long)b <= (long)val; b = (aas)((long)b + 1)) {
1439     if (cat[(long)b - (long)ala] == 3)
1440       cat[(long)b - (long)ala] = 2;
1441   }
1442 }  /* protdist_cats */
1443
1444
1445 void maketrans()
1446 {
1447   /* Make up transition probability matrix from code and category tables */
1448   long i, j, k, m, n, s, nb1, nb2;
1449   double x, sum;
1450   long sub[3], newsub[3];
1451   double f[4], g[4];
1452   aas b1, b2;
1453   double TEMP, TEMP1, TEMP2, TEMP3;
1454
1455   for (i = 0; i <= 19; i++) {
1456     pie[i] = 0.0;
1457     for (j = 0; j <= 19; j++)
1458       prob[i][j] = 0.0;
1459   }
1460   f[0] = freqt;
1461   f[1] = freqc;
1462   f[2] = freqa;
1463   f[3] = freqg;
1464   g[0] = freqc + freqt;
1465   g[1] = freqc + freqt;
1466   g[2] = freqa + freqg;
1467   g[3] = freqa + freqg;
1468   TEMP = f[0];
1469   TEMP1 = f[1];
1470   TEMP2 = f[2];
1471   TEMP3 = f[3];
1472   fracchange = xi * (2 * f[0] * f[1] / g[0] + 2 * f[2] * f[3] / g[2]) +
1473       xv * (1 - TEMP * TEMP - TEMP1 * TEMP1 - TEMP2 * TEMP2 - TEMP3 * TEMP3);
1474   sum = 0.0;
1475   for (i = 0; i <= 3; i++) {
1476     for (j = 0; j <= 3; j++) {
1477       for (k = 0; k <= 3; k++) {
1478         if (trans[i][j][k] != stop)
1479           sum += f[i] * f[j] * f[k];
1480       }
1481     }
1482   }
1483   for (i = 0; i <= 3; i++) {
1484     sub[0] = i + 1;
1485     for (j = 0; j <= 3; j++) {
1486       sub[1] = j + 1;
1487       for (k = 0; k <= 3; k++) {
1488         sub[2] = k + 1;
1489         b1 = trans[i][j][k];
1490         for (m = 0; m <= 2; m++) {
1491           s = sub[m];
1492           for (n = 1; n <= 4; n++) {
1493             memcpy(newsub, sub, sizeof(long) * 3L);
1494             newsub[m] = n;
1495             x = f[i] * f[j] * f[k] / (3.0 * sum);
1496             if (((s == 1 || s == 2) && (n == 3 || n == 4)) ||
1497                 ((n == 1 || n == 2) && (s == 3 || s == 4)))
1498               x *= xv * f[n - 1];
1499             else
1500               x *= xi * f[n - 1] / g[n - 1] + xv * f[n - 1];
1501             b2 = trans[newsub[0] - 1][newsub[1] - 1][newsub[2] - 1];
1502             if (b1 != stop) {
1503               nb1 = numaa[(long)b1 - (long)ala];
1504               pie[nb1 - 1] += x;
1505               if (b2 != stop) {
1506                 nb2 = numaa[(long)b2 - (long)ala];
1507                 if (cat[(long)b1 - (long)ala] != cat[(long)b2 - (long)ala]) {
1508                   prob[nb1 - 1][nb2 - 1] += x * ease;
1509                   prob[nb1 - 1][nb1 - 1] += x * (1.0 - ease);
1510                 } else
1511                   prob[nb1 - 1][nb2 - 1] += x;
1512               } else
1513                 prob[nb1 - 1][nb1 - 1] += x;
1514             }
1515           }
1516         }
1517       }
1518     }
1519   }
1520   for (i = 0; i <= 19; i++)
1521     prob[i][i] -= pie[i];
1522   for (i = 0; i <= 19; i++) {
1523     for (j = 0; j <= 19; j++)
1524       prob[i][j] /= sqrt(pie[i] * pie[j]);
1525   }
1526   /* computes pi^(1/2)*B*pi^(-1/2)  */
1527 }  /* maketrans */
1528
1529
1530 void givens(double (*a)[20], long i, long j, long n, double ctheta,
1531                         double stheta, boolean left)
1532 { /* Givens transform at i,j for 1..n with angle theta */
1533   long k;
1534   double d;
1535
1536   for (k = 0; k < n; k++) {
1537     if (left) {
1538       d = ctheta * a[i - 1][k] + stheta * a[j - 1][k];
1539       a[j - 1][k] = ctheta * a[j - 1][k] - stheta * a[i - 1][k];
1540       a[i - 1][k] = d;
1541     } else {
1542       d = ctheta * a[k][i - 1] + stheta * a[k][j - 1];
1543       a[k][j - 1] = ctheta * a[k][j - 1] - stheta * a[k][i - 1];
1544       a[k][i - 1] = d;
1545     }
1546   }
1547 }  /* givens */
1548
1549
1550 void coeffs(double x, double y, double *c, double *s, double accuracy)
1551 { /* compute cosine and sine of theta */
1552   double root;
1553
1554   root = sqrt(x * x + y * y);
1555   if (root < accuracy) {
1556     *c = 1.0;
1557     *s = 0.0;
1558   } else {
1559     *c = x / root;
1560     *s = y / root;
1561   }
1562 }  /* coeffs */
1563
1564
1565 void tridiag(double (*a)[20], long n, double accuracy)
1566 { /* Givens tridiagonalization */
1567   long i, j;
1568   double s, c;
1569
1570   for (i = 2; i < n; i++) {
1571     for (j = i + 1; j <= n; j++) {
1572       coeffs(a[i - 2][i - 1], a[i - 2][j - 1], &c, &s,accuracy);
1573       givens(a, i, j, n, c, s, true);
1574       givens(a, i, j, n, c, s, false);
1575       givens(eigvecs, i, j, n, c, s, true);
1576     }
1577   }
1578 }  /* tridiag */
1579
1580
1581 void shiftqr(double (*a)[20], long n, double accuracy)
1582 { /* QR eigenvalue-finder */
1583   long i, j;
1584   double approx, s, c, d, TEMP, TEMP1;
1585
1586   for (i = n; i >= 2; i--) {
1587     do {
1588       TEMP = a[i - 2][i - 2] - a[i - 1][i - 1];
1589       TEMP1 = a[i - 1][i - 2];
1590       d = sqrt(TEMP * TEMP + TEMP1 * TEMP1);
1591       approx = a[i - 2][i - 2] + a[i - 1][i - 1];
1592       if (a[i - 1][i - 1] < a[i - 2][i - 2])
1593         approx = (approx - d) / 2.0;
1594       else
1595         approx = (approx + d) / 2.0;
1596       for (j = 0; j < i; j++)
1597         a[j][j] -= approx;
1598       for (j = 1; j < i; j++) {
1599         coeffs(a[j - 1][j - 1], a[j][j - 1], &c, &s, accuracy);
1600         givens(a, j, j + 1, i, c, s, true);
1601         givens(a, j, j + 1, i, c, s, false);
1602         givens(eigvecs, j, j + 1, n, c, s, true);
1603       }
1604       for (j = 0; j < i; j++)
1605         a[j][j] += approx;
1606     } while (fabs(a[i - 1][i - 2]) > accuracy);
1607   }
1608 }  /* shiftqr */
1609
1610
1611 void qreigen(double (*prob)[20], long n)
1612 { /* QR eigenvector/eigenvalue method for symmetric matrix */
1613   double accuracy;
1614   long i, j;
1615
1616   accuracy = 1.0e-6;
1617   for (i = 0; i < n; i++) {
1618     for (j = 0; j < n; j++)
1619       eigvecs[i][j] = 0.0;
1620     eigvecs[i][i] = 1.0;
1621   }
1622   tridiag(prob, n, accuracy);
1623   shiftqr(prob, n, accuracy);
1624   for (i = 0; i < n; i++)
1625     eig[i] = prob[i][i];
1626   for (i = 0; i <= 19; i++) {
1627     for (j = 0; j <= 19; j++)
1628       prob[i][j] = sqrt(pie[j]) * eigvecs[i][j];
1629   }
1630   /* prob[i][j] is the value of U' times pi^(1/2) */
1631 }  /* qreigen */
1632
1633
1634 void jtteigen()
1635 { /* eigenanalysis for JTT matrix, precomputed */
1636   memcpy(prob,jttprobs,sizeof(jttprobs));
1637   memcpy(eig,jtteigs,sizeof(jtteigs));
1638   fracchange = 0.01;
1639 }  /* jtteigen */
1640
1641
1642 void pmbeigen()
1643 { /* eigenanalysis for PMB matrix, precomputed */
1644   memcpy(prob,pmbprobs,sizeof(pmbprobs));
1645   memcpy(eig,pmbeigs,sizeof(pmbeigs));
1646   fracchange = 1.0;
1647 }  /* pmbeigen */
1648
1649
1650 void pameigen()
1651 { /* eigenanalysis for PAM matrix, precomputed */
1652   memcpy(prob,pamprobs,sizeof(pamprobs));
1653   memcpy(eig,pameigs,sizeof(pameigs));
1654   fracchange = 0.01;
1655 }  /* pameigen */
1656
1657
1658 void predict(long nb1, long nb2, long cat)
1659 { /* make contribution to prediction of this aa pair */
1660   long m;
1661   double TEMP;
1662
1663   for (m = 0; m <= 19; m++) {
1664     if (gama || invar)
1665       elambdat = exp(-cvi*log(1.0-rate[cat-1]*tt*(eig[m]/(1.0-invarfrac))/cvi));
1666     else
1667       elambdat = exp(rate[cat-1]*tt * eig[m]);
1668     q = prob[m][nb1 - 1] * prob[m][nb2 - 1] * elambdat;
1669     p += q;
1670     if (!gama && !invar)
1671       dp += rate[cat-1]*eig[m] * q;
1672     else
1673       dp += (rate[cat-1]*eig[m]/(1.0-rate[cat-1]*tt*(eig[m]/(1.0-invarfrac))/cvi)) * q;
1674     TEMP = eig[m];
1675     if (!gama && !invar)
1676       d2p += TEMP * TEMP * q;
1677     else
1678       d2p += (rate[cat-1]*rate[cat-1]*eig[m]*eig[m]*(1.0+1.0/cvi)/
1679               ((1.0-rate[cat-1]*tt*eig[m]/cvi)
1680               *(1.0-rate[cat-1]*tt*eig[m]/cvi))) * q;
1681   }
1682   if (nb1 == nb2) {
1683     p *= (1.0 - invarfrac);
1684     p += invarfrac;
1685   }
1686   dp *= (1.0 - invarfrac);
1687   d2p *= (1.0 - invarfrac);
1688 }  /* predict */
1689
1690 void makedists()
1691 { /* compute the distances */
1692   long i, j, k, m, n, itterations, nb1, nb2, cat;
1693   double delta, lnlike, slope, curv;
1694   boolean neginfinity, inf, overlap;
1695   aas b1, b2;
1696
1697   if (!(printdata || similarity))
1698     fprintf(outfile, "%5ld\n", spp);
1699   if (progress)
1700     printf("Computing distances:\n");
1701   for (i = 1; i <= spp; i++) {
1702     if (progress)
1703       printf("  ");
1704     if (progress) {
1705       for (j = 0; j < nmlngth; j++)
1706         putchar(nayme[i - 1][j]);
1707     }
1708     if (progress) {
1709       printf("   ");
1710       fflush(stdout);
1711     }
1712     if (similarity)
1713       d[i-1][i-1] = 1.0;
1714     else
1715       d[i-1][i-1] = 0.0;
1716     for (j = 0; j <= i - 2; j++) {
1717       if (!(kimura || similarity)) {
1718         if (usejtt || usepmb || usepam)
1719           tt = 0.1/fracchange;
1720         else
1721           tt = 1.0;
1722         delta = tt / 2.0;
1723         itterations = 0;
1724         inf = false;
1725         do {
1726           lnlike = 0.0;
1727           slope = 0.0;
1728           curv = 0.0;
1729           neginfinity = false;
1730           overlap = false;
1731           for (k = 0; k < chars; k++) {
1732             if (oldweight[k] > 0) {
1733               cat = category[k];
1734               b1 = gnode[i - 1][k];
1735               b2 = gnode[j][k];
1736               if (b1 != stop && b1 != del && b1 != quest && b1 != unk &&
1737                   b2 != stop && b2 != del && b2 != quest && b2 != unk) {
1738                 overlap = true;
1739                 p = 0.0;
1740                 dp = 0.0;
1741                 d2p = 0.0;
1742                 nb1 = numaa[(long)b1 - (long)ala];
1743                 nb2 = numaa[(long)b2 - (long)ala];
1744                 if (b1 != asx && b1 != glx && b2 != asx && b2 != glx)
1745                   predict(nb1, nb2, cat);
1746                 else {
1747                   if (b1 == asx) {
1748                     if (b2 == asx) {
1749                       predict(3L, 3L, cat);
1750                       predict(3L, 4L, cat);
1751                       predict(4L, 3L, cat);
1752                       predict(4L, 4L, cat);
1753                     } else {
1754                       if (b2 == glx) {
1755                         predict(3L, 6L, cat);
1756                         predict(3L, 7L, cat);
1757                         predict(4L, 6L, cat);
1758                         predict(4L, 7L, cat);
1759                       } else {
1760                         predict(3L, nb2, cat);
1761                         predict(4L, nb2, cat);
1762                       }
1763                     }
1764                   } else {
1765                     if (b1 == glx) {
1766                       if (b2 == asx) {
1767                         predict(6L, 3L, cat);
1768                         predict(6L, 4L, cat);
1769                         predict(7L, 3L, cat);
1770                         predict(7L, 4L, cat);
1771                       } else {
1772                         if (b2 == glx) {
1773                           predict(6L, 6L, cat);
1774                           predict(6L, 7L, cat);
1775                           predict(7L, 6L, cat);
1776                           predict(7L, 7L, cat);
1777                         } else {
1778                           predict(6L, nb2, cat);
1779                           predict(7L, nb2, cat);
1780                         }
1781                       }
1782                     } else {
1783                       if (b2 == asx) {
1784                         predict(nb1, 3L, cat);
1785                         predict(nb1, 4L, cat);
1786                         predict(nb1, 3L, cat);
1787                         predict(nb1, 4L, cat);
1788                       } else if (b2 == glx) {
1789                         predict(nb1, 6L, cat);
1790                         predict(nb1, 7L, cat);
1791                         predict(nb1, 6L, cat);
1792                         predict(nb1, 7L, cat);
1793                       }
1794                     }
1795                   }
1796                 }
1797                 if (p <= 0.0)
1798                   neginfinity = true;
1799                 else {
1800                   lnlike += oldweight[k]*log(p);
1801                   slope += oldweight[k]*dp / p;
1802                   curv += oldweight[k]*(d2p / p - dp * dp / (p * p));
1803                 }
1804               }
1805             }
1806           }
1807           itterations++;
1808           if (!overlap){
1809             printf("\nWARNING: NO OVERLAP BETWEEN SEQUENCES %ld AND %ld; -1.0 WAS WRITTEN\n", i, j+1);
1810             tt = -1.0/fracchange;
1811             itterations = 20;
1812             inf = true;
1813           } else if (!neginfinity) {
1814             if (curv < 0.0) {
1815               tt -= slope / curv;
1816               if (tt > 10000.0) {
1817                 printf("\nWARNING: INFINITE DISTANCE BETWEEN SPECIES %ld AND %ld; -1.0 WAS WRITTEN\n", i, j+1);
1818                 tt = -1.0/fracchange;
1819                 inf = true;
1820                 itterations = 20;
1821               }
1822             }
1823             else {
1824               if ((slope > 0.0 && delta < 0.0) || (slope < 0.0 && delta > 0.0))
1825                 delta /= -2;
1826               tt += delta;
1827             }
1828           } else {
1829             delta /= -2;
1830             tt += delta;
1831           }
1832           if (tt < protepsilon && !inf)
1833             tt = protepsilon;
1834         } while (itterations != 20);
1835       } else {
1836         m = 0;
1837         n = 0;
1838         for (k = 0; k < chars; k++) {
1839           b1 = gnode[i - 1][k];
1840           b2 = gnode[j][k];
1841           if ((((long)b1 <= (long)val) || ((long)b1 == (long)ser))
1842            && (((long)b2 <= (long)val) || ((long)b2 == (long)ser))) {
1843             if (b1 == b2)
1844               m++;
1845             n++;
1846           }
1847         }
1848         p = 1 - (double)m / n;
1849         if (kimura) {
1850           dp = 1.0 - p - 0.2 * p * p;
1851           if (dp < 0.0) {
1852             printf(
1853 "\nDISTANCE BETWEEN SEQUENCES %3ld AND %3ld IS TOO LARGE FOR KIMURA FORMULA\n",
1854               i, j + 1);
1855             tt = -1.0;
1856           } else
1857             tt = -log(dp);
1858         } else {              /* if similarity */
1859             tt = 1.0 - p;
1860         }
1861       }
1862       d[i - 1][j] = fracchange * tt;
1863       d[j][i - 1] = d[i - 1][j];
1864       if (progress) {
1865         putchar('.');
1866         fflush(stdout);
1867       }
1868     }
1869     if (progress) {
1870       putchar('\n');
1871       fflush(stdout);
1872     }
1873   }
1874   if (!similarity) {
1875     for (i = 0; i < spp; i++) {
1876       for (j = 0; j < nmlngth; j++)
1877         putc(nayme[i][j], outfile);
1878       k = spp;
1879       for (j = 1; j <= k; j++) {
1880         fprintf(outfile, "%10.6f", d[i][j - 1]);
1881         if ((j + 1) % 7 == 0 && j < k)
1882           putc('\n', outfile);
1883       }
1884       putc('\n', outfile);
1885     }
1886   } else {
1887     for (i = 0; i < spp; i += 6) {
1888       if ((i+6) < spp)
1889         n = i+6;
1890       else
1891         n = spp;
1892       fprintf(outfile, "            ");
1893       for (j = i; j < n ; j++) {
1894         for (k = 0; k < (nmlngth-2); k++)
1895           putc(nayme[j][k], outfile);
1896         putc(' ', outfile);
1897         putc(' ', outfile);
1898       }
1899       putc('\n', outfile);
1900       for (j = 0; j < spp; j++) {
1901         for (k = 0; k < nmlngth; k++)
1902           putc(nayme[j][k], outfile);
1903         if ((i+6) < spp)
1904           n = i+6;
1905         else
1906           n = spp;
1907         for (k = i; k < n ; k++)
1908           fprintf(outfile, "%10.6f", d[j][k]);
1909         putc('\n', outfile);
1910       }
1911       putc('\n', outfile);
1912     }
1913   }
1914   if (progress)
1915     printf("\nOutput written to file \"%s\"\n\n", outfilename);
1916 }  /* makedists */
1917
1918
1919 int main(int argc, Char *argv[])
1920 {  /* ML Protein distances by PMB, JTT, PAM or categories model */
1921 #ifdef MAC
1922    argc = 1;   /* macsetup("Protdist",""); */
1923    argv[0] = "Protdist";
1924 #endif
1925   init(argc, argv);
1926   openfile(&infile,INFILE,"input file","r",argv[0],infilename);
1927   openfile(&outfile,OUTFILE,"output file","w",argv[0],outfilename);
1928   ibmpc = IBMCRT;
1929   ansi = ANSICRT;
1930   mulsets = false;
1931   datasets = 1;
1932   firstset = true;
1933   doinit();
1934   if (!(kimura || similarity))
1935     code();
1936   if (!(usejtt || usepmb || usepam ||  kimura || similarity)) {
1937     protdist_cats();
1938     maketrans();
1939     qreigen(prob, 20L);
1940   } else {
1941     if (kimura || similarity)
1942       fracchange = 1.0;
1943     else {
1944       if (usejtt)
1945         jtteigen();
1946       else {
1947         if (usepmb)
1948           pmbeigen();
1949         else
1950           pameigen();
1951         }
1952     }
1953   }
1954   if (ctgry)
1955     openfile(&catfile,CATFILE,"categories file","r",argv[0],catfilename);
1956   if (weights || justwts)
1957     openfile(&weightfile,WEIGHTFILE,"weights file","r",argv[0],weightfilename);
1958   for (ith = 1; ith <= datasets; ith++) {
1959     doinput();
1960     if (ith == 1)
1961       firstset = false;
1962     if ((datasets > 1) && progress)
1963       printf("\nData set # %ld:\n\n", ith);
1964     makedists();
1965   }
1966   FClose(outfile);
1967   FClose(infile);
1968 #ifdef MAC
1969   fixmacfile(outfilename);
1970 #endif
1971   return 0;
1972 }  /* Protein distances */
1973