.gitignore

+*.a
+*.o
+*.help
+
+bin/
+doc/man/
+

debian/changelog

-phast (1.4+dfsg-2) UNRELEASED; urgency=medium
+phast (1.5+dfsg-1) unstable; urgency=medium
 
+  * New upstream version
   * Fix homepage
   * debhelper 11
   * Point Vcs fields to salsa.debian.org
   * Standards-Version: 4.2.1
   * Run part or the test suite
-  TODO: Homepage has version 1.5, Github is lagging behind.
-        Download at homepage is lacking test suite
 
- -- Andreas Tille <tille@debian.org>  Thu, 18 Oct 2018 18:16:40 +0200
+ -- Andreas Tille <tille@debian.org>  Fri, 19 Oct 2018 09:10:34 +0200
 
 phast (1.4+dfsg-1) unstable; urgency=medium
 

include/external_libs.h

@@ -33,7 +33,7 @@
 #else
 
 #ifdef VECLIB
-#include <vecLib/clapack.h>
+#include <Accelerate/Accelerate.h>
 #define LAPACK_INT __CLPK_integer
 #define LAPACK_DOUBLE __CLPK_doublereal
 #else

include/misc.h

@@ -87,6 +87,17 @@ int int_pow(int x, int y) {
   return retval;
 }
 
+/** Take x,y & max as int and cast x & y as double
+ *  * raise x to the power y. return min(x^y,max) as int*/
+
+static PHAST_INLINE
+int pow_bounded(int x, int y, int max) {
+  double z = pow((double) x, (double) y);
+  return ( z > (double) max ? max : (int)z);
+}
+
+
+
 /** \name Log based calculation functions
 \{ */
 

src/lib/base/matrix.c

@@ -293,7 +293,8 @@ int mat_invert(Matrix *M_inv, Matrix *M) {
    same dimension, and C is diagonal.  C is described by a vector
    representing its diagonal elements.  */
 void mat_mult_diag(Matrix *A, Matrix *B, Vector *C, Matrix *D) {
-  int i, j, k;
+
+/*  int i, j, k;
   for (i = 0; i < C->size; i++) {
     for (j = 0; j < C->size; j++) {
       A->data[i][j] = 0;
@@ -302,7 +303,27 @@ void mat_mult_diag(Matrix *A, Matrix *B, Vector *C, Matrix *D) {
     }
   }
 }
+*/
+
+A->data[0][0] = B->data[0][0] * C->data[0] * D->data[0][0] + B->data[0][1] * C->data[1] * D->data[1][0] + B->data[0][2] * C->data[2] * D->data[2][0] + B->data[0][3] * C->data[3] * D->data[3][0];
+A->data[0][1] = B->data[0][0] * C->data[0] * D->data[0][1] + B->data[0][1] * C->data[1] * D->data[1][1] + B->data[0][2] * C->data[2] * D->data[2][1] + B->data[0][3] * C->data[3] * D->data[3][1];
+A->data[0][2] = B->data[0][0] * C->data[0] * D->data[0][2] + B->data[0][1] * C->data[1] * D->data[1][2] + B->data[0][2] * C->data[2] * D->data[2][2] + B->data[0][3] * C->data[3] * D->data[3][2];
+A->data[0][3] = B->data[0][0] * C->data[0] * D->data[0][3] + B->data[0][1] * C->data[1] * D->data[1][3] + B->data[0][2] * C->data[2] * D->data[2][3] + B->data[0][3] * C->data[3] * D->data[3][3];
+A->data[1][0] = B->data[1][0] * C->data[0] * D->data[0][0] + B->data[1][1] * C->data[1] * D->data[1][0] + B->data[1][2] * C->data[2] * D->data[2][0] + B->data[1][3] * C->data[3] * D->data[3][0];
+A->data[1][1] = B->data[1][0] * C->data[0] * D->data[0][1] + B->data[1][1] * C->data[1] * D->data[1][1] + B->data[1][2] * C->data[2] * D->data[2][1] + B->data[1][3] * C->data[3] * D->data[3][1];
+A->data[1][2] = B->data[1][0] * C->data[0] * D->data[0][2] + B->data[1][1] * C->data[1] * D->data[1][2] + B->data[1][2] * C->data[2] * D->data[2][2] + B->data[1][3] * C->data[3] * D->data[3][2];
+A->data[1][3] = B->data[1][0] * C->data[0] * D->data[0][3] + B->data[1][1] * C->data[1] * D->data[1][3] + B->data[1][2] * C->data[2] * D->data[2][3] + B->data[1][3] * C->data[3] * D->data[3][3];
+A->data[2][0] = B->data[2][0] * C->data[0] * D->data[0][0] + B->data[2][1] * C->data[1] * D->data[1][0] + B->data[2][2] * C->data[2] * D->data[2][0] + B->data[2][3] * C->data[3] * D->data[3][0];
+A->data[2][1] = B->data[2][0] * C->data[0] * D->data[0][1] + B->data[2][1] * C->data[1] * D->data[1][1] + B->data[2][2] * C->data[2] * D->data[2][1] + B->data[2][3] * C->data[3] * D->data[3][1];
+A->data[2][2] = B->data[2][0] * C->data[0] * D->data[0][2] + B->data[2][1] * C->data[1] * D->data[1][2] + B->data[2][2] * C->data[2] * D->data[2][2] + B->data[2][3] * C->data[3] * D->data[3][2];
+A->data[2][3] = B->data[2][0] * C->data[0] * D->data[0][3] + B->data[2][1] * C->data[1] * D->data[1][3] + B->data[2][2] * C->data[2] * D->data[2][3] + B->data[2][3] * C->data[3] * D->data[3][3];
+A->data[3][0] = B->data[3][0] * C->data[0] * D->data[0][0] + B->data[3][1] * C->data[1] * D->data[1][0] + B->data[3][2] * C->data[2] * D->data[2][0] + B->data[3][3] * C->data[3] * D->data[3][0];
+A->data[3][1] = B->data[3][0] * C->data[0] * D->data[0][1] + B->data[3][1] * C->data[1] * D->data[1][1] + B->data[3][2] * C->data[2] * D->data[2][1] + B->data[3][3] * C->data[3] * D->data[3][1];
+A->data[3][2] = B->data[3][0] * C->data[0] * D->data[0][2] + B->data[3][1] * C->data[1] * D->data[1][2] + B->data[3][2] * C->data[2] * D->data[2][2] + B->data[3][3] * C->data[3] * D->data[3][2];
+A->data[3][3] = B->data[3][0] * C->data[0] * D->data[0][3] + B->data[3][1] * C->data[1] * D->data[1][3] + B->data[3][2] * C->data[2] * D->data[2][3] + B->data[3][3] * C->data[3] * D->data[3][3];
 
+
+}
 int mat_equal(Matrix *A, Matrix *B) {
   int i, j;
   if (A->nrows != B->nrows || A->ncols != B->ncols) return 0;

src/lib/msa/maf.c

@@ -139,12 +139,10 @@ MSA *maf_read_cats_subset(FILE *F, FILE *REFSEQF, int tuple_size,
       msa->alloc_len = msa->length = refseqlen;  //this may still not be big enough because of gaps in refseq
     else 
       msa->alloc_len = 50000;
-    max_tuples =  max(1000000, (int)pow(strlen(msa->alphabet)+strlen(msa->missing)+1, 
-				   2 * msa->nseqs * tuple_size));
+    max_tuples = pow_bounded(strlen(msa->alphabet)+strlen(msa->missing)+1, 2 * msa->nseqs * tuple_size, 1000000);
   }
   else 
-    max_tuples = min(1000000,
-		     (int)pow(strlen(msa->alphabet)+strlen(msa->missing)+1, 2 * msa->nseqs * tuple_size));
+    max_tuples = pow_bounded(strlen(msa->alphabet)+strlen(msa->missing)+1, 2 * msa->nseqs * tuple_size, 1000000);
 
   if (max_tuples > 10000000 || max_tuples < 0) max_tuples = 10000000;
   if (max_tuples < 1000000) max_tuples = 1000000;
@@ -596,15 +594,13 @@ MSA *maf_read_unsorted(FILE *F, FILE *REFSEQF, int tuple_size, char *alphabet,
 
   if (store_order) {
     msa->length = map != NULL ? map->msa_len : refseqlen;
-    max_tuples = min(msa->length,
-                     (int)pow(strlen(msa->alphabet)+strlen(msa->missing)+1, msa->nseqs * tuple_size));
+    max_tuples = pow_bounded(strlen(msa->alphabet)+strlen(msa->missing)+1, msa->nseqs * tuple_size, msa->length);
     if (max_tuples < 0) max_tuples = msa->length;
     if (max_tuples > 1000000) max_tuples = 1000000; 
   }
   else {
     msa->length = 0;
-    max_tuples = min(50000,
-                     (int)pow(strlen(msa->alphabet)+strlen(msa->missing)+1, msa->nseqs * tuple_size));
+    max_tuples = pow_bounded(strlen(msa->alphabet)+strlen(msa->missing)+1, msa->nseqs * tuple_size, 50000);
     if (max_tuples < 0) max_tuples = 50000;
   }
 

src/lib/msa/maf_block.c

@@ -13,7 +13,7 @@
 #include <maf_block.h>
 #include <hashtable.h>
 #include <ctype.h>
-
+#include <assert.h>
 
 MafBlock *mafBlock_new() {
   MafBlock *block = smalloc(sizeof(MafBlock));
@@ -201,7 +201,10 @@ void mafBlock_add_iLine(String *line, MafSubBlock *sub) {
   
   for (i=0; i<6; i++) str_free((String*)lst_get_ptr(l, i));
   lst_free(l);
+  if (sub->numLine >= 4) die("Error: bad MAF file");
   sub->lineType[sub->numLine++] = 'i';
+
+
 }
 
 
@@ -247,7 +250,9 @@ void mafBlock_add_qLine(String *line, MafSubBlock *sub) {
    
   for (i=0; i<2; i++) str_free((String*)lst_get_ptr(l, i));
   lst_free(l);
+  if (sub->numLine >= 4) die("Error: bad MAF file");
   sub->lineType[sub->numLine++] = 'q';
+
 }
 
 
@@ -710,7 +715,8 @@ void mafSubBlock_strip_iLine(MafSubBlock *sub) {
   for (i=0; i<sub->numLine; i++) 
     if (sub->lineType[i]=='i') break;
   if (i < sub->numLine) {
-    for (j=i+1; j<sub->numLine; j++) {
+    assert(i < 4);
+    for (j=i+1; j<sub->numLine && j<4; j++) {
       sub->lineType[j-1] = sub->lineType[j];
       if (sub->lineType[j] == 'i')
 	die("ERROR mafSubBlock_strip_iLine: sub->lineType[%i]=%c\n",

src/lib/msa/sufficient_stats.c

@@ -117,8 +117,8 @@ void ss_from_msas(MSA *msa, int tuple_size, int store_order,
       upper_bound = source_msa->ss->ntuples;
     else if (source_msa != NULL) upper_bound = source_msa->length;
     else upper_bound = min(msa->length, MAX_NTUPLE_ALLOC);
-    max_tuples = min((int)pow(strlen(msa->alphabet)+ (int)strlen(msa->missing) + 1, 
-                         msa->nseqs * tuple_size),
+    max_tuples = pow_bounded(strlen(msa->alphabet)+ (int)strlen(msa->missing) + 1,
+                         msa->nseqs * tuple_size,   
                      upper_bound);
     if (max_tuples < 0) max_tuples = upper_bound;  //protect against overflow
     ss_new(msa, tuple_size, max_tuples, do_cats, store_order);
@@ -139,8 +139,8 @@ void ss_from_msas(MSA *msa, int tuple_size, int store_order,
       upper_bound = msa->ss->ntuples + source_msa->ss->ntuples;
     else
       upper_bound = msa->ss->ntuples + source_msa->length;
-    max_tuples = min((int)pow(strlen(msa->alphabet) + (int)strlen(msa->missing) + 1, 
-                         msa->nseqs * tuple_size),
+      max_tuples = pow_bounded(strlen(msa->alphabet) + (int)strlen(msa->missing) + 1,
+                         msa->nseqs * tuple_size,
                      upper_bound);
     if (max_tuples < 0) max_tuples = upper_bound;
     ss_realloc(msa, tuple_size, max_tuples, do_cats, store_order);

src/make-include.mk

@@ -133,7 +133,7 @@ endif
 # vecLib
 ifdef VECLIB
 CFLAGS += -DVECLIB
-LIBS = -lphast -framework vecLib -lc -lm
+LIBS = -lphast -framework Accelerate -lc -lm
 
 # CLAPACK
 else

src/util/clean_genes.c

@@ -34,16 +34,16 @@
 typedef enum {OKAY, BAD_REF, BAD_REF_START, BAD_REF_STOP, BAD_REF_5_SPLICE,
               BAD_REF_3_SPLICE, BAD_REF_ORF, BAD_REF_INDEL_STRICT_FAIL,
               NO_ALN, BAD_START, BAD_STOP, BAD_5_SPLICE, BAD_3_SPLICE, 
-              BAD_5_SPLICE_UTR, BAD_3_SPLICE_UTR, NONSENSE, FSHIFT, 
-              BAD_INTRON, TOO_MANY_Ns, WARN_FSHIFT, WARN_Ns, NTYPES} 
+              BAD_5_SPLICE_UTR, BAD_3_SPLICE_UTR, NONSENSE, FRAMESHIFT, 
+              BAD_INTRON, TOO_MANY_Ns, WARN_FRAMESHIFT, WARN_Ns, NTYPES} 
 status_type;                    /* NTYPES marks the total number */
 
-/* possible gap types for cds exon -- frame-shift gaps (FSHIFT_BAD),
-   no apparent frame shift (FSHIFT_OK), all gaps multiple of 3 in len
+/* possible gap types for cds exon -- frame-shift gaps (FRAMESHIFT_BAD),
+   no apparent frame shift (FRAMESHIFT_OK), all gaps multiple of 3 in len
    (CLN_GAPS), gaps mult. of 3 and nonoverlapping (NOVRLP_CLN_GAPS),
    and no gaps present (NGAPS).  Defined in order of least to most
    stringent criterion (ordering gets used) */
-typedef enum {FSHIFT_BAD, FSHIFT_OK, CLN_GAPS, 
+typedef enum {FRAMESHIFT_BAD, FRAMESHIFT_OK, CLN_GAPS, 
               NOVRLP_CLN_GAPS, NGAPS, NGAP_TYPES} 
 cds_gap_type;                   /* NGAP_TYPES marks the total number */
 
@@ -102,7 +102,7 @@ typedef struct Problem {
   status_type status;
   int start;
   int end;
-  cds_gap_type cds_gap;  /* if status if FSHIFT */
+  cds_gap_type cds_gap;  /* if status if FRAMESHIFT */
 } Problem;
 
 /* create a new problem.  feat can be null for whole gene.
@@ -311,7 +311,7 @@ int gap_compare(const void *ptr1, const void* ptr2) {
   return g1->start - g2->start;
 }
 /* scans a cds for gaps.  Returns CLN_GAPS, NOVRLP_CLN_GAPS, NO_GAPS,
-   or FSHIFT_BAD; doesn't try to check for compensatory indels, which
+   or FRAMESHIFT_BAD; doesn't try to check for compensatory indels, which
    is more complicated (this is left for the special-purpose function
    below) */
 int scan_for_gaps(GFF_Feature *feat, MSA *msa, Problem **problem) {
@@ -322,7 +322,7 @@ int scan_for_gaps(GFF_Feature *feat, MSA *msa, Problem **problem) {
   cds_gap_type retval = NGAPS;
   List *gaps = lst_new_ptr(10);
 
-  for (j = 0; retval != FSHIFT_BAD && j < msa->nseqs; j++) {
+  for (j = 0; retval != FRAMESHIFT_BAD && j < msa->nseqs; j++) {
     for (i = msa_start; i <= msa_end; i++) {
       if (ss_get_char_pos(msa, i, j, 0) == GAP_CHAR) {
         int gap_start, gap_end;
@@ -336,9 +336,9 @@ int scan_for_gaps(GFF_Feature *feat, MSA *msa, Problem **problem) {
         gap_end--;              /* inclusive */
 
         if ((gap_end - gap_start + 1) % 3 != 0) {
-          retval = FSHIFT_BAD;
-          *problem = problem_new(feat, FSHIFT, gap_start, gap_end);
-          (*problem)->cds_gap = FSHIFT_BAD;
+          retval = FRAMESHIFT_BAD;
+          *problem = problem_new(feat, FRAMESHIFT, gap_start, gap_end);
+          (*problem)->cds_gap = FRAMESHIFT_BAD;
           break;
         }
 
@@ -449,18 +449,18 @@ int is_fshift_okay(GFF_Feature *feat, MSA *msa) {
    "clean" gaps (all multiples of 3 in length; CLEAN_GAPS)
    non-overlapping clean gaps (NOVRLP_CLN_GAPS), "okay" gaps (only
    temporary frame shifts, corrected by compensatory indels;
-   FSHIFT_OK), or real frame-shift gaps (FSHIFT_BAD) */
+   FRAMESHIFT_OK), or real frame-shift gaps (FRAMESHIFT_BAD) */
 cds_gap_type get_cds_gap_type(GFF_Feature *feat, MSA *msa, List *problems) {
   Problem *problem = NULL;
   cds_gap_type retval = scan_for_gaps(feat, msa, &problem);
 
-  if (retval == FSHIFT_BAD && is_fshift_okay(feat, msa)) {
-    retval = FSHIFT_OK;
+  if (retval == FRAMESHIFT_BAD && is_fshift_okay(feat, msa)) {
+    retval = FRAMESHIFT_OK;
                                 /* most of the time the call to
                                    is_fshift_okay won't be
                                    necessary */
-    problem->status = WARN_FSHIFT;
-    problem->cds_gap = FSHIFT_OK;
+    problem->status = WARN_FRAMESHIFT;
+    problem->cds_gap = FRAMESHIFT_OK;
   }
   if (problem != NULL) {
     lst_push_ptr(problems, problem);
@@ -669,9 +669,9 @@ void write_log(FILE *logf, GFF_FeatureGroup *group, status_type status,
     case NONSENSE:
       reason = "Nonsense mutation";
       break;
-    case FSHIFT:
+    case FRAMESHIFT:
       {
-        if (problem->cds_gap ==  FSHIFT_OK) 
+        if (problem->cds_gap ==  FRAMESHIFT_OK) 
           reason = "Frame-shift gap [gaps not clean]";
         else if (problem->cds_gap == CLN_GAPS) 
           reason = "Frame-shift gap [gaps clean but overlapping/near boundary]";
@@ -689,7 +689,7 @@ void write_log(FILE *logf, GFF_FeatureGroup *group, status_type status,
       reason = "Too many Ns";
       print_alignment = FALSE;
       break;
-    case WARN_FSHIFT:
+    case WARN_FRAMESHIFT:
       reason = "Frame shift gap (with compensatory gap)";
       break;
     case WARN_Ns:
@@ -741,13 +741,13 @@ char *status_type_str(status_type status) {
     return "bad_3_splice_utr";
   case NONSENSE:
     return "nonsense";
-  case FSHIFT:
+  case FRAMESHIFT:
     return "frameshift";
   case BAD_INTRON:
     return "bad_intron";
   case TOO_MANY_Ns:
     return "too_many_Ns";
-  case WARN_FSHIFT:
+  case WARN_FRAMESHIFT:
     return "warning_fshift";
   case WARN_Ns:
     return "warning_Ns";
@@ -980,7 +980,7 @@ int main(int argc, char *argv[]) {
     *discards=NULL, *intron_splice = lst_new_ptr(10);
   char *rseq_fname = NULL;
   FILE *logf = NULL, *mlogf = NULL, *statsf = NULL, *discardf = NULL;
-  cds_gap_type fshift_mode = FSHIFT_BAD;
+  cds_gap_type fshift_mode = FRAMESHIFT_BAD;
   char *groupby = "transcript_id";
   msa_coord_map *map;
   int *countNs, *countCDSs;
@@ -1029,11 +1029,11 @@ int main(int argc, char *argv[]) {
       break;
     case 'f':
       check_alignment = 1;
-      fshift_mode = FSHIFT_OK;
+      fshift_mode = FRAMESHIFT_OK;
       break;
     case 'c':
       check_alignment = check_start = check_stop = check_splice = check_nonsense = 1;
-      if (fshift_mode < FSHIFT_OK) fshift_mode = FSHIFT_OK;
+      if (fshift_mode < FRAMESHIFT_OK) fshift_mode = FRAMESHIFT_OK;
       break;
     case 'N':
       Nfrac = get_arg_dbl_bounds(optarg, 0, 1);
@@ -1166,7 +1166,7 @@ int main(int argc, char *argv[]) {
     List *gfeatures = group->features;
     GFF_Feature *feat;
     status_type status = OKAY;
-    cds_gap_type gt = FSHIFT_BAD;
+    cds_gap_type gt = FRAMESHIFT_BAD;
     problems_clear(problems);
 
     /* make sure have frame info for CDSs */
@@ -1295,9 +1295,9 @@ int main(int argc, char *argv[]) {
 
           else if (str_equals_charstr(feat->feature, GFF_CDS_TYPE)) {
  
-            if (fshift_mode > FSHIFT_BAD 
+            if (fshift_mode > FRAMESHIFT_BAD 
 		&& (gt = get_cds_gap_type(feat, msa, problems)) < fshift_mode) {
-              if (status == OKAY || status == NONSENSE) status = FSHIFT;
+              if (status == OKAY || status == NONSENSE) status = FRAMESHIFT;
             }
 
             if (check_nonsense && !is_nonsense_clean(feat, msa, problems)) {
@@ -1336,8 +1336,8 @@ int main(int argc, char *argv[]) {
             for (j = 0; j < lst_size(problems); j++) {
               struct Problem *problem = lst_get_ptr(problems, j);
               status_type ftype = problem->status;
-              if ((ftype == FSHIFT || ftype == NONSENSE) && 
-                  status != FSHIFT && status != NONSENSE)
+              if ((ftype == FRAMESHIFT || ftype == NONSENSE) && 
+                  status != FRAMESHIFT && status != NONSENSE)
                 continue;       /* don't count secondary frame shifts
                                    and nonsense mutations */ 
 
@@ -1350,7 +1350,7 @@ int main(int argc, char *argv[]) {
 
           /* also keep track of the total number of "conserved exons", and
              the number having each kind of gap */
-          if ((status == OKAY || (status == FSHIFT && gt >= FSHIFT_OK))) {
+          if ((status == OKAY || (status == FRAMESHIFT && gt >= FRAMESHIFT_OK))) {
             nconserved_exons++;
             nce_gap_type[gt]++;     /* number of conserved exons having
                                        given type of gaps */
@@ -1429,10 +1429,10 @@ int main(int argc, char *argv[]) {
             nconsid[BAD_5_SPLICE], nfail[BAD_3_SPLICE], nconsid[BAD_3_SPLICE],
             nfail[BAD_5_SPLICE_UTR], nconsid[BAD_5_SPLICE_UTR],
             nfail[BAD_3_SPLICE_UTR], nconsid[BAD_3_SPLICE_UTR], 
-            nfail[BAD_INTRON], nfail[NONSENSE], nfail[FSHIFT], 
+            nfail[BAD_INTRON], nfail[NONSENSE], nfail[FRAMESHIFT], 
             nfail[TOO_MANY_Ns], nconserved_exons, nce_gap_type[NGAPS], 
             nce_gap_type[NOVRLP_CLN_GAPS], nce_gap_type[CLN_GAPS], 
-            nce_gap_type[FSHIFT_OK]);
+            nce_gap_type[FRAMESHIFT_OK]);
     fprintf(statsf, "%s", STATS_DESCRIPTION);
   }
 

src/util/msa_view.c

@@ -627,7 +627,7 @@ int main(int argc, char* argv[]) {
       if (f->frame == GFF_NULL_FRAME) f->frame = 0;
       if (fourD_refseq == NULL) fourD_refseq = str_new_charstr(f->seqname->chars);
       else if (!str_equals(fourD_refseq, f->seqname))
-	die("--4d requires all features have same source column");
+	die("--4d requires all features are on the same chromosome");
       if (str_equals_charstr(f->feature, "CDS") && f->strand != '-')
 	str_cpy_charstr(f->feature, "CDSplus");
       else if (str_equals_charstr(f->feature, "CDS") && f->strand == '-')

test/test_phast.sh

@@ -37,39 +37,39 @@ tree_doctor --scale 3.0 hpmr.mod > hpmr_fast.mod
 ******************** phyloP ********************
 
 phyloFit hmrc.ss --tree "(human, (mouse,rat), cow)" --quiet
-@phyloP --null 10 phyloFit.mod
+@phyloP  --seed 123 --null 10 phyloFit.mod
 msa_view -o SS --end 100 hmrc.ss > hmrc_short.ss
-@phyloP phyloFit.mod hmrc_short.ss
-@phyloP --method LRT --base-by-base phyloFit.mod hmrc.ss
-@phyloP --method LRT --mode CONACC --wig-scores phyloFit.mod hmrc.ss
-@phyloP -d 12345 --method SCORE --mode NNEUT --base-by-base phyloFit.mod hmrc.ss
-@phyloP --method GERP --mode ACC --wig-scores phyloFit.mod hmrc.ss
-@phyloP --method GERP --base-by-base phyloFit.mod hmrc.ss
-@phyloP -d 12345 --method SCORE --wig-scores phyloFit.mod hmrc.ss
-@phyloP --method LRT --wig-scores phyloFit.mod hmrc.ss
-@phyloP --method LRT --base-by-base phyloFit.mod hmrc.ss
-@phyloP -d 12345 --method SCORE --wig-scores --refidx 2 phyloFit.mod hmrc.ss
+@phyloP  --seed 123 phyloFit.mod hmrc_short.ss
+@phyloP  --seed 123 --method LRT --base-by-base phyloFit.mod hmrc.ss
+@phyloP  --seed 123 --method LRT --mode CONACC --wig-scores phyloFit.mod hmrc.ss
+@phyloP  --seed 123 -d 12345 --method SCORE --mode NNEUT --base-by-base phyloFit.mod hmrc.ss
+@phyloP  --seed 123 --method GERP --mode ACC --wig-scores phyloFit.mod hmrc.ss
+@phyloP  --seed 123 --method GERP --base-by-base phyloFit.mod hmrc.ss
+@phyloP  --seed 123 -d 12345 --method SCORE --wig-scores phyloFit.mod hmrc.ss
+@phyloP  --seed 123 --method LRT --wig-scores phyloFit.mod hmrc.ss
+@phyloP  --seed 123 --method LRT --base-by-base phyloFit.mod hmrc.ss
+@phyloP  --seed 123 -d 12345 --method SCORE --wig-scores --refidx 2 phyloFit.mod hmrc.ss
 echo -e "chr1\t0\t10\nchr1\t50\t100\nchr1\t200\t300" > temp.bed
-@phyloP --method LRT --mode CONACC --features temp.bed phyloFit.mod hmrc.ss
-@phyloP --method SCORE --features temp.bed -g phyloFit.mod hmrc.ss
+@phyloP  --seed 123 --method LRT --mode CONACC --features temp.bed phyloFit.mod hmrc.ss
+@phyloP  --seed 123 --method SCORE --features temp.bed -g phyloFit.mod hmrc.ss
 tree_doctor --name-ancestors phyloFit.mod > phyloFit-named.mod
-@phyloP --method LRT --mode CONACC --subtree mouse-rat --base-by-base phyloFit-named.mod hmrc.ss
-@phyloP --method LRT --mode ACC --branch mouse-rat -w phyloFit-named.mod hmrc.ss
-@phyloP --posterior phyloFit.mod hmrc_short.ss
-@phyloP --fit-model -w --subtree mouse-rat phyloFit-named.mod hmrc_short.ss
-@phyloP --epsilon 0.0001 phyloFit-named.mod hmrc_short.ss
-@phyloP --confidence-interval 0.05 phyloFit-named.mod hmrc_short.ss
-@phyloP --quantiles --null 100 phyloFit-named.mod
+@phyloP  --seed 123 --method LRT --mode CONACC --subtree mouse-rat --base-by-base phyloFit-named.mod hmrc.ss
+@phyloP  --seed 123 --method LRT --mode ACC --branch mouse-rat -w phyloFit-named.mod hmrc.ss
+@phyloP  --seed 123 --posterior phyloFit.mod hmrc_short.ss
+@phyloP  --seed 123 --fit-model -w --subtree mouse-rat phyloFit-named.mod hmrc_short.ss
+@phyloP  --seed 123 --epsilon 0.0001 phyloFit-named.mod hmrc_short.ss
+@phyloP  --seed 123 --confidence-interval 0.05 phyloFit-named.mod hmrc_short.ss
+@phyloP  --seed 123 --quantiles --null 100 phyloFit-named.mod
 
 
 # want to try subtree with LRT, SCORE, SPH, with both --base-by-base and --features
 # note that bugs in SCORE method mean tests will fail against revisions < 6563
-@phyloP --method LRT --subtree mouse-rat --mode CONACC --base-by-base phyloFit-named.mod hmrc.ss
-@phyloP --method LRT --subtree mouse-rat --mode CONACC --features temp.bed phyloFit-named.mod hmrc.ss
-@phyloP --method SCORE --subtree mouse-rat --mode CONACC --base-by-base phyloFit-named.mod hmrc.ss
-@phyloP --method SCORE --subtree mouse-rat --mode CONACC --features temp.bed phyloFit-named.mod hmrc.ss
-@phyloP --method SPH --subtree mouse-rat --mode CONACC --base-by-base phyloFit-named.mod hmrc.ss
-@phyloP --method SPH --subtree mouse-rat --mode CONACC --features temp.bed phyloFit-named.mod hmrc.ss
+@phyloP  --seed 123 --method LRT --subtree mouse-rat --mode CONACC --base-by-base phyloFit-named.mod hmrc.ss
+@phyloP  --seed 123 --method LRT --subtree mouse-rat --mode CONACC --features temp.bed phyloFit-named.mod hmrc.ss
+@phyloP  --seed 123 --method SCORE --subtree mouse-rat --mode CONACC --base-by-base phyloFit-named.mod hmrc.ss
+@phyloP  --seed 123 --method SCORE --subtree mouse-rat --mode CONACC --features temp.bed phyloFit-named.mod hmrc.ss
+@phyloP  --seed 123 --method SPH --subtree mouse-rat --mode CONACC --base-by-base phyloFit-named.mod hmrc.ss
+@phyloP  --seed 123 --method SPH --subtree mouse-rat --mode CONACC --features temp.bed phyloFit-named.mod hmrc.ss
 
 rm -f hmrc_short.ss phyloFit.mod phyloFit-named.mod temp.bed
 

version

-v1.3
+v1.4