TNT_Script/guoyi.run

macro=;

log %7.tnt.log;

/*Arguments*/
if ( (argnumber == 0) || (argnumber >7) )
  silent -console;
  quote 
/----------------------------------------------------\
|        GUOYI TNT SCRIPT 2022-2023 MIT              |
|        You need to give your filename              |
| shell> tnt run guoyi.tnt filename, (Linux & Mac)   |
| shell> tnt run guoyi.run filename(semicolon) (Win) |
|        ==============================              |
|        === Parameters     Details ===              |
|        ==============================              |
|        filename type weight K cons resamp prefix   |
|          - type should be 32, dna, prot, num       |
|            num=number, dna=DNA, prot=protein       |
|            32=max number allowed (default)         |
|          - weight should be iw, ew, eiw            |
|            iw=implied weight, ew=equal weight      |
|            eiw=extended implied weight (default)   |
|            N.B. K of ew must followed 0 (=NA)      |
|          - K is 12 (default) following Goloboff    |
|            et al. 2017 (Cladistics 34: 407–437)    |
|            it must more than 0                     |
|          - cons should be str, mjr, hlf            |
|            mjr=majority rule, hlf=half             |
|            str=strict (default)                    |
|          - resamp should be sum of what you want   |
|            jak=1, boot=2, sym=4 relative-bremer    |
|            (rbr)=0.1, bremer(br)=0.2               |
|            i.e rbr+br+jak+boot+sym=7.3 (default)   |
|          - prefix can be empty, or a string        |
|            default is empty                        |
\----------------------------------------------------/;
  proc/; 
end

/*Basic settings*/
taxname+1000;
taxname=;
mxram 10240;
if(argnumber>=2)
  nstates %2;
else
  nstates 32;
end
nstates NOGAPS;

/*Set K*/
if(argnumber>2 )
  if (eqstring [ %3 ew ])
    quote Equal weighting is used;
    if ((argnumber>=4) && (%4 != 0))
      errmsg equal weighting %4 should be fullfilled with 0;
    end	
  else
  if ((eqstring [ %3 iw ]) || (eqstring [ %3 eiw ]))
    if (argnumber>=4)
      piwe=%4;
    else
      piwe=12;
    end
  else
    errmsg %3 must be ew iw or eiw;
  end end
else
  piwe=12;
end

/*Reopen tnt*/
procedure %1;
hold 10000;

/*Report what will be done*/
quote 
/-----------------------------------------------\;
if(argnumber>2)
  if (eqstring [ %3 iw ])
    if (argnumber>3)
      quote
|    Implied weighting will be used, K is %4.    |;
    else
      quote
|    Implied weighting will be used, K is 12.    |;
    end
  else
  if (eqstring [ %3 eiw ])
    if (argnumber>3)
      quote
|    Ex-implied weighting will be used, K is %4. |;
    else
      quote
|    Ex-implied weighting will be used, K is 12. |;
    end
  else
  if (eqstring [ %3 ew ])
    quote
|    Equal weighting will be used.               |;
  else
    errmsg %3 must be ew iw or eiw;
  end end end
else
  quote
|    Ex-implied weighting will be used, K is 12. |;
end
if (ntax<=25)
  quote
|    Implicit enumeration will be performed.     |;
else
if ((ntax>25) && (ntax < 75))
  quote
|    TBR Mult will be performed.                 |;
else
if (ntax>=75)
  quote
|    Xmult will be performed.                    |;
end end end

if (argnumber<=5)
  quote
|    Strict consensus will be used.              |;
else
if (( %6 == 0.1 )||( %6 == 0.2 )||( %6 == 0.3 )||( %6==1.1 )||( %6 ==1.2 )||( %6 ==1.3 )||( %6 == 2.1 )||( %6 == 2.2 )||( %6 == 2.3 )||( %6 == 3.1 )||( %6 == 3.2 )||( %6 == 3.3 )||( %6 == 4.1 )||( %6 == 4.2 )||( %6 == 4.3 )||( %6 == 5.1 )||( %6 == 5.2 )||( %6 == 5.3 )||( %6 == 6.1 )||( %6 == 6.2 )||( %6 == 6.3 )||( %6 == 7.1 )||( %6 == 7.2 )||( %6 == 7.3 ))
  if (!(eqstring [ %5 str ]))
    errmsg Bremer support or relative bremer support must use strict consensus;
  end
  quote
|    Strict consensus will be used.              |;
else
  if (eqstring [ %5 mjr ])
    quote
|    Majority-rule consensus will be used.       |;
  else
  if (eqstring [ %5 hlf ])
    quote
|    Half strict consensus will be used.         |;
  else
  if(eqstring [ %5 str ])
    quote
|    Strict consensus will be used.              |;
  else
    errmsg %5 must be mjr hlf or str;
  end end end
end end

if (argnumber>=6)
  if ( (%6>0) && (%6<1) )
    if (%6==0.1 || %6==0.3)
      quote
|    relative bremer support                     |;
    else
    if (%6==0.2 || %6==0.3)
      quote
|    bremer support                              |;
    end end
  else
  if ((%6>=1) && (%6<2))
    quote
|    Jackknifing                                 |;
    if (%6==1.1 || %6==1.3)
      quote
|    relative bremer support                     |;
    else
    if (%6==1.2 || %6==1.3)
      quote
|    bremer support                              |;
    end end
  else
  if ((%6>=2) && (%6<3))
    quote
|    Bootstrap                                   |;
    if (%6==2.1 || %6==2.3)
      quote
|    relative bremer support                     |;
    else
    if (%6==2.2 || %6==2.3)
      quote
|    bremer support                              |;
    end end
  else
  if ((%6>=3) && (%6<4))
    quote
|    Jackknifing and bootstrap                   |;
    if (%6==3.1 || %6==3.3)
      quote
|    relative bremer support                     |;
    else
    if (%6==3.2 || %6==3.3)
      quote
|    bremer support                              |;
    end end
  else
  if ((%6>=4) && (%6<5))
    quote
|    Symmetric resampling                        |;
    if (%6==4.1 || %6==4.3)
      quote
|    relative bremer support                     |;
    else
    if (%6==4.2 || %6==4.3)
      quote
|    bremer support                              |;
    end end
  else
  if ((%6>=5) && (%6<6))
    quote
|    Jackknifing and symmetric resampling        |;
    if (%6==5.1 || %6==5.3)
      quote
|    relative bremer support                     |;
    else
    if (%6==5.2 || %6==5.3)
      quote
|    bremer support                              |;
    end end
  else
  if ((%6>=6) && (%6<7))
    quote
|    Bootstrap and symmetric resampling          |;
    if (%6==6.1 || %6==6.3)
      quote
|    relative bremer support                     |;
    else
    if (%6==6.2 || %6==6.3)
      quote
|    bremer support                              |;
    end end
  else
  if ((%6>=7) && (%6<8))
    quote
|    Jackknifing, bootstrap and symmetric        |
|    resampling                                  |;
    if (%6==7.1 || %6==7.3)
      quote
|    relative bremer support                     |;
    else
    if (%6==7.2 || %6==7.3)
      quote
|    bremer support                              |;
    end end 
  else
    errmsg Mistake %6 value;
  end end end end end end end end
end
if (argnumber<6)
  quote
|        relative bremer support, bremer support |
|        jackknifing, bootstrap, symmetric       |
|        resampling                              |;
end
quote
|        will be shown on the resample.svg.      |
|    Apomorphic characters mapping will be shown |
|        on the apo.svg and saved to apo*.tre.   |
|    TL, CI and RI will be calculated finally.   |
\------------------------------------------------/;

/*Implied weighting settings*/
if(argnumber>2)
  if (eqstring [ %3 iw ])
    piwe&;
  else
  if (eqstring [ %3 eiw ])
    xpiwe(*;
    log %7.eiw.log;
    piwe&;
    log/;
    log + %7.tnt.log;
  else
  if (eqstring [ %3 ew ])
  else
    errmsg %3 must be ew iw or eiw;
  end end end
else
  xpiwe(*;
  log %7.eiw.log;
  piwe&;
  log/;
  log+ %7.tnt.log;
end

/*Search trees*/
if(argnumber>2)
  if(eqstring [ %3 ew ] || eqstring [ %3 iw ] )
    if (ntax<=25)
      ienum;
    else
      if (ntax<=74)
        mult=replic 1000 tbr hold 10;
        bbreak=tbr fill;
      else
      if (ntax>74)
        sect: slack 40; 
        xmult=hit 50 replications 20 drift 10 ratchet 10 fuse 10 hold 1 keepall;
      end end
    end
  else
    if (ntax<=74)
      mult=replic 1000 tbr hold 10;
      bbreak=tbr fill;
    else
      sect: slack 40;
      xmult=hit 50 replications 20 drift 20 ratchet 10 fuse 10 hold 1 keepall;
    end
  end  
else
  mult=replic 1000 tbr hold 10;
  bbreak=tbr fill;
  if (ntax>74)
    sect: slack 40;
    xmult=hit 50 replications 20 drift 10 ratchet 10 fuse 10 hold 1 keepall;
  end
end

/*Export trees*/
export= %7.trees.tre;
taxname-;
export= %7.trees_no.tre;
tsave *= %7.trees.tnt.tre;
taxname=;
tsave = %7.trees.ctf;

/*Get npars number*/
var: npars;
set npars ntrees;

/*Get rbr/br/jak/boot/sym support and get consensus tree*/
ttags=;
ttags]; /*in one line*/

/*    jak=1,		boot=2,		sym=4,		*/
/*    jak+boot=3,					*/
/*    jak+sym=5,	boot+sym=6,			*/
/*    jak+boot+sym=7,  					*/

/*    bremer=0.1,	relative-bremer=0.2		*/
/*    bremer+relative-bremer=0.3			*/

if (argnumber<=5)

  /*set value of suboptimal from most parsimony tree */
  sub: 0;

  /*save strict consensus tree for bremer*/
  hold+1;
  nelsen * 0.'npars';

  /*get the consensus tree number*/
  var: contree;
  set contree ntrees;

  /*calculate relative bremer support*/
  bs ]!! 0;

  /*set the br value float*/
  macfloat 1;

  /*calculate bremer support*/
  bsupport !! 0;

  /*choose final tree*/
  tchoose/;

  /*adjust to zero following Pablo Goloboff*/
  sub 0;

  /*jakknifing*/
  resample jak replications 1000 from 0; /*from 0 will orphan other trees*/
  /*bootstrap*/
  resample boot replications 1000 from 0;
  /*symmetric resampling*/
  resample sym replications 1000 from 0;
else

  if(( %6 == 0.1 )||( %6 == 0.2 )||( %6 == 0.3 )||( %6==1.1 )||( %6 ==1.2 )||( %6 ==1.3 )||( %6 == 2.1 )||( %6 == 2.2 )||( %6 == 2.3 )||( %6 == 3.1 )||( %6 == 3.2 )||( %6 == 3.3 )||( %6 == 4.1 )||( %6 == 4.2 )||( %6 == 4.3 )||( %6 == 5.1 )||( %6 == 5.2 )||( %6 == 5.3 )||( %6 == 6.1 )||( %6 == 6.2 )||( %6 == 6.3 )||( %6 == 7.1 )||( %6 == 7.2 )||( %6 == 7.3 )) /*bremer support related can only use strict consensus*/

    /*set value of suboptimal from most parsimony tree*/
    sub: 0;

    /*save strict consensus tree for bremer*/
    hold+1;
    nelsen * 0.'npars';

    /*get the consensus tree number*/
    var: contree;
    set contree ntrees;

    if(( %6 == 0.1 )||( %6 == 0.3 )||( %6==1.1 )||( %6 ==1.3 )||( %6 == 2.1 )||( %6 == 2.3 )||( %6 == 3.1 )||( %6 == 3.3 )||( %6 == 4.1 )||( %6 == 4.3 )||( %6 == 5.1 )||( %6 == 5.3 )||( %6 == 6.1 )||( %6 == 6.3 )||( %6 == 7.1 )||( %6 == 7.3 ))
      /*calculate relative bremer support*/
      bs ]!! 0;
    end

    if(( %6 == 0.2 )||( %6 == 0.3 )||( %6==1.2 )||( %6 ==1.3 )||( %6 == 2.2 )||( %6 == 2.3 )||( %6 == 3.2 )||( %6 == 3.3 )||( %6 == 4.2 )||( %6 == 4.3 )||( %6 == 5.2 )||( %6 == 5.3 )||( %6 == 6.2 )||( %6 == 6.3 )||( %6 == 7.2 )||( %6 == 7.3 ))

      /*set the br value float*/
      macfloat 1;

      /*calculate bremer support*/
      bsupport !! 0;
    end

    /*choose final tree*/
    tchoose/;

    /*adjust to zero following Pablo Goloboff*/
    sub 0;

    /*jackknifing*/
    if((%6 == 1.1) || (%6 == 1.2) || (%6 == 1.3) || (%6 == 3.1) || (%6 == 3.2) || (%6 == 3.3) || (%6 == 5.1) || (%6 == 5.2) || (%6 == 5.3) || (%6 == 7.1) || (%6 == 7.2) || (%6 == 7.3))
      resample jak replications 1000 from 0;
    end

    /*bootstrap*/
    if((%6 == 2.1) || (%6 == 2.2) || (%6 == 2.3) || (%6 == 3.1) || (%6 == 3.2) || (%6 == 3.3) || (%6 == 6.1) || (%6 == 6.2) || (%6 == 6.3) || (%6 == 7.1) || (%6 == 7.2) || (%6 == 7.3))
      resample boot replications 1000 from 0;
    end

    /*symmetric resampling*/
    if((%6 == 4.1) || (%6 == 4.2) || (%6 == 4.3) || (%6 == 5.1) || (%6 == 5.2) || (%6 == 5.3) || (%6 == 6.1) || (%6 == 6.2) || (%6 == 6.3) || (%6 == 7.1) || (%6 == 7.2) || (%6 == 7.3))
      resample sym replications 1000 from 0;
    end

  else /*non-bremer related should be consensus firstly*/

    if (eqstring [ %5 str ])
      nelsen * 0.'npars';
    else
    if (eqstring [ %5 mjr ])
      majority * 0.'npars';
    else
    if (eqstring [ %5 hlf ])
      comcomp * 0.'npars';
    else
      errmsg %5 must be mjr hlf or str;
    end end end

    /*choose the final tree: consenus tree*/
    tchoose/;

    /*jackkinfing*/
    if((%6==1) || (%6==3) || (%6==5) || (%6==7))
      resample jak replications 1000 from 0;
    end

    /*bootstrap*/
    if((%6==2) || (%6==3) || (%6==6) || (%6==7))
      resample boot replications 1000 from 0;
    end

    /*symmetric resampling*/
    if((%6==4) || (%6==5) || (%6==6) || (%6==7))
      resample sym replications 1000 from 0;
    end

  end

  /*Tree vault store*/
  hold /+0; 
  tvault >/;

end

/*Export consensus tree with supports*/
ttags & %7.resample.svg thickness 7 italics fontsize 15;
log %7.resample.log;
quote /-------------resample tags start---------------\;
ttags/;
quote \-------------resample tags stop----------------/;
log/;
log + %7.tnt.log;
export < %7.resample.tre;
taxname-;
export - %7.resample_no.tre;
ttags-;
taxname=;

/*Export consensus tree*/
tchoose /;
export - %7.original.tre;
taxname-;
export - %7.original_no.tre;
tsave = %7.original.ctf;
tsave *= %7.original.tnt.tre;

/*Apomorphic characters*/
export = %7.winclada.tre;
taxname =;
ttags=;
apo >0;
log %7.apo.log;
quote /-----------apomorphy tags start --------------\;
ttags/;
quote \-----------apomorphy tags stop ---------------/;
log/;
log + %7.tnt.log;
ttags & %7.apo.svg thickness 7 italics fontsize 15;
export < %7.apo.tre;
taxname-;
export - %7.apo_no.tre;
ttags-;

/*Character Analysis*/
log %7.homo.log;
chomo;
cscores;
log/;
log + %7.tnt.log;

/*Caulculate TL/CI/RI score*/
report-;
var: themin themax CI RI TL ;
set themin minsteps;
set themax maxsteps;
set TL length[0];
set CI 'themin'/'TL'; /*CI=1 means no homoplasy*/
set RI ('themax'-'TL')/('themax'-'themin'); /*RI=1 character fits perfetcly*/
report=;

/*Report CI/RI/TL */
log %7.report.log;
macfloat 3;
quote Consistency Index (CI) is 'CI';
quote Retention Index (RI) is 'RI';
quote Tree Length (TL) is 'TL';
log/;
log + %7.tnt.log;

/*Report*/
quote 
/----------------------------------------------\
|      The analysis has been finished.         |
|      The file `tnt.log` contains             |
|          K, TL, CI and RI                    |
|      The file `trees*.tre` contain           |
|          trees found by mult and xmult       |
|      The file `original*.tre` contain        |
|          strict consensus tree without label |
|      The file `resample*.tre` contain        |
|          consensus tree with support         |
|      The file `apo*.tre` contain             |
|          tree with apomorphic character      |
|      The file `*_no.tre` contain             |
|          tree with `taxname-`                |
|      The file `*.ctf` tree file is           |
|          only readable for TNT               |
|      The file `*.tnt.tre` contain            |
|          is the ctf file with taxname        |
|      The file `resample.svg` contain         |;


if (argnumber<=5)
  quote
|          strict consensus tree with          |;
else
if ( (argnumber >=6) && (( %6 == 0.1 )||( %6 == 0.2 )||( %6 == 0.3 )||( %6==1.1 )||( %6 ==1.2 )||( %6 ==1.3 )||( %6 == 2.1 )||( %6 == 2.2 )||( %6 == 2.3 )||( %6 == 3.1 )||( %6 == 3.2 )||( %6 == 3.3 )||( %6 == 4.1 )||( %6 == 4.2 )||( %6 == 4.3 )||( %6 == 5.1 )||( %6 == 5.2 )||( %6 == 5.3 )||( %6 == 6.1 )||( %6 == 6.2 )||( %6 == 6.3 )||( %6 == 7.1 )||( %6 == 7.2 )||( %6 == 7.3 )) )
  quote
|          strict consensus tree with          |;
else
  if (eqstring [ %5 mjr ])
    quote
|          majority-rule consensus tree with   |;
  else
  if (eqstring [ %5 hlf ])
    quote
|          half strict consensus tree with     |;
  else
  if(eqstring [ %5 str ])
    quote
|          strict consensus tree with          |;
  else
    errmsg %5 must be mjr hlf or str;
    errmsg Bremer support or relative bremer support must use strict consensus;
  end end end
end end

if (argnumber>=6)
  if ( (%6>0) && (%6<1) )
    if (%6==0.1 || %6==0.3)
      quote
|          relative bremer support             |;
    else
    if (%6==0.2 || %6==0.3)
      quote    
|          bremer support                      |;
    end end
  else
  if ((%6>=1) && (%6<2))
    if (%6==1.1 || %6==1.3)
      quote
|          relative bremer support             |;
    else
    if (%6==1.2 || %6==1.3)
      quote
|          bremer support                      |;
    end end
      quote
|          jackknifing                         |;
  else
  if ((%6>=2) && (%6<3))
    if (%6==2.1 || %6==2.3)
      quote
|          relative bremer support             |;
    else
    if (%6==2.2 || %6==2.3)
      quote
|          bremer support                      |;
    end end
      quote
|          bootstrap                           |;
  else
  if ((%6>=3) && (%6<4))
    if (%6==3.1 || %6==3.3)
      quote
|          relative bremer support             |;
    else
    if (%6==3.2 || %6==3.3)
      quote
|          bremer support                      |;
    end end
    quote
|          jackknifing, bootstrap              |;
  else
  if ((%6>=4) && (%6<5))
    if (%6==4.1 || %6==4.3)
      quote
|          relative bremer support             |;
    else
    if (%6==4.2 || %6==4.3)
      quote
|          bremer support                      |;
    end end
    quote
|          symmetric resampling                |;
  else
  if ((%6>=5) && (%6<6))
    if (%6==5.1 || %6==5.3)
      quote
|          relative bremer support             |;
    else
    if (%6==5.2 || %6==5.3)
      quote
|          bremer support                      |;
    end end
    quote
|          jackknifing, symmetric resampling   |; 
  else
  if ((%6>=6) && (%6<7))
    if (%6==6.1 || %6==6.3)
      quote
|          relative bremer support             |;
    else
    if (%6==6.2 || %6==6.3)
      quote
|          bremer support                      |;
    end end
    quote
|          bootstrap, symmetric resampling     |;
  else
  if ((%6>=7) && (%6<8))
    if (%6==7.1 || %6==7.3)
      quote
|          relative bremer support             |;
    else
    if (%6==7.2 || %6==7.3)
      quote
|          bremer support                      |;
    end end
    quote
|          jackknifing, bootstrap, symmetric   |
|          resampling                          |;
  else
    errmsg Mistake %6 value;
  end end end end end end end end
else
  quote
|          relative bremer support, bremer     |
|          support, jackknifing, bootstrap,    |
|          symmetric resampling                |;
end

quote
|      The file `apo.svg` contains the         |
|          tree with apomorphy mapping         |
|      The file `report.log` contains the      |
|          CI RI TL publish-needed info        |
|      The file `resample/apo.log` contain the |
|          tree tags in text                   |
|      The file `homo.log` contain the report  |
|          of character homoplasy              |
|      The file `eiw.log` contain the report   |
|          of character concavities            |
|      The file `winclada.tre` can be          |
|          converted by tnt2winclada           |
\----------------------------------------------/;

/*Quit*/
zzz;