使用IQTree和INDELible生成层次树模拟数据集

# 1. 使用MAFFT等多序列比对软件将数据比对好（IQTree的输入数据是对齐好的序列）
# 生成碱基替换模型的数据集
$ mafft example_all.fas > example_all_aligned.fas
# 生成进化树的数据集
$ mafft example_100.fas > example_100_aligned.fas

# 2. 运行IQTree
# -s 输入数据， -B bootstrap次数， -m 测试模型， -T 并行数， --prefix 输出结果的前缀
# 不同版本的参数略有区别，具体使用方法查看./iqtree -h
# 生成碱基替换模型
$ ./iqtree -s PATH/example_all_aligned.fas -B 1000 -m TEST -T 20 --prefix example_model
# 生成进化树
$ ./iqtree -s PATH/example_100_aligned.fas -B 1000 -m TEST -T 20 --prefix example_tree

# 3. 查看example_model.iqtree文件确定碱基替换模型
# example_model.iqtree文件中的重要信息（模型和参数，这些需要输入到INDELible的control.txt中）：
# Model list: 确定碱基替换模型
# General matrix: 确定替换模型的参数（a b c d e f）
# Base frequencies: A: xxx  C: xxx  G: xxx  T: xxx
# Proportion of invariable sites: xxx
# Gamma shape alpha: xxx

# 4. 将example_tree.treefile文件中的冒号和冒号后的数字删除得到进化树

////////////////////////////////////////////////////////////////////////////////////
//                                                                                 //
//  INDELible V1.03 control file - basicnucleotide.txt                             //
//                                                                                 //
//      A basic introduction to the structure of the INDELible control file.       //
//                                                                                 //
/////////////////////////////////////////////////////////////////////////////////////

// It is useful to know that anything on a line after two forward slashes is ignored.

/*
   Another useful thing to know is that anything after a forward slash and star
   is ignored until INDELible sees a star followed by a forward slash later on.
*/

[TYPE] NUCLEOTIDE 1 //  EVERY control file must begin with a [TYPE] command.
            //  The number after "NUCLEOTIDE" can be 1 or 2 and chooses the
            //  algorithm that INDELible uses (see manuscript). Both give
            //  identical results but in some cases one is quicker.
            //  Other blocks and commands following this statement
            //  can come in any order you like.

[SETTINGS]
  [ancestralprint]           FALSE     // NEW, SAME or FALSE
  [output]                   FASTA    // FASTA, NEXUS, PHYLIP or PHYLIPT
  [phylipextension]          phy       // any alpha-numeric string
  [nexusextension]           nex       // any alpha-numeric string
  [fastaextension]           fas       // any alpha-numeric string
  [printrates]               FALSE     // FALSE or TRUE
  [insertaslowercase]        FALSE      // FALSE or TRUE
  [markdeletedinsertions]    FALSE     // FALSE or TRUE
  [printcodonsasaminoacids]  FALSE     // FALSE or TRUE
  [fileperrep]               TRUE     // FALSE or TRUE

[MODEL]  TVM_mt
  [submodel]     TVM 0.01 0.01 0.04 0.04        //  TVM: b=0.01, c=0.01, d=0.04, e=0.04,a=f=1
  [statefreq]    0.25 0.31 0.31 0.13                //  pi_T=0.25, pi_C=0.31, pi_A=0.31, pi_G=0.13
  [rates]        0.84 1.03 4                      //  pinv=0.84, alpha=1.03, ngamcat=4
  [indelmodel]   LAV  5 50                      //  Lavalette indel length distribution (a=5, M=50)
  [insertrate]   0.01                           //  insertion rate = 0.01 relative to substitution rate of 1
  [deleterate]   0.1                            //  deletion rate = 0.1 relative to substitution rate of 1

[TREE] mt100 (sq001,(((((((sq002,(sq044,sq053)),((sq024,sq077),sq071)),(sq052,sq065)),((sq014,(sq040,sq066)),sq076)),((((((((((((sq003,sq034),sq028),((((sq008,(sq018,sq029)),(sq037,sq056)),(((sq022,sq033),sq050),sq080)),(sq055,sq099))),((((((sq016,sq081),sq085),sq084),(((sq026,sq059),sq035),(sq047,sq069))),sq051),sq060)),sq091),(((sq010,sq041),sq082),sq094)),((sq012,sq062),(sq025,sq075))),((sq007,(sq046,sq087)),(((((((sq011,(sq021,sq089)),(sq067,sq096)),((sq039,sq079),sq093)),sq063),sq057),(sq015,sq054)),(sq032,(sq073,sq095))))),(((sq045,sq070),sq078),sq072)),((sq013,(((sq030,sq048),sq036),sq088)),(sq064,sq090))),(sq042,sq100)),sq058)),(((((((sq004,sq009),sq043),sq023),sq083),sq038),sq017),(((sq005,sq092),(((sq006,sq098),(sq031,sq097)),sq074)),sq049))),((sq019,sq086),((sq020,sq061),sq027))),sq068);
  [treelength] 0.22
  [branchlengths] NON-ULTRAMETRIC


[PARTITIONS] partitionname             //  [PARTITIONS] blocks say which models go with
  [mt100 TVM_mt 16000]            //  which trees and define the length of the
                                       //  sequence generated at the root (1000 here).

[EVOLVE] partitionname 10 tl0022  //  This will generate 100 replicate datasets
                                       //  from the [PARTITIONS] block named above.

// The true alignment will be output in a file named outputname_TRUE.phy
// The unaligned sequences will be output in a file named outputname.fas
// To learn how to implement more complicated simulations (or different
// models) please consult the manual or the other example control files.s

# INDELible使用方法
$ ./indelible PATH/control.txt