TNT_Script/guoyi.run

macro=;

log %7.tnt.log;

/*primary notification*/
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   |
|          - filename extension must be lower case   |
|            fasta (fas/fasta), tnt (tnt/ss), nexus  |
|            (nex/nexus) format are allowed          |
|          - 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+jak+boot+sym=7.1 (default)      |
|          - prefix can be empty, or a string        |
|            default is empty                        |
\----------------------------------------------------/;
  proc/; 
end

/*Var for handling args*/
var:
  dojak doboot dosym dobrs dorbrs rsmp dobremer
  dostr domjr dohlf contype
  doew doiw doeiw dowt wtstring wttype
  kvalue
  search
  isfas istnt isnex input;

set dojak 0; set doboot 0; set dosym 0; set dobrs 0; set dorbrs 0; set rsmp 0; set dobremer 0;
set dostr 0; set domjr 0; set dohlf 0;
set doew 0; set doiw 0; set doeiw 0; set dowt 0; set wttype 0;
set kvalue 12;
set search 0;
set isfas 0; set istnt 0; set isnex 0;

goto=%0;

/*Get prefix and extension of input*/
set input $%1;

/*confirm input format*/
if ((eqstring [ $input>. fas ]) || (eqstring [ $input>. fasta ]))
  set isfas 1;
else
if ((eqstring [ $input>. tnt ]) || (eqstring [ $input>. ss ]))
  set istnt 1;
else
if ((eqstring [ $input>. nex ]) || (eqstring [ $input>. nexus ]))
  set isnex 1;	
else
  errmsg extension name of input file must be fas/fasta (for fas format), tnt/ss (for Hennig86/NONA/TNT format), nex/nexus (for nex format);
end end end

/*handle weighting type*/
if (argnumber>=3)
  set wtstring $%3;
  if (eqstring [ $wtstring ew ])
    set doew 1;
    set wttype 1;
  else
  if (eqstring [ $wtstring iw ])
    set doiw 1;
    set wttype 2;
  else
  if (eqstring [ $wtstring eiw ])
    set doeiw 1;
    set wttype 3;
  else
   goto wttypeerr;
  end end end
else
  set doeiw 1;
  set wttype 3;
end

if ('doiw' || 'doeiw')
  set dowt 1;
end

/*handle k value*/
if (argnumber>=4)
  set kvalue %4;
  if ( 'doew' && ('kvalue' != 0))
    goto kvalueerr;
  end
end

/*handle resmaple type*/
if (argnumber<6)
  set dojak 1; set doboot 1; set dosym 1; set dorbrs 1;
else
  set rsmp %6;
  loop 1 5
    if ('rsmp' >= 4)
      set dosym 1;
      set rsmp 'rsmp'-4;
    else
    if ('rsmp' >= 2)
      set doboot 1;
      set rsmp 'rsmp'-2;
    else
    if ('rsmp' >= 1)
      set dojak 1;
      set rsmp 'rsmp'-1;
    else
    if ('rsmp' >= 0.2)
      set dobrs 1;
      set rsmp 'rsmp'-0.2;
    else
    if ('rsmp' >= 0.1)
      set dorbrs 1;
      set rsmp 'rsmp'-0.1;
    end end end end end
  stop
end

if ('rsmp' != 0)
  errmsg Input 'rsmp' is illegal;
end

if ('dobrs' || 'dorbrs')
  set dobremer 1;
end

/*handle consensus type*/
if (argnumber>=6)
  set contype $%5;
  if ('dobremer')
    set dostr 1;
    if (eqstring [ $contype str ])
    else
    if (eqstring [ $contype mjr ])
       goto bremererr;
    else
    if (eqstring [ $contype hlf ])
      goto bremererr;
    else
      goto contypeerr;
    end end end
  else
    if (eqstring [ $contype str ])
      set dostr 1;
    else
    if (eqstring [ $contype mjr ])
      set domjr 1;
    else
    if (eqstring [ $contype hlf ])
      set dohlf 1;
    else
      goto contypeerr;
    end end end
  end
else
if (argnumber==5)
  set contype $%5;
  if ((eqstring [ $contype mjr ]) || (eqstring [ $contype hlf ]))
    goto bremererr;
  else
  if (eqstring [ $contype str ])
    set dostr 1;
  else
    goto contypeerr;
  end end
else
  set dostr 1;
end end

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

/*Set K*/ 
if ('dowt')
  piwe='kvalue';
end

/*Reopen tnt*/
if ('istnt' || 'isnex')
  procedure $input;
else
if ('isfas')
  procedure & $input;
end end
hold 10000;

/*handle search*/
if (ntax<=24) /*taxa lower than 25 (including 25)*/
  if ('doeiw')
    set search 2;
  else
    set search 1;
  end
else
  if (ntax<74) /*taxa lower than 75 do mult*/
    set search 2;
  else
    set search 3;
  end
end

/*Report what will be done*/
quote 
/-----------------------------------------------\;
if ('wttype'==1)
  quote
|    Equal weighting will be used.               |;
else
if ('wttype'==2)
  quote
|    Implied weighting will be used, K is 'kvalue'.    |;
else
if ('wttype'==3)
  quote
|    Ex-implied weighting will be used, K is 'kvalue'. |;
end end end

if ('search'==1)
  quote
|    Implicit enumeration will be performed.     |;
else
if ('search'==2)
  quote
|    TBR Mult will be performed.                 |;
else
if ('search'==3)
  quote
|    Xmult will be performed.                    |;
end end end

if ('dostr')
  quote
|    Strict consensus will be used.              |;
else
if ('domjr')
  quote
|    Majority-rule consensus will be used.       |;
else
if ('dohlf')
  quote
|    Half strict consensus will be used.         |;
end end end

quote
|    resample.svg will contain a tree with       |;
if ('dorbrs')
  quote
|        relative bremer support                 |;
end
if ('dobrs')
  quote
|        bremer support                          |;
end
if ('dojak')
  quote
|        jackknifing                             |;
end
if ('doboot')
  quote
|        bootstrap                               |;
end
if ('dosym')
  quote
|        symmetric resampling                    |;
end

quote
|    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 ('doiw')
  piwe&;
else
if ('doeiw')
  xpiwe(*;
  log %7.eiw.log;
  piwe&;
  log/;
  log + %7.tnt.log;
end end

/*Search trees*/
if ('search'==1)
  ienum;
else
if ('search'==2)
  mult=replic 1000 tbr hold 10;
  bbreak=tbr fill;        
else
if ('search'==3)
  sect: slack 40; 
  xmult=hit 50 replications 20 drift 10 ratchet 10 fuse 10 hold 1 keepall;
  bbreak=tbr fill;
end 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 consensus tree*/
if ('dostr')
  nelsen * 0.'npars';
else
if ('domjr')
  majority * 0.'npars';
else
if ('dohlf')
  comcomp * 0.'npars';
end end end

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

/*Store consensus tree to tree vault*/
hold /+0; 
tv>/;
tchoose 0.'npars';

/*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,  					*/

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

if ('dobremer')

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

  /*calculate relative bremer support*/
  if ('dorbrs')
  bs ]!! 0;
  end

  /*calculate bremer support*/
  if ('dobrs')
  macfloat 1; /*set the br value float*/
  bsupport !! 0;
  end

end
  
/*Choose final tree*/
tv<;
tchoose/;

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

/*jakknifing*/
if ('dojak')
  resample jak replications 1000 from 0; /*from 0 will orphan other trees*/
end

/*bootstrap*/
if ('doboot')
  resample boot replications 1000 from 0;
end

/*symmetric resampling*/
if ('dosym')
  resample sym replications 1000 from 0;
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*/
if ('themax' != 'themin')
	set RI ('themax'-'TL')/('themax'-'themin'); /*RI=1 character fits perfetcly*/
else
	set RI 999;
end

/*Report CI/RI/TL */
log %7.report.log;
macfloat 3;
quote Consistency Index (CI) is 'CI';
if ('RI' == 999)
	quote Retention Index (RI) is NA;
else
	quote Retention Index (RI) is 'RI';
end
quote Tree Length (TL) is 'TL';
log/;

watch-;
/*Generate ss file for winclada*/
log %7.winclada.ss;
/*get basic info*/
var: taxnum chanum;
set taxnum ntax+1;
set chanum nchar+1;
/*output head*/
quote
xread
'
Data saved from TNT
'
'chanum' 'taxnum';
/*outhead body*/
xread!;
/*output split symnol*/
runc!
tntprintf(";\n\n");
!
/*output the tread*/
taxname-;
tp*;
/*output the end*/
runc!
tntprintf("proc/;\n");
!
log/;
log + %7.tnt.log;

/*Report*/
quote 
/----------------------------------------------\
|      The analysis has been finished.         |
|      The file `tnt.log` contains             |
|          the performing details              |
|      The file `trees*.tre` contain           |
|          trees found by mult and xmult       |
|      The file `original*.tre` contain        |
|          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 tre file without taxname     |
|      The file `resample.svg` contain         |;

if ('dostr')
  quote
|          strict consensus tree with          |;
else
if ('domjr')
  quote
|          majority-rule consensus tree with   |;
else
if ('dohlf')
  quote
|          half strict consensus tree with     |;
end end end

if ('dorbrs')
  quote
|          relative bremer support             |;
end
if ('dobrs')
  quote    
|          bremer support                      |;
end
if ('dojak')
  quote
|          jackknifing                         |;
end
if ('doboot')
  quote
|          bootstrap                           |;
end
if ('dosym')
  quote
|          symmetric resampling                |;
end

if ('doew')
  quote
|          under equal weighting               |;
else
if ('doiw')
  quote
|          under implied weighting,            |
|          K value is 'kvalue'                       |;
else
if ('doeiw')
  quote
|          under extended implied weighting,   |
|          K value is 'kvalue'                       |;
end end 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` contains    |
|          the tree tags text                  |
|      The file `homo.log` contains the report |
|          of character homoplasy              |;

if ('doeiw')
quote
|      The file `eiw.log` contains the report  |
|          of character concavities            |;
end

quote
|      The file `winclada.tre` can be          |
|          converted by tnt2winclada           |
|      The file `winclada.ss` can be read by   |
|          winclada directly                   |
\----------------------------------------------/;

/*Quit*/
zzz;

/*consensus type error*/
label contypeerr
errmsg Consensus type must be str mjr or hlf;
proc/;

/*weighting type error*/
label wttypeerr
errmsg Weighting type must be ew, iw or eiw;
proc/;

/*equal weighting k value error*/
label kvalueerr
errmsg Equal weighting (ew) must be followed with 0 instead of 'kvalue';
proc/;

/*bremer support not strict consensus error*/
label bremererr
errmsg Bremer support or any bremer support variations must use strict consensus;
proc/;