本文章包括：raxml-ng "for 循环" & "paralle" 批量建树。
正想摸索一下怎么做，发现之前已经做过了，不过做出来的结果并没有派上用场所以忘记了，我习惯给所有用过的软件都建一个文件夹，并下载 Tutorial 文件，加上 txt 写上使用方法。
Phylosuite 也可以用 FastTree 批量建树直接导入 ASTRAL，不过 FastTree 可能不如 RAxML 准确，用 Phylosuite 里面的 RAxML 也可以，但是基因数目较多，很费内存和时间。
##1. conda 安装 raxml-ng
conda install -c bioconda raxml-ng
##############查看帮助
raxml-ng -h
RAxML-NG v. 1.2.1 released on 22.12.2023 by The Exelixis Lab.
Developed by: Alexey M. Kozlov and Alexandros Stamatakis.
Contributors: Diego Darriba, Tomas Flouri, Benoit Morel, Sarah Lutteropp, Ben Bettisworth, Julia Haag, Anastasis Togkousidis.
Latest version: https://github.com/amkozlov/raxml-ng
Questions/problems/suggestions? Please visit: https://groups.google.com/forum/#!forum/raxml
System: Intel(R) Xeon(R) Gold 6240 CPU @ 2.60GHz, 72 cores, 1510 GB RAM
Usage: raxml-ng [OPTIONS]
Commands (mutually exclusive):
--help display help information
--version display version information
--evaluate evaluate the likelihood of a tree (with model+brlen optimization)
--search ML tree search (default: 10 parsimony + 10 random starting trees)
--bootstrap bootstrapping (default: use bootstopping to auto-detect #replicates)
--all all-in-one (ML search + bootstrapping)
--support compute bipartition support for a given reference tree (e.g., best ML tree)
and a set of replicate trees (e.g., from a bootstrap analysis)
--bsconverge test for bootstrapping convergence using autoMRE criterion
--bsmsa generate bootstrap replicate MSAs
--terrace check whether a tree lies on a phylogenetic terrace
--check check alignment correctness and remove empty columns/rows
--parse parse alignment, compress patterns and create binary MSA file
--start generate parsimony/random starting trees and exit
--rfdist compute pair-wise Robinson-Foulds (RF) distances between trees
--consense [ STRICT | MR | MR<n> | MRE ] build strict, majority-rule (MR) or extended MR (MRE) consensus tree (default: MR)
eg: --consense MR75 --tree bsrep.nw
--ancestral ancestral state reconstruction at all inner nodes
--sitelh print per-site log-likelihood values
Command shortcuts (mutually exclusive):
--search1 Alias for: --search --tree rand{1}
--loglh Alias for: --evaluate --opt-model off --opt-branches off --nofiles --log result
--rf Alias for: --rfdist --nofiles --log result
Input and output options:
--tree rand{N} | pars{N} | FILE starting tree: rand(om), pars(imony) or user-specified (newick file)
N = number of trees (default: rand{10},pars{10})
--msa FILE alignment file
--msa-format VALUE alignment file format: FASTA, PHYLIP, CATG or AUTO-detect (default)
--data-type VALUE data type: DNA, AA, BIN(ary) or AUTO-detect (default)
--tree-constraint FILE constraint tree
--prefix STRING prefix for output files (default: MSA file name)
--log VALUE log verbosity: ERROR,WARNING,RESULT,INFO,PROGRESS,DEBUG (default: PROGRESS)
--redo overwrite existing result files and ignore checkpoints (default: OFF)
--nofiles do not create any output files, print results to the terminal only
--precision VALUE number of decimal places to print (default: 6)
--outgroup o1,o2,..,oN comma-separated list of outgroup taxon names (it's just a drawing option!)
--site-weights FILE file with MSA column weights (positive integers only!)
General options:
--seed VALUE seed for pseudo-random number generator (default: current time)
--pat-comp on | off alignment pattern compression (default: ON)
--tip-inner on | off tip-inner case optimization (default: OFF)
--site-repeats on | off use site repeats optimization, 10%-60% faster than tip-inner (default: ON)
--threads VALUE number of parallel threads to use (default: 72)
--workers VALUE number of tree searches to run in parallel (default: 1)
--simd none | sse3 | avx | avx2 vector instruction set to use (default: auto-detect).
--rate-scalers on | off use individual CLV scalers for each rate category (default: ON for >2000 taxa)
--force [ <CHECKS> ] disable safety checks (please think twice!)
Model options:
--model <name>+G[n]+<Freqs> | FILE model specification OR partition file
--brlen linked | scaled | unlinked branch length linkage between partitions (default: scaled)
--blmin VALUE minimum branch length (default: 1e-6)
--blmax VALUE maximum branch length (default: 100)
--blopt nr_fast | nr_safe branch length optimization method (default: nr_fast)
nr_oldfast | nr_oldsafe
--opt-model on | off ML optimization of all model parameters (default: ON)
--opt-branches on | off ML optimization of all branch lengths (default: ON)
--prob-msa on | off use probabilistic alignment (works with CATG and VCF)
--lh-epsilon VALUE log-likelihood epsilon for optimization/tree search (default: 0.1)
Topology search options:
--spr-radius VALUE SPR re-insertion radius for fast iterations (default: AUTO)
--spr-cutoff VALUE | off relative LH cutoff for descending into subtrees (default: 1.0)
--lh-epsilon-triplet VALUE log-likelihood epsilon for branch length triplet optimization (default: 1000)
Bootstrapping options:
--bs-trees VALUE number of bootstraps replicates
--bs-trees autoMRE{N} use MRE-based bootstrap convergence criterion, up to N replicates (default: 1000)
--bs-trees FILE Newick file containing set of bootstrap replicate trees (with --support)
--bs-cutoff VALUE cutoff threshold for the MRE-based bootstopping criteria (default: 0.03)
--bs-metric fbp | tbe branch support metric: fbp = Felsenstein bootstrap (default), tbe = transfer distance
--bs-write-msa on | off write all bootstrap alignments (default: OFF)
EXAMPLES:
1. Perform tree inference on DNA alignment
(10 random + 10 parsimony starting trees, general time-reversible model, ML estimate of substitution rates and
nucleotide frequencies, discrete GAMMA model of rate heterogeneity with 4 categories):
./raxml-ng --msa testDNA.fa --model GTR+G
2. Perform an all-in-one analysis (ML tree search + non-parametric bootstrap)
(10 randomized parsimony starting trees, fixed empirical substitution matrix (LG),
empirical aminoacid frequencies from alignment, 8 discrete GAMMA categories,
200 bootstrap replicates):
./raxml-ng --all --msa testAA.fa --model LG+G8+F --tree pars{10} --bs-trees 200
3. Optimize branch lengths and free model parameters on a fixed topology
(using multiple partitions with proportional branch lengths)
./raxml-ng --evaluate --msa testAA.fa --model partitions.txt --tree test.tree --brlen scaled
###########################
##主要需要用到的参数
--msa（multiple sequence alignment）输入序列文件
--outgroup 外类群多个外类群以英文逗号分开
--all ML+bootstrap
--model 模型
--bs-trees 自举值bootstrap次数
--prefix 输出文件名前缀
--tree pars{10} rand{10} 搜索树次数默认20次，可自行增加，增加会明显增加运行时间。
#建一棵树
raxml-ng --msa 10IGSconcatenation.fas --model GTR+I+G4 --outgroup Polystachya1,Polystachya2 --all --bs-trees 1000 --prefix 10IGS
##nohup后台运行，不受断网影响
nohup raxml-ng --msa 10IGSconcatenation.fas --model GTR+I+G4 --outgroup Polystachya1,Polystachya2 --all --bs-trees 1000 --prefix 10IGS 2> 10IGS.raxml.log &
#查看后台运行进程：ps -u username -eo pid,comm,etime | grep raxml
##2. 批量构建基因树，将结果放入 Trees 文件夹，得到的树文件的后缀是.support文件
#进入基因序列所在位置
cd /home/XXX/Astral/mycds
#新建一个Trees文件夹放结果
mkdir Trees
# for 循环使用 raxml-ng 给每个基因序列建树
for i in *.fas; do raxml-ng --msa $i --model GTR+I+G4 --all --bs-trees 1000 -prefix ./Trees/${i%.fas}; done
##注意检查一下是不是所有的基因都跑出树来了，每个基因有10个结果文件。我跑的时候遇到了报错：
terminate called after throwing an instance of 'std::runtime_error' what(): ERROR in branch length optimization (LIBPLL-2240): BL opt converged to a worse likelihood score by -0.912976208965802 units Aborted (core dumped)
解决方法：1. 更新raxml-ng 2. 在命令中加上 --blopt nr_safe，参考https://github.com/amkozlov/raxml-ng/issues/86
##3. 用于 ASTRAL 建物种树，需要将所有的树合并为一个文件
cat *.support > sp_genes_tree.tre
##※注意：ASTRAL 得到的都是无根树，需为 Newick 格式，输入的树需要是无根树，输出的树也是无根数，所以要用于 ASTRAL 并联法建树的话，上一个批量建树那里就不要设置外类群。
#更新####################################
#看到别人用parallel批量运行任务，于是去了解了一下
# https://hackmd.io/@mossmatters/ry9PP6_2u#IQTree
#conda安装 parallel
conda install conda-forge::parallel
#官方使用说明书 https://www.gnu.org/software/parallel/parallel_tutorial.html
#查看帮助
parallel -h
Usage:
parallel [options] [command [arguments]] < list_of_arguments
parallel [options] [command [arguments]] (::: arguments|:::: argfile(s))...
cat ... | parallel --pipe [options] [command [arguments]]
-j n Run n jobs in parallel
-k Keep same order
-X Multiple arguments with context replace
--colsep regexp Split input on regexp for positional replacements
{} {.} {/} {/.} {#} {%} {= perl code =} Replacement strings
{3} {3.} {3/} {3/.} {=3 perl code =} Positional replacement strings
With --plus: {} = {+/}/{/} = {.}.{+.} = {+/}/{/.}.{+.} = {..}.{+..} =
{+/}/{/..}.{+..} = {...}.{+...} = {+/}/{/...}.{+...}
-S sshlogin Example: foo@server.example.com
--slf .. Use ~/.parallel/sshloginfile as the list of sshlogins
--trc {}.bar Shorthand for --transfer --return {}.bar --cleanup
--onall Run the given command with argument on all sshlogins
--nonall Run the given command with no arguments on all sshlogins
--pipe Split stdin (standard input) to multiple jobs.
--recend str Record end separator for --pipe.
--recstart str Record start separator for --pipe.
GNU Parallel can do much more. See 'man parallel' for details
Academic tradition requires you to cite works you base your article on.
If you use programs that use GNU Parallel to process data for an article in a
scientific publication, please cite:
Tange, O. (2024, July 22). GNU Parallel 20240722 ('Assange').
Zenodo. https://doi.org/10.5281/zenodo.12789352
This helps funding further development; AND IT WON'T COST YOU A CENT.
If you pay 10000 EUR you should feel free to use GNU Parallel without citing.
#批量运行，首先需要一个文件list
ls -1 /home/XXX/Astral/mycds > faslist.txt
#检查一下list把不用的删掉。
ulimit -n 4096 # 增加打开文件数限制
mkdir output #新建一个文件夹，结果放入里面
parallel -j 40 "raxml-ng --msa {} --model GTR+I+G4 --all --bs-trees 1000 -prefix output/{/.}" :::: faslist.txt
#-j 并行运行任务数目，我第一次跑的时候报错了:
ERROR: Resource temporarily unavailable
似乎是同时启动太多的进程，检查并增加打开文件限制、内存CPU限制，也可以降低 -j
ulimit -a #查看文件限制
core file size (blocks, -c) 0
data seg size (kbytes, -d) unlimited
scheduling priority (-e) 0
file size (blocks, -f) unlimited
pending signals (-i) 6184327
max locked memory (kbytes, -l) 64
max memory size (kbytes, -m) unlimited
open files (-n) 1024
pipe size (512 bytes, -p) 8
POSIX message queues (bytes, -q) 819200
real-time priority (-r) 0
stack size (kbytes, -s) 8192
cpu time (seconds, -t) unlimited
max user processes (-u) 4096
virtual memory (kbytes, -v) unlimited
file locks (-x) unlimited
#注意：ulimit -n 4096是临时的改动，只对当前会话有效。
对比环节：for 循环和 parallel 运行时间对比。
我设置的 -j 40，parallel 的运行时间不到 for 循环批量的一半，68个基因，for循环运行了48min，parallel 只运行了23min，效率翻倍哇！
