add: Bremer support

This commit is contained in:
kuoi 2023-02-27 01:24:22 +08:00
parent 53bfd967b2
commit 0a1821cb11
4 changed files with 174 additions and 117 deletions

View File

@ -1,30 +1,36 @@
# TNT Script used by Guoyi
These scripts follow MIT except for setk.run which belongs to Salvador Arias (Instituto Miguel Lillo, San Miguel de Tucumán, Argentina).
These scripts follow MIT, part of script is from setk.run belonging to Salvador Arias (Instituto Miguel Lillo, San Miguel de Tucumán, Argentina).
## Usage
### Method1
- Place this script and your tnt matrix under the folder that you call `tnt` or `exe` file is placed
- OS: Arch Linux;
- Enter `tnt`
- `tnt run guoyi.run file.tnt`, `file.tnt` is your data;
- Enter command `guoyi filename;`
- Type `;` and enter;
## Functions
### Method2
- Estimated implied weighting K value.
- Replace `%1` with your filename;
- Search trees via TBR Mult and Xmult.
- Run `tnt p guoyi.run` directly;
- Perform Strict consensus.
- Calculate Relative Bremer support, jackknifing, and bootstrap.
- Calculate TL, CI, and RI.
## Options
- `trees.tre`, `bt.tre` are trees with taxaname;
- Results instructions are at the end of `tnt.log`.
- `trees_no.tre`, `bt_no.tre` are trees without taxaname which can be put into Winclada with the `file.tnt`;
- `con.tre`, `con.tre` are trees with taxaname.
- `majority` can be replaced by `nelsen`;
- `trees_no.tre`, `bt_no.tre` are trees without taxaname which can be put into Winclada with the `file.tnt`.
- `boot` can be replaced by `jack`;
- `nelsen` can be replaced by `majority`.
- `xmult` and `mult` replications and hold trees number can be adjusted.

167
guoyi.run
View File

@ -5,29 +5,174 @@ mxram 10240;
nstates 32;
nstates NOGAPS;
procedure %1;
setk.run;
/*Arguments*/
if ( argnumber != 1 )
silent -console;
quote - /-----------------------------------------------\;
quote - | GUOYI TNT SCRIPT 2022-2023 MIT |;
quote - | You need to give your filename |;
quote - | tnt*> guoyi filename; |;
quote - \-----------------------------------------------/;
proc/;
end
/*Report what will be done*/
quote - /-----------------------------------------------\;
quote - | Implied weighting will be estimated. |;
quote - | TBR Mult and Xmult will be performed. |;
quote - | Strict consensus will be used. |;
quote - | Relative bremer support, jackknifing and |;
quote - | bootstrap will be shown on the svg. |;
quote - | TL, CI and RI will be calculated finally. |;
quote - \------------------------------------------------/;
/*Set K*/
procedure %1;
hold 1000;
var:
actK
minK
maxK
gVal
a
b
Nval
maxIts
;
if (ntax<0)
quote NO HAY DATOS!;
proc/;
end
set gVal 1;
/** busca el valor de g **/
loop 0 nchar
if (!isinfo[#1]) continue; end
if (!isact[#1]) continue; end
set a=maxsteps[#1];
set b=minsteps[#1];
if (('a'-'b')>'gVal')
set gVal 'a'-'b';
end
stop
quote valor de g = 'gVal';
set actK 10;
set minK 0;
set maxK 500; /** Asume 500 como el maximo posible valor de K **/
set maxIts 0;
/** busca el mejor valor de K **/
loop 0 1
quote actual valor de k = 'actK';
set a 1-('actK'/('actK'+1));
quote 'a';
set b ('actK'/('actK'+'gVal'-1))-('actK'/('actK'+'gVal'));
quote 'b';
set Nval 'a'/'b';
quote valor de N = 'Nval';
if (('Nval'>14.8)&&('Nval'<15.2)) endloop; end /** N esta en el rango **/
if ('Nval'>15.2)
set minK 'actK';
else
set maxK 'actK';
end
set a ('maxK'-'minK')/2;
set actK 'minK' + 'a';
set maxIts ++;
if ('maxIts'==100) endloop; end /** salida de emergencia **/
setloop 0;
stop
/** Si salio de emergencia **/
if ('maxIts'==100)
quote NO SE TERMINARON LAS ITERACIONES;
quote mejor K encontrado: 'actK';
proc /;
end
quote Valor de K: 'actK' (N='Nval');
piwe='actK';
/*Reopen tnt*/
procedure %1;
hold 5000;
/*Search trees*/
xpiwe(*;
mult=replic 1000 tbr hold 1;
mult=replic 1000 tbr hold 10;
sect: slack 1000;
xmult=rss hit 1000 drift 10 ratchet 10 fuse 10;
piwe&;
/*Export trees*/
export= trees.tre;
taxname-;
export= trees_no.tre;
taxname=;
/*Consensus tree*/
bbreak=tbr;
majority *;
nelsen *;
/*Get Bremer/jak/boot support*/
ttags=;
sub 1; hold +1000; bbreak=fill;
sub 2; hold +2000; bbreak=fill;
sub 3; hold +3000; bbreak=fill;
sub 4; hold +4000; bbreak=fill;
sub 5; hold +5000; bbreak=fill;
sub 6; hold +6000; bbreak=fill;
sub 7; hold +7000; bbreak=fill;
sub 8; hold +8000; bbreak=fill;
bsupport[;
resample jak replications 1000;
resample boot replications 1000;
ttags & bt.svg thickness 7 italics fontsize 15;
export < bt.tre;
taxname-;
export - bt_no.tre;
/*Export consensus tree with supports*/
ttags & con.svg thickness 7 italics fontsize 15;
export < con.tre;
taxname-;
export - con_no.tre;
/*Caulculate TL/CI/RI score*/
length *;
stats.run;
quit ;
report-;
var-;
var =
0 themin
1 themax
+ CIs [ ( ntrees + 1 ) ]
+ RIs [ ( ntrees + 1 ) ]
+ this
;
set themin minsteps;
set themax maxsteps;
loop 0 ntrees
progress #1 (ntrees+1) Calculating indices...;
set this length[#1];
set CIs[#1] 'themin'/'this'; /*CI=1 means no homoplasy*/
set RIs[#1] ('themax'-'this')/('themax'-'themin'); /*RI=1 character fits perfetcly*/
stop
progress/;
report=;
macfloat 3;
maketable CIs Consistency index;
maketable RIs Retention index;
/*Report*/
quote - /----------------------------------------------\;
quote - | The analysis has been finished. |;
quote - | The file `tnt.log` contains |;
quote - | K, TL, CI and RI |;
quote - | The file `trees*.tre` contain |;
quote - | trees found by mult and xmult |;
quote - | The file `con*.tre` contain |;
quote - | strict consensus tree |;
quote - | The file `*_no.tre` contain |;
quote - | tree with `taxaname-` |;
quote - | The file `con.svg` contain |;
quote - | strict consensus tree with |;
quote - | relative bremer, jak and boot |;
quote - | support on the tree |;
quote - \----------------------------------------------/;
/*Quit*/
zzz;

View File

@ -1,71 +0,0 @@
macro=
var:
actK
minK
maxK
gVal
a
b
Nval
maxIts
;
if (ntax<0)
quote NO HAY DATOS!;
proc/;
end
set gVal 1;
/* busca el valor de g */
loop 0 nchar
if (!isinfo[#1]) continue; end
if (!isact[#1]) continue; end
set a=maxsteps[#1];
set b=minsteps[#1];
if (('a'-'b')>'gVal')
set gVal 'a'-'b';
end
stop
quote valor de g = 'gVal';
set actK 10;
set minK 0;
set maxK 500; /* Asume 500 como el maximo posible valor de K */
set maxIts 0;
/* busca el mejor valor de K */
loop 0 1
quote actual valor de k = 'actK';
set a 1-('actK'/('actK'+1));
quote 'a';
set b ('actK'/('actK'+'gVal'-1))-('actK'/('actK'+'gVal'));
quote 'b';
set Nval 'a'/'b';
quote valor de N = 'Nval';
if (('Nval'>14.8)&&('Nval'<15.2)) endloop; end /* N esta en el rango */
if ('Nval'>15.2)
set minK 'actK';
else
set maxK 'actK';
end
set a ('maxK'-'minK')/2;
set actK 'minK' + 'a';
set maxIts ++;
if ('maxIts'==100) endloop; end /* salida de emergencia */
setloop 0;
stop
/* Si salio de emergencia */
if ('maxIts'==100)
quote NO SE TERMINARON LAS ITERACIONES;
quote mejor K encontrado: 'actK';
proc /;
end
quote Valor de K: 'actK' (N='Nval');
piwe='actK';
proc/;

View File

@ -1,23 +0,0 @@
macro= ;
report- ;
var =
0 themin
1 themax
+ CIs [ ( ntrees + 1 ) ]
+ RIs [ ( ntrees + 1 ) ]
+ this
;
set themin minsteps ;
set themax maxsteps ;
loop 0 ntrees
progress #1 (ntrees+1) Calculating indices... ;
set this length[#1] ;
set CIs[#1] 'themin'/'this' ;
set RIs[#1] ('themax'-'this')/('themax'-'themin') ;
stop
progress/ ;
report= ;
macfloat 3 ;
maketable CIs Consistency index ;
maketable RIs Retention index ;
proc/ ;