被子植物·分子系统发育分析·实用命令

（1） 序列比对：MAFFT

  主要设置分析模式：精度优先（3种）、速度优先（7种）、组对组对齐（Group-to-group alignments，1种）。其余默认。
     组对组对齐：两个分别进行过序列比对的数据集利用该模式进行合并对齐。
     依数据判断选择何种模式。
  Manual: https://mafft.cbrc.jp/alignment/software/manual/manual.html#lbAG

（2） 序列修剪：Trimal

  例：trimal -in input.fasta -out output.fasta -automated1（足够使用）；-nogaps参数能将所有含有gap的对齐位点剔除。

（3） 模型预测：PartitionFinder2

  例：Python PartitionFinder.py <Input.cfg> -p 6 #-p限制用于计算的处理器数目，Python可省，其余默认
    运行脚本（Input.cfg）示例：

## ALIGNMENT FILE ##
alignment = input.phy;   # input.phy需与本脚本文件放在同一文件夹
## BRANCHLENGTHS: linked | unlinked ##
branchlengths = unlinked; # 取决于后续系统分析软件
## MODELS OF EVOLUTION: all | allx | mrbayes | beast | gamma | gammai | <list> ##
models = mrbayes; # 取决于后续系统分析软件与希望达到的预测精度
# MODEL SELECCTION: AIC | AICc | BIC #
model_selection = aicc; # Manual建议选aicc
## DATA BLOCKS: **see manual for how to define** ##
[data_blocks]
ITS = 1-201; # 关键，可针对多基因联合的数据集
## SCHEMES, search: all | user | greedy | rcluster | rclusterf | kmeans ##
[schemes]
search = greedy; # 根据数据集决定使用何种算法

（4） 构建基因树、蛋白树

  ① 极大似然法：iqtree
  例：iqtree -s input_file -m MFP -bb 1000 -nt 16 -pre its.ml1 -quiet # 模型预测、树构建一起干，1000次超快自举检验，16个线程，输出文件前缀设置its.ml1，运行过程不打印
  对于包含外类群的矩阵，iqtree无需指定外类群。运行结果得到无根树，只需将所得结果在Figtree之类的软件中打开，使用“reroot”将已知的外类群拉出即可。
  ② 贝叶斯推断：MrBayes
  例：mpirun -np 4 mb input.nexus #四线程版本并行MrBayes，mb只能绝对路径或相对路径
    运行脚本（input.nexus）示例：

#NEXUS 
BEGIN TAXA;
       DIMENSIONS NTAX=n;    #NTAX数目等于n
       TAXLABELS
          [taxa_list]    #与character.block的物种顺序相同
;
END;

BEGIN TREES;
    tree mytree = [tree.newick];   #如果有根树，在[tree.newick]前加[&R]代号示意。
END;
# tree_block必须在taxa_block之后。且不能在data_block之前，data_block会重新定义taxa与data，
# 因此当要使用引导树时，只能用character_block装矩阵，不能用data_block。

BEGIN CHARACTERS;
    Dimensions nchar=n;    #nchar: character numble, towards DNA, i.e. base number
    Format datatype=dna missing=N interleave=yes gap=-;
    Matrix
[interleave matrix]
;
End

BEGIN MRBAYES;
Outgroup name1;
Outgroup name2;
    [partition info] or such as below;
Lset nst=2 rates=gamma; #Lset设置核苷酸替代模型，nst设置核苷酸转换与颠换的模式，rates设置核苷酸替换速率模式。
Prset statefreqpr=dirichlet(1,1,1,1); #Prset设置核苷酸替代的先验概率，尽管核苷酸替代模型的估计往往表明样本核苷酸拥有不同的替代速率，然而一般的分析保守地假定核苷酸频率是1:1:1:1（如同说明书中表述）。

mcmcp nperts=3 ngen=2000000 printfreq=1000 samplefreq=100 nchains=4 savebrlens=yes filename=output_name; 
# 马尔可夫链蒙特卡罗链运行命令，ngen代数，printfreq屏幕打印频率，
# samplefreq样本取样频率，nchains链数目，savebrlens说明书没记载但许多网络教程都有该项。
# nperts 略微改变引导树，避免不同运行难以检测收敛。
startvals tau=mytree;    #设置引导树，必须在mcmcp之后
mcmc;
sumt relburnin=yes burninfrac=0.25 filename= output_name Contype=Allcompat;
#输出树结果与参数值，relburnin=yes burninfrac=0.25与mcmcp设置的一致，
#此处即默认值，说明书解释了该处不能默认不写的原因。
sump burninfrac=0.25 filename= output_name; #输出样本参数值（样本即碱基）
END;

## burnin=5000，burnin/(ngen/samplefreq)=burninfrac：表示样本前burninfrac的数据被去掉，
## 因为这部分数据统计是不可靠的。默认为前25%的样本被去掉。该比例可以控制在10~25%。
## mcmc命令burnin缺失时，默认去掉25%的样品统计结果信息。
## 默认设置是精良的，一般不用变。需要改变的即前面脚本中的数值。

  运行MrBayes时，有时候会中途断掉，或者跑完迭代数（generation number）后出现了错误，为了避免重新从头开始跑，可以使用Mrbaye的checkpoint文件（.ckp，默认中断运行后输出），只需要在mcmcp命令后增加append=yes，到原目录重新启动程序即可接续分析。
  有教程提到MrBayes分析的线程要设置≤Nruns（独立搜索次数，默认2）×nchains（马尔科夫链数，默认4），因此可以设置-np 4或8。ps. 多线程经测试相当不稳定，分析经常中断，且四线程的速度远没有单线程的四倍，对于大数据集不推荐多线程运行。
  未设置的参数默认，具体参考manual。
  如果用来建树的序列太长，建议将序列改成interleave形式（交错排列式）的连续.nex文件。此时只需在.nex文件BEGIN CHARACTERS; FORMAT增加interleave=yes即可。
  将non-interleave nexus转换成interleave nenux可通过paup实现，在paup中进入待转换的nexus文件目录，输入命令export file=waiting_conversion.nexus format=nexus interleaved=yes;运行，建议使用原文件的备份件，因为该操作将覆盖原始的文件（我不太清楚paup如何进行重定向）。

（5） ……（待续）

参考文献：

[1] Partitionfinder-Manual
[2] MrBayes-Manual
[3] MrBayes建树的关键命令脚本与注意事项
[4] 基于PartitionFinder做基因分区、模型选择，并用MrBayes建树案例与参数解析（详细流程）
[5] 多线程MrBayes安装教程及使用（Mpi_MrBayes）
[6] Mrbayes构建贝叶斯树复习
[7] LD_LIBRARY_PATH详解