Problem with setting fixed topology for MrBayes

Question

I need to use partitioned DNA input (concatenated multiple genes) and fixed topology (species tree from ASTRAL) to get branch lengths in MrBayes. I have input NEXUS like:

#NEXUS
BEGIN TAXA;
  DIMENSIONS NTAX = 318;
  TAXLABELS
    Achillea_nobilis_SAMN11585351
    Anthemis_arvensis_SAMN11585352
    ...
    Haplopappus_sp_SAMN11585375
  ;
END
BEGIN TREES;
  TRANSLATE
    1       Achillea_nobilis_SAMN11585351,
    2       Anthemis_arvensis_SAMN11585352,
    ...
    318     Haplopappus_sp_SAMN11585375
  ;
  TREE treeshrink_exons_good_nobrl = [&R] (1,((2,3...);
END;
BEGIN DATA;
dimensions ntax=318 nchar=169024;
format datatype=dna missing=? gap=-;
MATRIX
Achillea_nobilis_SAMN11585351               TTATGCAA...
Anthemis_arvensis_SAMN11585352              TTATGCAA...
...
Haplopappus_sp_SAMN11585375                 TTATGCAA...
;
END;
BEGIN MRBAYES;
  charset At1g09700 = 1-213;
  charset At1g10740 = 214-456;
  ...
  charset At5g66280 = 168530-168757;
  charset At5g66680 = 168758-169024;
  outgroup Achillea_nobilis_SAMN11585351;
  partition by_gene=467:At1g09700,At1g10740,...,At5g66280,At5g66680;
  set partition=by_gene;
  prset topologypr=fixed(treeshrink_exons_good_nobrl);
  prset applyto=(all) ratepr=variable;
  lset applyto=(all) nst=mixed rates=gamma ngammacat=10;
  unlink tratio=(all) revmat=(all) omega=(all) statefreq=(all) shape=(all) pinvar=(all) correlation=(all) ratemultiplier=(all);
  mcmcp ngen=10000000 nruns=4 nchains=4 samplefreq=5000;
  ss;
  sump relburnin=yes burninfrac=0.25 nruns=4;
  sumt relburnin=yes burninfrac=0.25 nruns=4;
END;

I'm not 100% sure about the MRBAYES block, but generally I think it's correct. When I launch it (MrBayes 3.2.7a) I get an error which I'm unable to resolve:

$ mb sp_treeshrink_exons_good.mb.nex
...
   Executing file "sp_treeshrink_exons_good.mb.nex"
   UNIX line termination
   Longest line length = 169069
   Parsing file
   Expecting NEXUS formatted file
   Reading taxa block
      Allocated taxon set
      Defining new set of 318 taxa
   Exiting taxa block
   Reading trees block
      Successfully read tree 'treeshrink_exons_good_nobrl' # NOTE
   Exiting trees block
   Deleting user trees
   Deleting previously defined taxa
   Reading data block
      Allocated taxon set
      Allocated matrix
      Defining new matrix with 318 taxa and 169024 characters
      Data is Dna
      Missing data coded as ?
      Gaps coded as -
      Taxon   1 -> Achillea_nobilis_SAMN11585351
      Taxon   2 -> Anthemis_arvensis_SAMN11585352
...
      Taxon 318 -> Haplopappus_sp_SAMN11585375
      Successfully read matrix
      Setting default partition (does not divide up characters)
      Setting model defaults
      Seed (for generating default start values) = 1654610114
      Setting output file names to "sp_treeshrink_exons_good.mb.nex.run<i>.<p|t>"
   Exiting data block
   Reading mrbayes block
      Defining charset called 'At1g09700'
      Defining charset called 'At1g10740'
...
      Defining charset called 'At5g66280'
      Defining charset called 'At5g66680'
      Defining partition called 'by_gene'
      Setting by_gene as the partition, dividing characters into 467 parts.
      Setting model defaults
      Seed (for generating default start values) = 1901076251
      Setting Topologypr to Fixed for partition 1
      Setting Topologypr to Fixed for partition 2
...
      Setting Topologypr to Fixed for partition 466
      Setting Topologypr to Fixed for partition 467
      Could not find tree 'treeshrink_exons_good_nobrl' # NOTE
      Could not set fixed topology from user tree 'treeshrink_exons_good_nobrl' # NOTE
      Error when setting parameter "Topologypr" (2)
      The error occurred when reading char. 26-52 on line 1445
         in the file 'sp_treeshrink_exons_good.mb.nex'
...

So it correctly reads tree treeshrink_exons_good_nobrl and after while it complains it couldn't find tree treeshrink_exons_good_nobrl...? I can't find where this weird error comes from.

Answer 1

The answer is pretty simple. What you are requesting is fixing the tree prior at a tree you are specifying. The getter/setter is complaining it cannot fit the data to the specified tree because of 26 bp at 26-52 on line 1445 of the nexus file.

The first thing needed is to post output to:

sed -n '1445p' ./sp_treeshrink_exons_good.mb.nex

What I suspect you will find are non-nucleotide characters associated with "missing data" for example xxxxx, or 0000 what the nexus file wants is ????? for missing data. It is possible part of the degenerate nucleotide code is not being observed (long while since I used MB), e.g. RYB I think its unlikely.

The second approach, if the output of the sed command is pretty normal is to remove the prior setting the tree, i.e. the tree you are supplying is violating e.g. Jukes-Cantor correction.

Finally and a related issue is there's an alignment error, i.e. the alignment at this region is out of frame for these residues and Jukes-Cantor is violated so it refuses to proceed. This would be really unusual however, because I've used MB where there are clear pockets of saturation and it hasn't complained.

Personally, I'm impressed by the level of foresight in MB under written into this getter/setter (its checking your input). This algorithm is rarely used but is clearly well written.

---

Okay there is no line 1445. What I recommend is placing the alignment in either Sequtron or Clustal X and scrolling down at the position 26-52 and if not successful 'optical analysis' of the alignment.

Personally, I would use Biopython AlignIO for this task and Counter from the Collections library to interrogate the bug. However, thats coding stuff and optical analysis does the same thing. Essentially, you loop through each sequence and Counter will make a frequency count of everything in each sequence. You then say if you are not AGCT- tell me what taxa this is. If you code in Python thats a good follow-up question (separate question).