source: branches/lib/GDE/RAxML/parsePartitions.c

Last change on this file was 10409, checked in by aboeckma, 12 years ago

Updated raxml to current version
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 34.6 KB
Line 
1/*  RAxML-VI-HPC (version 2.2) a program for sequential and parallel estimation of phylogenetic trees
2 *  Copyright August 2006 by Alexandros Stamatakis
3 *
4 *  Partially derived from
5 *  fastDNAml, a program for estimation of phylogenetic trees from sequences by Gary J. Olsen
6 * 
7 *  and
8 *
9 *  Programs of the PHYLIP package by Joe Felsenstein.
10 *
11 *  This program is free software; you may redistribute it and/or modify its
12 *  under the terms of the GNU General Public License as published by the Free
13 *  Software Foundation; either version 2 of the License, or (at your option)
14 *  any later version.
15 *
16 *  This program is distributed in the hope that it will be useful, but
17 *  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 *  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
19 *  for more details.
20 *
21 *
22 *  For any other enquiries send an Email to Alexandros Stamatakis
23 *  Alexandros.Stamatakis@epfl.ch
24 *
25 *  When publishing work that is based on the results from RAxML-VI-HPC please cite:
26 *
27 *  Alexandros Stamatakis:"RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models".
28 *  Bioinformatics 2006; doi: 10.1093/bioinformatics/btl446
29 */
30
31
32#ifndef WIN32
33#include <sys/times.h>
34#include <sys/types.h>
35#include <sys/time.h>
36#include <unistd.h>
37#endif
38
39#include <math.h>
40#include <time.h>
41#include <stdlib.h>
42#include <stdio.h>
43#include <ctype.h>
44#include <string.h>
45#include <strings.h>
46
47#include "axml.h"
48
49/*****************************FUNCTIONS FOR READING MULTIPLE MODEL SPECIFICATIONS************************************************/
50
51
52extern char modelFileName[1024];
53extern char excludeFileName[1024];
54extern char secondaryStructureFileName[1024];
55
56
57extern char seq_file[1024];
58
59extern char *protModels[NUM_PROT_MODELS];
60
61static boolean lineContainsOnlyWhiteChars(char *line)
62{
63  int i, n = strlen(line);
64
65  if(n == 0)
66    return TRUE;
67
68  for(i = 0; i < n; i++)
69    {
70      if(!whitechar(line[i]))
71        return FALSE;
72    }
73  return TRUE;
74}
75
76
77static int isNum(char c)
78{
79 
80  return (c == '0' || c == '1' || c == '2' || c == '3' || c == '4' ||
81          c == '5' || c == '6' || c == '7' || c == '8' || c == '9');
82}
83
84
85static void skipWhites(char **ch)
86{
87  while(**ch == ' ' || **ch == '\t')
88    *ch = *ch + 1;
89}
90
91static void analyzeIdentifier(char **ch, int modelNumber, tree *tr)
92{
93  char 
94    ident[2048] = "",
95    model[2048] = "",
96    thisModel[2048] = "";
97 
98  int 
99    i = 0, 
100    n, 
101    j,
102    containsComma = 0;
103
104  while(**ch != '=')
105    {
106      if(**ch != ' ' && **ch != '\t')
107        {
108          ident[i] = **ch;     
109          i++;
110        }
111      *ch = *ch + 1;
112    }
113 
114  n = i;
115  i = 0;
116 
117  for(i = 0; i < n; i++)
118    if(ident[i] == ',') 
119      containsComma = 1;
120
121  if(!containsComma)
122    {
123      printf("Error, model file must have format: DNA or AA model, then a comma, and then the partition name\n");
124      exit(-1);
125    }
126  else
127    {
128      boolean
129        useProteinSubstitutionFile = FALSE,
130        found = FALSE;
131     
132      int 
133        openBracket = 0,
134        closeBracket = 0,
135        openPos = 0,
136        closePos = 0;
137           
138      i = 0;
139     
140      while(ident[i] != ',')
141        {       
142          if(ident[i] == '[')
143            {
144              openPos = i;
145              openBracket++;
146            }
147          if(ident[i] == ']')
148            {
149              closePos = i;
150              closeBracket++;
151            }
152          model[i] = ident[i];
153          i++;
154        }     
155
156      if(closeBracket > 0 || openBracket > 0)
157        {
158          if((closeBracket == 1) && (openBracket == 1) && (openPos < closePos))
159            useProteinSubstitutionFile = TRUE;
160          else
161            {
162              printf("\nError: Apparently you want to specify a user-defined protein substitution model that shall be read from file\n");
163              printf("It must be enclosed in opening and closing bracktes like this: [fileName]\n\n");
164              printf("you specified: %s\n\n", model);
165              exit(-1);
166            }
167        }
168     
169      if(useProteinSubstitutionFile)
170        {
171          char 
172            protFileName[2048] = "";
173
174          int 
175            pos,
176            k,
177            lower = 0,
178            upper = i - 1;
179
180          while(model[lower] == '[' || model[lower] == ' ')
181            lower++;
182
183          while(model[upper] == ']' || model[upper] == ' ')
184            upper--;
185         
186          assert(lower < upper);
187
188          for(k = lower, pos = 0; k <= upper; k++, pos++)
189            protFileName[pos] = model[k];
190         
191          protFileName[pos] = '\0';
192
193         
194
195          if(!filexists(protFileName))
196            {
197              printf("\n\ncustom protein substitution file [%s] you want to use does not exist!\n", protFileName);
198              printf("you need to specify the full path\n");
199              printf("the file name shall not contain blanks!\n\n");
200              exit(0);
201            }
202         
203
204          strcpy(tr->initialPartitionData[modelNumber].proteinSubstitutionFileName, protFileName);
205          /*printf("%s \n", tr->initialPartitionData[modelNumber].proteinSubstitutionFileName);*/
206         
207          tr->initialPartitionData[modelNumber].protModels = PROT_FILE;           
208          tr->initialPartitionData[modelNumber].usePredefinedProtFreqs  = TRUE;
209          tr->initialPartitionData[modelNumber].dataType   = AA_DATA;
210        }
211      else
212        {
213          /* AA */
214         
215          for(i = 0; i < NUM_PROT_MODELS && !found; i++)
216            {   
217              strcpy(thisModel, protModels[i]);
218             
219              if(strcasecmp(model, thisModel) == 0)
220                {                     
221                  tr->initialPartitionData[modelNumber].protModels = i;           
222                  tr->initialPartitionData[modelNumber].usePredefinedProtFreqs = TRUE;
223                  tr->initialPartitionData[modelNumber].dataType   = AA_DATA;             
224                  found = TRUE;
225                }
226             
227
228              if(i != GTR && i != GTR_UNLINKED)
229                {
230                  strcpy(thisModel, protModels[i]);
231                  strcat(thisModel, "F");
232                 
233                  if(strcasecmp(model, thisModel) == 0)
234                    {         
235                      tr->initialPartitionData[modelNumber].protModels = i;               
236                      tr->initialPartitionData[modelNumber].usePredefinedProtFreqs  = FALSE;
237                      tr->initialPartitionData[modelNumber].dataType   = AA_DATA;                   
238                      found = TRUE;
239                    }
240                }
241             
242              if(found && (tr->initialPartitionData[modelNumber].protModels == GTR || tr->initialPartitionData[modelNumber].protModels == GTR_UNLINKED))
243                tr->initialPartitionData[modelNumber].usePredefinedProtFreqs  = FALSE;                         
244            }
245         
246          if(!found)
247            {                     
248              if(strcasecmp(model, "DNA") == 0)
249                {                     
250                  tr->initialPartitionData[modelNumber].protModels = -1;                 
251                  tr->initialPartitionData[modelNumber].usePredefinedProtFreqs  = FALSE;
252                  tr->initialPartitionData[modelNumber].dataType   = DNA_DATA;
253                 
254                  found = TRUE;
255                }         
256              else
257                {                 
258                  if(strcasecmp(model, "BIN") == 0)
259                    {                 
260                      tr->initialPartitionData[modelNumber].protModels = -1;             
261                      tr->initialPartitionData[modelNumber].usePredefinedProtFreqs  = FALSE;
262                      tr->initialPartitionData[modelNumber].dataType   = BINARY_DATA;
263                     
264                      found = TRUE;
265                    }
266                  else
267                    {
268                      if(strcasecmp(model, "MULTI") == 0)
269                        {                     
270                          tr->initialPartitionData[modelNumber].protModels = -1;                 
271                          tr->initialPartitionData[modelNumber].usePredefinedProtFreqs  = FALSE;
272                          tr->initialPartitionData[modelNumber].dataType   = GENERIC_32;
273                         
274                          found = TRUE;
275                        }
276                      else
277                        {
278                          if(strcasecmp(model, "CODON") == 0)
279                            {                 
280                              tr->initialPartitionData[modelNumber].protModels = -1;             
281                              tr->initialPartitionData[modelNumber].usePredefinedProtFreqs  = FALSE;
282                              tr->initialPartitionData[modelNumber].dataType   = GENERIC_64;
283                             
284                              found = TRUE;
285                            }
286                        }
287                    }
288                 
289                }
290            }
291
292          if(!found)
293            {
294              printf("ERROR: you specified the unknown model %s for partition %d\n", model, modelNumber);
295              exit(-1);
296            }
297        }
298         
299      i = 0;
300      while(ident[i++] != ',');     
301         
302      tr->initialPartitionData[modelNumber].partitionName = (char*)rax_malloc((n - i + 1) * sizeof(char));         
303         
304      j = 0;
305      while(i < n)     
306        tr->initialPartitionData[modelNumber].partitionName[j++] =  ident[i++];
307     
308      tr->initialPartitionData[modelNumber].partitionName[j] = '\0';                     
309    }
310}
311
312
313
314static void setModel(int model, int position, int *a)
315{
316  if(a[position] == -1)
317    a[position] = model;
318  else
319    {
320      printf("ERROR trying to assign model %d to position %d \n", model, position);
321      printf("while already model %d has been assigned to this position\n", a[position]);
322      exit(-1);
323    }     
324}
325
326
327static int myGetline(char **lineptr, int *n, FILE *stream)
328{
329  char *line, *p;
330  int size, copy, len;
331  int chunkSize = 256 * sizeof(char);
332
333   if (*lineptr == NULL || *n < 2) 
334    {
335      line = (char *)rax_realloc(*lineptr, chunkSize, FALSE);
336      if (line == NULL)
337        return -1;
338      *lineptr = line;
339      *n = chunkSize;
340    }
341
342   line = *lineptr;
343   size = *n;
344 
345   copy = size;
346   p = line;
347   
348   while(1)
349     {
350       while (--copy > 0)
351         {
352           register int c = getc(stream);
353           if (c == EOF)
354             goto lose;
355           else
356             {
357               *p++ = c;
358               if(c == '\n' || c == '\r')       
359                 goto win;
360             }
361         }
362
363       /* Need to enlarge the line buffer.  */
364       len = p - line;
365       size *= 2;
366       line = rax_realloc (line, size, FALSE);
367       if (line == NULL)
368         goto lose;
369       *lineptr = line;
370       *n = size;
371       p = line + len;
372       copy = size - len;
373     }
374   
375 lose:
376  if (p == *lineptr)
377    return -1;
378  /* Return a partial line since we got an error in the middle.  */
379 win:
380  *p = '\0';
381  return p - *lineptr;
382}
383
384
385
386void parsePartitions(analdef *adef, rawdata *rdta, tree *tr)
387{
388  FILE *f; 
389  int numberOfModels = 0; 
390  int nbytes = 0;
391  char *ch;
392  char *cc = (char *)NULL;
393  char **p_names;
394  int n, i, l;
395  int lower, upper, modulo;
396  char buf[256];
397  int **partitions;
398  int pairsCount;
399  int as, j;
400  int k; 
401
402  f = myfopen(modelFileName, "rb");   
403
404 
405  while(myGetline(&cc, &nbytes, f) > -1)
406    {     
407      if(!lineContainsOnlyWhiteChars(cc))
408        {
409          numberOfModels++;
410        }
411      if(cc)
412        rax_free(cc);
413      cc = (char *)NULL;
414    }     
415 
416  rewind(f);
417     
418  p_names = (char **)rax_malloc(sizeof(char *) * numberOfModels);
419  partitions = (int **)rax_malloc(sizeof(int *) * numberOfModels);
420     
421 
422 
423  tr->initialPartitionData = (pInfo*)rax_malloc(sizeof(pInfo) * numberOfModels);
424
425     
426  for(i = 0; i < numberOfModels; i++) 
427    {     
428      tr->initialPartitionData[i].protModels = adef->proteinMatrix;
429      tr->initialPartitionData[i].usePredefinedProtFreqs  = adef->protEmpiricalFreqs;
430      tr->initialPartitionData[i].dataType   = -1;
431    }
432
433  for(i = 0; i < numberOfModels; i++)   
434    partitions[i] = (int *)NULL;
435   
436  i = 0;
437  while(myGetline(&cc, &nbytes, f) > -1)
438    {         
439      if(!lineContainsOnlyWhiteChars(cc))
440        {
441          n = strlen(cc);       
442          p_names[i] = (char *)rax_malloc(sizeof(char) * (n + 1));
443          strcpy(&(p_names[i][0]), cc); 
444          i++;
445        }
446      if(cc)
447        rax_free(cc);
448      cc = (char *)NULL;
449    }         
450
451 
452
453  for(i = 0; i < numberOfModels; i++)
454    {         
455     
456      ch = p_names[i];     
457      pairsCount = 0;
458      skipWhites(&ch);
459     
460      if(*ch == '=')
461        {
462          printf("Identifier missing prior to '=' in %s\n", p_names[i]);
463          exit(-1);
464        }
465     
466      analyzeIdentifier(&ch, i, tr);
467      ch++;
468           
469    numberPairs:
470      pairsCount++;
471      partitions[i] = (int *)rax_realloc((void *)partitions[i], (1 + 3 * pairsCount) * sizeof(int), FALSE);
472      partitions[i][0] = pairsCount;
473      partitions[i][3 + 3 * (pairsCount - 1)] = -1;     
474     
475      skipWhites(&ch);
476     
477      if(!isNum(*ch))
478        {
479          printf("%c Number expected in %s\n", *ch, p_names[i]);
480          exit(-1);
481        }   
482     
483      l = 0;
484      while(isNum(*ch))         
485        {
486          /*printf("%c", *ch);*/
487          buf[l] = *ch;
488          ch++; 
489          l++;
490        }
491      buf[l] = '\0';
492      lower = atoi(buf);
493      partitions[i][1 + 3 * (pairsCount - 1)] = lower;   
494     
495      skipWhites(&ch);
496     
497      /* NEW */
498     
499      if((*ch != '-') && (*ch != ','))
500        {
501          if(*ch == '\0' || *ch == '\n' || *ch == '\r')
502            {
503              upper = lower;
504              goto SINGLE_NUMBER;
505            }
506          else
507            {
508              printf("'-' or ',' expected in %s\n", p_names[i]);
509              exit(-1);
510            }
511        }       
512     
513      if(*ch == ',')
514        {           
515          upper = lower;
516          goto SINGLE_NUMBER;
517        }
518     
519      /* END NEW */
520     
521      ch++;   
522     
523      skipWhites(&ch);
524     
525      if(!isNum(*ch))
526        {
527          printf("%c Number expected in %s\n", *ch, p_names[i]);
528          exit(-1);
529        }   
530     
531      l = 0;
532      while(isNum(*ch))
533        {   
534          buf[l] = *ch;
535          ch++; 
536          l++;
537        }
538      buf[l] = '\0';
539      upper = atoi(buf);     
540    SINGLE_NUMBER:
541      partitions[i][2 + 3 * (pairsCount - 1)] = upper;           
542     
543      if(upper < lower)
544        {
545          printf("Upper bound %d smaller than lower bound %d for this partition: %s\n", upper, lower,  p_names[i]);
546          exit(-1);
547        }
548     
549      skipWhites(&ch);
550     
551      if(*ch == '\0' || *ch == '\n' || *ch == '\r') /* PC-LINEBREAK*/
552        {   
553          goto parsed;
554        }
555     
556      if(*ch == ',')
557        {       
558          ch++;
559          goto numberPairs;
560        }
561     
562      if(*ch == '\\')
563        {
564          ch++;
565          skipWhites(&ch);
566         
567          if(!isNum(*ch))
568            {
569              printf("%c Number expected in %s\n", *ch, p_names[i]);
570              exit(-1);
571            }     
572         
573          l = 0;
574          while(isNum(*ch))
575            {
576              buf[l] = *ch;
577              ch++;     
578              l++;
579            }
580          buf[l] = '\0';
581          modulo = atoi(buf);     
582          partitions[i][3 + 3 * (pairsCount - 1)] = modulo;     
583         
584          skipWhites(&ch);
585          if(*ch == '\0' || *ch == '\n' || *ch == '\r')
586            {       
587              goto parsed;
588            }
589          if(*ch == ',')
590            {         
591              ch++;
592              goto numberPairs;
593            }
594        } 
595     
596      if(*ch == '/')
597        {
598          printf("\nRAxML detected the character \"/\" in your partition file.\n");
599          printf("Did you mean to write something similar to this: \"DNA, p1=1-100\\3\" ?\n");
600          printf("It's actually a backslash, not a slash, the program will exit now with an error!\n\n");
601        }   
602     
603      assert(0);
604       
605    parsed:
606      i = i;
607    }
608 
609  fclose(f);
610 
611  /*********************************************************************************************************************/ 
612
613  for(i = 0; i <= rdta->sites; i++)
614    tr->model[i] = -1;
615 
616  for(i = 0; i < numberOfModels; i++)
617    {   
618      as = partitions[i][0];     
619     
620      for(j = 0; j < as; j++)
621        {
622          lower = partitions[i][1 + j * 3];
623          upper = partitions[i][2 + j * 3]; 
624          modulo = partitions[i][3 + j * 3];   
625         
626          if(modulo == -1)
627            {
628              for(k = lower; k <= upper; k++)
629                setModel(i, k, tr->model);
630            }
631          else
632            {
633              for(k = lower; k <= upper; k += modulo)
634                {
635                  if(k <= rdta->sites)
636                    setModel(i, k, tr->model);       
637                }
638            }
639        }       
640    }
641
642
643  for(i = 1; i < rdta->sites + 1; i++)
644    {
645     
646      if(tr->model[i] == -1)
647        {
648          printf("ERROR: Alignment Position %d has not been assigned any model\n", i);
649          exit(-1);
650        }     
651    } 
652
653  for(i = 0; i < numberOfModels; i++)
654    {
655      rax_free(partitions[i]);
656      rax_free(p_names[i]);
657    }
658 
659  rax_free(partitions);
660  rax_free(p_names);   
661   
662  tr->NumberOfModels = numberOfModels;     
663 
664  if(adef->perGeneBranchLengths)
665    {
666      if(tr->NumberOfModels > NUM_BRANCHES)
667        {
668          printf("You are trying to use %d partitioned models for an individual per-gene branch length estimate.\n", tr->NumberOfModels);
669          printf("Currently only %d are allowed to improve efficiency.\n", NUM_BRANCHES);
670          printf("\n");
671          printf("In order to change this please replace the line \"#define NUM_BRANCHES   %d\" in file \"axml.h\" \n", NUM_BRANCHES);
672          printf("by \"#define NUM_BRANCHES   %d\" and then re-compile RAxML.\n", tr->NumberOfModels);
673          exit(-1);
674        }
675      else
676        {
677          tr->multiBranch = 1;
678          tr->numBranches = tr->NumberOfModels;
679        }
680    }
681}
682
683/*******************************************************************************************************************************/
684
685void handleExcludeFile(tree *tr, analdef *adef, rawdata *rdta)
686{
687  FILE *f; 
688  char buf[256];
689  int
690    ch,
691    j, value, i,
692    state = 0,
693    numberOfModels = 0,
694    l = -1,
695    excludeRegion   = 0,
696    excludedColumns = 0,
697    modelCounter    = 1;
698  int
699    *excludeArray, *countArray, *modelList;
700  int
701    **partitions;
702
703  printf("\n\n");
704
705  f = myfopen(excludeFileName, "rb");   
706
707  while((ch = getc(f)) != EOF)
708    {
709      if(ch == '-')
710        numberOfModels++;
711    } 
712
713  excludeArray = (int*)rax_malloc(sizeof(int) * (rdta->sites + 1));
714  countArray   = (int*)rax_malloc(sizeof(int) * (rdta->sites + 1));
715  modelList    = (int *)rax_malloc((rdta->sites + 1)* sizeof(int));
716
717  partitions = (int **)rax_malloc(sizeof(int *) * numberOfModels); 
718  for(i = 0; i < numberOfModels; i++)
719    partitions[i] = (int *)rax_malloc(sizeof(int) * 2);
720
721  rewind(f);
722 
723  while((ch = getc(f)) != EOF)
724    {     
725      switch(state)
726        {
727        case 0: /* get first number */
728          if(!whitechar(ch))
729            {
730              if(!isNum(ch))
731                {
732                  printf("exclude file must have format: number-number [number-number]*\n");
733                  exit(-1);
734                }
735              l = 0;
736              buf[l++] = ch;
737              state = 1;
738            }
739          break;
740        case 1: /*get the number or detect - */
741          if(!isNum(ch) && ch != '-')
742            {
743              printf("exclude file must have format: number-number [number-number]*\n");
744              exit(-1);
745            }
746          if(isNum(ch))
747            {
748              buf[l++] = ch;
749            }
750          else
751            {
752              buf[l++] = '\0';       
753              value = atoi(buf);
754              partitions[excludeRegion][0] = value;
755              state = 2;
756            }
757          break;
758        case 2: /*get second number */
759          if(!isNum(ch))
760            {
761              printf("exclude file must have format: number-number [number-number]*\n");
762              exit(-1);
763            }
764          l = 0;
765          buf[l++] = ch;
766          state = 3;
767          break;
768        case 3: /* continue second number or find end */         
769          if(!isNum(ch) && !whitechar(ch))
770            {
771              printf("exclude file must have format: number-number [number-number]*\n");
772              exit(-1);
773            }
774          if(isNum(ch))
775            {
776              buf[l++] = ch;
777            }
778          else
779            {         
780              buf[l++] = '\0';       
781              value = atoi(buf);
782              partitions[excludeRegion][1] = value;
783              excludeRegion++;
784              state = 0;
785            }
786          break;
787        default:
788          assert(0);
789        }
790    }
791     
792  if(state == 3)
793    {
794      buf[l++] = '\0';     
795      value = atoi(buf);
796      partitions[excludeRegion][1] = value;
797      excludeRegion++;
798    }
799 
800  assert(excludeRegion == numberOfModels);
801
802  for(i = 0; i <= rdta->sites; i++)
803    {
804      excludeArray[i] = -1;
805      countArray[i] = 0;     
806      modelList[i] = -1;
807    } 
808
809  for(i = 0; i < numberOfModels; i++)
810    {
811      int lower = partitions[i][0];
812      int upper = partitions[i][1];
813
814      if(lower > upper)
815        {
816          printf("Misspecified exclude region %d\n", i);
817          printf("lower bound %d is greater than upper bound %d\n", lower, upper);
818          exit(-1);
819        }
820
821      if(lower == 0)
822        {
823          printf("Misspecified exclude region %d\n", i);
824          printf("lower bound must be greater than 0\n");
825          exit(-1);
826        }
827
828      if(upper > rdta->sites)
829        {
830          printf("Misspecified exclude region %d\n", i);
831          printf("upper bound %d must be smaller than %d\n", upper, (rdta->sites + 1));
832          exit(-1);
833        }       
834      for(j = lower; j <= upper; j++)
835        {
836          if(excludeArray[j] != -1)
837            {
838              printf("WARNING: Exclude regions %d and %d overlap at position %d (already excluded %d times)\n", 
839                     excludeArray[j], i, j, countArray[j]);
840            }
841          excludeArray[j] = i;
842          countArray[j]   =  countArray[j] + 1; 
843        }
844    }
845
846  for(i = 1; i <= rdta->sites; i++)
847    {
848      if(excludeArray[i] != -1)
849        excludedColumns++;
850      else
851        {
852          modelList[modelCounter] = tr->model[i];
853          modelCounter++;
854        }
855    }
856
857  printf("You have excluded %d out of %d columns\n", excludedColumns, rdta->sites);
858
859  if(excludedColumns == rdta->sites)
860    {
861      printf("Error: You have excluded all sites\n");
862      exit(-1);
863    }
864
865  if(adef->useSecondaryStructure && (excludedColumns > 0))
866    {
867      char mfn[2048];
868      int countColumns;
869      FILE *newFile;
870
871      assert(adef->useMultipleModel);
872
873      strcpy(mfn, secondaryStructureFileName);
874      strcat(mfn, ".");
875      strcat(mfn, excludeFileName);
876
877      newFile = myfopen(mfn, "wb");
878
879      printBothOpen("\nA secondary structure file with analogous structure assignments for non-excluded columns is printed to file %s\n", mfn);                     
880                 
881      for(i = 1, countColumns = 0; i <= rdta->sites; i++)
882        {                 
883          if(excludeArray[i] == -1)
884            fprintf(newFile, "%c", tr->secondaryStructureInput[i - 1]);
885          else
886            countColumns++;
887        }
888                 
889      assert(countColumns == excludedColumns);
890                 
891      fprintf(newFile,"\n");
892                 
893      fclose(newFile);
894    }
895
896
897  if(adef->useMultipleModel && (excludedColumns > 0))
898    {     
899      char mfn[2048];
900      FILE *newFile;
901
902      strcpy(mfn, modelFileName);
903      strcat(mfn, ".");
904      strcat(mfn, excludeFileName);
905
906      newFile = myfopen(mfn, "wb");
907
908      printf("\nA partition file with analogous model assignments for non-excluded columns is printed to file %s\n", mfn);     
909             
910      for(i = 0; i < tr->NumberOfModels; i++)
911        {
912          boolean modelStillExists = FALSE;
913                 
914          for(j = 1; (j <= rdta->sites) && (!modelStillExists); j++)
915            {
916              if(modelList[j] == i)
917                modelStillExists = TRUE;
918            }
919
920          if(modelStillExists)
921            {                 
922              int k = 1;
923              int lower, upper;
924              int parts = 0;
925
926              switch(tr->partitionData[i].dataType)
927                {
928                case AA_DATA:                               
929                  {
930                    char AAmodel[1024];
931                   
932                    strcpy(AAmodel, protModels[tr->partitionData[i].protModels]);
933                    if(tr->partitionData[i].usePredefinedProtFreqs == FALSE)
934                      strcat(AAmodel, "F");               
935                   
936                    fprintf(newFile, "%s, ", AAmodel);
937                  }
938                  break;
939                case DNA_DATA:
940                  fprintf(newFile, "DNA, ");
941                  break;
942                case BINARY_DATA:
943                  fprintf(newFile, "BIN, ");
944                  break;
945                case GENERIC_32:
946                  fprintf(newFile, "MULTI, ");
947                  break;
948                case GENERIC_64:
949                  fprintf(newFile, "CODON, ");
950                  break;
951                default:
952                  assert(0);
953                }
954
955              fprintf(newFile, "%s = ", tr->partitionData[i].partitionName);
956             
957              while(k <= rdta->sites)
958                {
959                  if(modelList[k] == i)
960                    {
961                      lower = k;
962                      while((modelList[k + 1] == i) && (k <= rdta->sites))                                     
963                        k++;
964                      upper = k;
965                     
966                      if(lower == upper)                 
967                        {
968                          if(parts == 0)
969                            fprintf(newFile, "%d", lower);
970                          else
971                            fprintf(newFile, ",%d", lower);
972                        }
973                      else
974                        {
975                          if(parts == 0)
976                            fprintf(newFile, "%d-%d", lower, upper);
977                          else
978                            fprintf(newFile, ",%d-%d", lower, upper);
979                        }                 
980                      parts++;
981                    }
982                  k++;
983                }
984              fprintf(newFile, "\n");
985            }             
986        }       
987      fclose(newFile);
988    }
989
990 
991  {
992    FILE *newFile;
993    char mfn[2048];
994   
995
996    strcpy(mfn, seq_file);
997    strcat(mfn, ".");
998    strcat(mfn, excludeFileName);
999   
1000    newFile = myfopen(mfn, "wb");
1001   
1002    printf("\nAn alignment file with excluded columns is printed to file %s\n\n\n", mfn);
1003   
1004    fprintf(newFile, "%d %d\n", tr->mxtips, rdta->sites - excludedColumns);
1005   
1006    for(i = 1; i <= tr->mxtips; i++)
1007      {   
1008        unsigned char *tipI =  &(rdta->y[i][1]);
1009        fprintf(newFile, "%s ", tr->nameList[i]);
1010       
1011        for(j = 0; j < rdta->sites; j++)
1012          {
1013            if(excludeArray[j + 1] == -1)             
1014              fprintf(newFile, "%c", getInverseMeaning(tr->dataVector[j + 1], tipI[j]));                 
1015          }
1016       
1017        fprintf(newFile, "\n");
1018      }
1019   
1020    fclose(newFile);
1021  }
1022
1023 
1024  fclose(f);
1025  for(i = 0; i < numberOfModels; i++)
1026    rax_free(partitions[i]);
1027  rax_free(partitions); 
1028  rax_free(excludeArray);
1029  rax_free(countArray);
1030  rax_free(modelList);
1031}
1032
1033
1034void parseProteinModel(double *externalAAMatrix, char *fileName)
1035{
1036  FILE 
1037    *f = myfopen(fileName, "rb"); 
1038 
1039  int 
1040    doublesRead = 0,
1041    result = 0,
1042    i, 
1043    j;
1044 
1045  double 
1046    acc = 0.0;
1047
1048  printf("\nParsing user-defined protein model from file %s\n", fileName); 
1049
1050  while(doublesRead < 420)
1051    {     
1052      result = fscanf(f, "%lf", &(externalAAMatrix[doublesRead++]));           
1053
1054      if(result == EOF)
1055        {
1056          printf("Error protein model file must consist of exactly 420 entries \n");
1057          printf("The first 400 entries are for the rates of the AA matrix, while the\n");
1058          printf("last 20 should contain the empirical base frequencies\n");
1059          printf("Reached End of File after %d entries\n", (doublesRead - 1));
1060          exit(-1);
1061        }   
1062    }
1063       
1064  fclose(f);
1065
1066  /* CHECKS */
1067  for(i = 0; i < 20; i++)
1068    for(j = 0; j < 20; j++)
1069      {
1070        if(i != j)
1071          {
1072            if(externalAAMatrix[i * 20 + j] != externalAAMatrix[j * 20 + i])
1073              {
1074                printf("Error user-defined Protein model matrix must be symmetric\n");
1075                printf("Entry P[%d][%d]=%f at position %d is not equal to P[%d][%d]=%f at position %d\n", 
1076                       i, j,  externalAAMatrix[i * 20 + j], (i * 20 + j),
1077                       j, i,  externalAAMatrix[j * 20 + i], (j * 20 + i));
1078                exit(-1);
1079              }
1080          }
1081      }
1082
1083  acc = 0.0;
1084
1085  for(i = 400; i < 420; i++)   
1086    acc += externalAAMatrix[i];         
1087
1088  if((acc > 1.0 + 1.0E-6) || (acc <  1.0 - 1.0E-6))
1089    {
1090      printf("Base frequencies in user-defined AA substitution matrix do not sum to 1.0\n");
1091      printf("the sum is %1.80f\n", acc);
1092      exit(-1);
1093    }
1094}
1095
1096
1097
1098
1099void parseSecondaryStructure(tree *tr, analdef *adef, int sites)
1100{
1101  if(adef->useSecondaryStructure)
1102    {
1103      FILE *f = myfopen(secondaryStructureFileName, "rb");
1104
1105      int
1106        i,
1107        k,
1108        countCharacters = 0,
1109        ch,
1110        *characters,
1111        **brackets,
1112        opening,
1113        closing,
1114        depth,
1115        numberOfSymbols,
1116        numSecondaryColumns;     
1117
1118      unsigned char bracketTypes[4][2] = {{'(', ')'}, {'<', '>'},  {'[', ']'},  {'{', '}'}};         
1119
1120      numberOfSymbols = 4;     
1121
1122      tr->secondaryStructureInput = (char*)rax_malloc(sizeof(char) * sites);
1123
1124      while((ch = fgetc(f)) != EOF)
1125        {
1126          if(ch == '(' || ch == ')' || ch == '<' || ch == '>' || ch == '[' || ch == ']' || ch == '{' || ch == '}' || ch == '.')
1127            countCharacters++;
1128          else
1129            {
1130              if(!whitechar(ch))
1131                {
1132                  printf("Secondary Structure file %s contains character %c at position %d\n", secondaryStructureFileName, ch, countCharacters + 1);
1133                  printf("Allowed Characters are \"( ) < > [ ] { } \" and \".\" \n");
1134                  errorExit(-1);
1135                }
1136            }
1137        }
1138     
1139      if(countCharacters != sites)
1140        {
1141          printf("Error: Alignment length is: %d, secondary structure file has length %d\n", sites, countCharacters);
1142          errorExit(-1);
1143        }
1144   
1145      characters = (int*)rax_malloc(sizeof(int) * countCharacters); 
1146
1147      brackets = (int **)rax_malloc(sizeof(int*) * numberOfSymbols);
1148     
1149      for(k = 0; k < numberOfSymbols; k++)       
1150        brackets[k]   = (int*)rax_calloc(countCharacters, sizeof(int));
1151
1152      rewind(f);
1153
1154      countCharacters = 0;
1155      while((ch = fgetc(f)) != EOF)
1156        {
1157          if(!whitechar(ch)) 
1158            {
1159              tr->secondaryStructureInput[countCharacters] = ch;
1160              characters[countCharacters++] = ch;
1161            }
1162        }
1163     
1164      assert(countCharacters == sites);
1165
1166      for(k = 0; k < numberOfSymbols; k++)
1167        {
1168          for(i = 0, opening = 0, closing = 0, depth = 0; i < countCharacters; i++)
1169            {
1170              if((characters[i] == bracketTypes[k][0] || characters[i] == bracketTypes[k][1]) && 
1171                 (tr->extendedDataVector[i+1] == AA_DATA || tr->extendedDataVector[i+1] == BINARY_DATA ||
1172                  tr->extendedDataVector[i+1] == GENERIC_32 || tr->extendedDataVector[i+1] == GENERIC_64))
1173                {
1174                  printf("Secondary Structure only for DNA character positions \n");
1175                  printf("I am at position %d of the secondary structure file and this is not part of a DNA partition\n", i+1);
1176                  errorExit(-1);
1177                }
1178             
1179              if(characters[i] == bracketTypes[k][0])
1180                {             
1181                  depth++;
1182                  /*printf("%d %d\n", depth, i);*/
1183                  brackets[k][i] = depth;
1184                  opening++;
1185                }
1186              if(characters[i] == bracketTypes[k][1])
1187                {         
1188                  brackets[k][i] = depth; 
1189                  /*printf("%d %d\n", depth, i);  */
1190                  depth--;     
1191                 
1192                  closing++;
1193                }                         
1194             
1195              if(closing > opening)
1196                {
1197                  printf("at position %d there is a closing bracket too much\n", i+1);
1198                  errorExit(-1);
1199                }
1200            }   
1201
1202          if(depth != 0)
1203            {
1204              printf("Problem: Depth: %d\n", depth);
1205              printf("Your secondary structure file may be missing a closing or opening paraenthesis!\n");
1206            }
1207          assert(depth == 0);
1208         
1209          if(countCharacters != sites)
1210            {
1211              printf("Problem: sec chars: %d sites: %d\n",countCharacters, sites);
1212              printf("The number of sites in the alignment does not match the length of the secondary structure file\n");
1213            }
1214          assert(countCharacters == sites);
1215       
1216     
1217          if(closing != opening)
1218            {
1219              printf("Number of opening brackets %d should be equal to number of closing brackets %d\n", opening, closing);
1220              errorExit(-1);
1221            }
1222        }
1223     
1224      for(i = 0, numSecondaryColumns = 0; i < countCharacters; i++)
1225        {
1226          int checkSum = 0;
1227
1228          for(k = 0; k < numberOfSymbols; k++)
1229            {
1230              if(brackets[k][i] > 0)
1231                {
1232                  checkSum++;
1233                 
1234                  switch(tr->secondaryStructureModel)
1235                    {
1236                    case SEC_16:
1237                    case SEC_16_A:
1238                    case SEC_16_B:
1239                    case SEC_16_C:
1240                    case SEC_16_D:
1241                    case SEC_16_E:
1242                    case SEC_16_F:
1243                    case SEC_16_I:
1244                    case SEC_16_J:
1245                    case SEC_16_K:
1246                      tr->extendedDataVector[i+1] = SECONDARY_DATA;
1247                      break;
1248                    case SEC_6_A:
1249                    case SEC_6_B:
1250                    case SEC_6_C:
1251                    case SEC_6_D:
1252                    case SEC_6_E:
1253                      tr->extendedDataVector[i+1] = SECONDARY_DATA_6;
1254                      break;
1255                    case SEC_7_A:
1256                    case SEC_7_B:
1257                    case SEC_7_C:
1258                    case SEC_7_D:
1259                    case SEC_7_E:
1260                    case SEC_7_F:
1261                      tr->extendedDataVector[i+1] = SECONDARY_DATA_7;
1262                      break;
1263                    default:
1264                      assert(0);
1265                    }
1266                 
1267                  numSecondaryColumns++;
1268                }
1269            }
1270          assert(checkSum <= 1);
1271        }
1272     
1273      assert(numSecondaryColumns % 2 == 0);
1274     
1275      /*printf("Number of secondary columns: %d merged columns: %d\n", numSecondaryColumns, numSecondaryColumns / 2);*/
1276
1277      tr->numberOfSecondaryColumns = numSecondaryColumns;
1278      if(numSecondaryColumns > 0)
1279        {
1280          int model = tr->NumberOfModels;
1281          int countPairs;
1282          pInfo *partBuffer = (pInfo*)rax_malloc(sizeof(pInfo) * tr->NumberOfModels);
1283
1284          for(i = 1; i <= sites; i++)
1285            {
1286              for(k = 0; k < numberOfSymbols; k++)
1287                {
1288                  if(brackets[k][i-1] > 0)
1289                    tr->model[i] = model;
1290                }
1291
1292            }
1293
1294          /* now make a copy of partition data */
1295
1296         
1297          for(i = 0; i < tr->NumberOfModels; i++)
1298            {
1299              partBuffer[i].partitionName = (char*)rax_malloc((strlen(tr->extendedPartitionData[i].partitionName) + 1) * sizeof(char));
1300              strcpy(partBuffer[i].partitionName, tr->extendedPartitionData[i].partitionName);
1301              strcpy(partBuffer[i].proteinSubstitutionFileName, tr->extendedPartitionData[i].proteinSubstitutionFileName);
1302              partBuffer[i].dataType =  tr->extendedPartitionData[i].dataType;
1303              partBuffer[i].protModels=  tr->extendedPartitionData[i].protModels;
1304              partBuffer[i].usePredefinedProtFreqs=  tr->extendedPartitionData[i].usePredefinedProtFreqs;             
1305            }
1306
1307          for(i = 0; i < tr->NumberOfModels; i++)
1308            rax_free(tr->extendedPartitionData[i].partitionName);
1309          rax_free(tr->extendedPartitionData);
1310         
1311          tr->extendedPartitionData = (pInfo*)rax_malloc(sizeof(pInfo) * (tr->NumberOfModels + 1));
1312         
1313          for(i = 0; i < tr->NumberOfModels; i++)
1314            {
1315              tr->extendedPartitionData[i].partitionName = (char*)rax_malloc((strlen(partBuffer[i].partitionName) + 1) * sizeof(char));
1316              strcpy(tr->extendedPartitionData[i].partitionName, partBuffer[i].partitionName);
1317              strcpy(tr->extendedPartitionData[i].proteinSubstitutionFileName, partBuffer[i].proteinSubstitutionFileName);
1318              tr->extendedPartitionData[i].dataType =  partBuffer[i].dataType;
1319              tr->extendedPartitionData[i].protModels= partBuffer[i].protModels;
1320              tr->extendedPartitionData[i].usePredefinedProtFreqs=  partBuffer[i].usePredefinedProtFreqs;             
1321              rax_free(partBuffer[i].partitionName);
1322            }
1323          rax_free(partBuffer);
1324
1325          tr->extendedPartitionData[i].partitionName = (char*)rax_malloc(64 * sizeof(char));
1326
1327          switch(tr->secondaryStructureModel)
1328            {
1329            case SEC_16:
1330            case SEC_16_A:
1331            case SEC_16_B:
1332            case SEC_16_C:
1333            case SEC_16_D:
1334            case SEC_16_E:
1335            case SEC_16_F:
1336            case SEC_16_I:
1337            case SEC_16_J:
1338            case SEC_16_K:
1339              strcpy(tr->extendedPartitionData[i].partitionName, "SECONDARY STRUCTURE 16 STATE MODEL");
1340              tr->extendedPartitionData[i].dataType = SECONDARY_DATA;
1341              break;
1342            case SEC_6_A:
1343            case SEC_6_B:
1344            case SEC_6_C:
1345            case SEC_6_D:
1346            case SEC_6_E:
1347              strcpy(tr->extendedPartitionData[i].partitionName, "SECONDARY STRUCTURE 6 STATE MODEL");
1348              tr->extendedPartitionData[i].dataType = SECONDARY_DATA_6;
1349              break;
1350            case SEC_7_A:
1351            case SEC_7_B:
1352            case SEC_7_C:
1353            case SEC_7_D:
1354            case SEC_7_E:
1355            case SEC_7_F:
1356              strcpy(tr->extendedPartitionData[i].partitionName, "SECONDARY STRUCTURE 7 STATE MODEL");
1357              tr->extendedPartitionData[i].dataType = SECONDARY_DATA_7;
1358              break;
1359            default:
1360              assert(0);
1361            }
1362
1363          tr->extendedPartitionData[i].protModels= -1;
1364          tr->extendedPartitionData[i].usePredefinedProtFreqs=  FALSE;   
1365
1366          tr->NumberOfModels++; 
1367         
1368          if(adef->perGeneBranchLengths)
1369            {
1370              if(tr->NumberOfModels > NUM_BRANCHES)
1371                {
1372                  printf("You are trying to use %d partitioned models for an individual per-gene branch length estimate.\n", tr->NumberOfModels);
1373                  printf("Currently only %d are allowed to improve efficiency.\n", NUM_BRANCHES);
1374                  printf("Note that the number of partitions has automatically been incremented by one to accomodate secondary structure models\n");
1375                  printf("\n");
1376                  printf("In order to change this please replace the line \"#define NUM_BRANCHES   %d\" in file \"axml.h\" \n", NUM_BRANCHES);
1377                  printf("by \"#define NUM_BRANCHES   %d\" and then re-compile RAxML.\n", tr->NumberOfModels);
1378                  exit(-1);
1379                }
1380              else
1381                {
1382                  tr->multiBranch = 1;
1383                  tr->numBranches = tr->NumberOfModels;
1384                }
1385            }
1386         
1387          assert(countCharacters == sites);
1388
1389          tr->secondaryStructurePairs = (int*)rax_malloc(sizeof(int) * countCharacters);
1390          for(i = 0; i < countCharacters; i++)
1391            tr->secondaryStructurePairs[i] = -1;
1392          /*
1393            for(i = 0; i < countCharacters; i++)
1394            printf("%d", brackets[i]);
1395            printf("\n");
1396          */
1397          countPairs = 0;
1398
1399          for(k = 0; k < numberOfSymbols; k++)
1400            {
1401              i = 0;
1402             
1403         
1404              while(i < countCharacters)
1405                {
1406                  int 
1407                    j = i,
1408                    bracket = -1,
1409                    openBracket,
1410                    closeBracket;
1411                 
1412                  while(j < countCharacters && ((bracket = brackets[k][j]) == 0))
1413                    {
1414                      i++;
1415                      j++;
1416                    }
1417
1418                  assert(bracket >= 0);
1419
1420                  if(j == countCharacters)
1421                    {
1422                      assert(bracket == 0);
1423                      break;
1424                    }
1425           
1426                  openBracket = j;
1427                  j++;
1428                 
1429                  while(bracket != brackets[k][j] && j < countCharacters)
1430                    j++;
1431                  assert(j < countCharacters);
1432                  closeBracket = j;
1433                 
1434                  assert(closeBracket < countCharacters && openBracket < countCharacters);
1435
1436                  assert(brackets[k][closeBracket] > 0 && brackets[k][openBracket] > 0);
1437                 
1438                  /*printf("%d %d %d\n", openBracket, closeBracket, bracket);*/
1439                  brackets[k][closeBracket] = 0;
1440                  brackets[k][openBracket]  = 0;                     
1441                  countPairs++;
1442                 
1443                  tr->secondaryStructurePairs[closeBracket] = openBracket;
1444                  tr->secondaryStructurePairs[openBracket] = closeBracket;
1445                }
1446
1447              assert(i == countCharacters);
1448            }
1449             
1450          assert(countPairs == numSecondaryColumns / 2);
1451
1452         
1453          /*for(i = 0; i < countCharacters; i++)
1454            printf("%d ", tr->secondaryStructurePairs[i]);
1455            printf("\n");*/
1456         
1457
1458          adef->useMultipleModel = TRUE;
1459
1460        }
1461
1462     
1463      for(k = 0; k < numberOfSymbols; k++)       
1464        rax_free(brackets[k]);
1465      rax_free(brackets);
1466      rax_free(characters);
1467
1468      fclose(f);
1469    }
1470}
Note: See TracBrowser for help on using the repository browser.