source: trunk/GDE/PHYLIP/protdist.c

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