系统发育分析教程

1、系统发育分析教程大致流程：1. 从18个mtDNA基因组中提取rRNA基因12S、16S和蛋白质基因ND1、ND2、CytB2. 分别进行序列比对，并进行比对精制3. 将精制比对结果串联成一个独立的分析文件，记录基因位置4. NJ分析(MEGA)5. MP分析(PAUP)6. ML分析（RAXML）7. 贝叶斯分析（MRBAYES）1.安装DNASTAR软件（又名Lasergene），软件内包含很多组件。2.例子中有18个转录组的数据，ctrl+A,点住第一个文件拖到DNASTAR的MegAlign里。确保MegAlign左侧的序列名称完全按照英文字母顺序来排。3.双击第一条序列，在出来的选框

2、中选取12S序列，点击NEXT。不断重复，直至将所有物种的12S序列挑出来。4.然后ctrl+A全选，点击OPTION下面的Genetic Codes,选择编码方式，根据基因来选，这里选择Vertebrate Mito。点击Align下面的By Clustal w Method等待程序对齐完成。这时的序列应该已经对齐了。5.将结果存为12S.MSF,MSF格式可以同时保存多个序列文件。6.重复2-5步，分别挑出16S、ND1、ND2、CytB，存为相应的名称。7.安装GeneStudioPro软件8. 打开GeneStudioPro的SeqVerter软件。点击Import sequences

3、导入序列，保留gaps全选序列，点击右侧Merge为一个Fasta序列。点击Clear清空，如此将所有序列处理完，将文件的后缀改为fas9.将改好名的文件复制入GBlocks的目录底下。10.打开GBlock.exe，输入o,回车输入上一步的文件名，回车输入t,回车，直到第一项t项为所选的序列类型输入g，回车，这时出现了两个文件重命名文件将-gb移动到.fas之前重复此步，将所有序列处理完，注意所选序列类型要正确。检查所有序列是否已切整齐，且为3的倍数。新建一个txt，命名为5genes打开txt，输入：序列类型，序列名称=起始位-终止位，基因按照特定顺序排列打开第一个序列，记录终止位置选择A

4、ppend alignment,按之前的顺序将序列全部导入，并记录下每个基因分布，即起始与终止位置，输入txt中将串联好的序列存为5genes.fas用mega打开序列选择分析然后选择核酸序列选择遗传密码选择distances/compute overall mean选择替代模型为nucleotide/jukes-cantor点击compute遗传距离为0.396，在0<x<1之间，适合建NJ树选择建NJ树选择对所有位点进行计算，假如蛋白编码的基因第三位替代过饱和，就选择1和2替代模型选择maximum composite likelihood这个是默认设置也可以改为另一个，填入之

5、前预测的模型的gamma参数再改bootstrap运行树已建好用seqverter将序列转为nex格式打开paup参数设置外类群设置搜索次数设置bootstrap次数设置brlens次数保存并退出打开paup，载入转换好的序列保存操作命令在操作行逐条输入并运行命令Outgroup 外类群Bootstrap nreps=1000 keepallContreeDescribetreesSavetrees from=1 to=1000或者直接打开刚刚做好的参数直接运行，等程序运行完产生了六个文件把树拖进treeview查看查看分数，一致性指数CI完全一致时为1，如果存在趋同进化或平行进化，则接近0.

6、保留指数RI与CI类似。将文件另存为phy4格式准备好phy格式的序列文件和txt格式的注释，复制到RAXML文件夹底下。打开RAXML下的AutoRun.txt文件，修改好相应参数，将后缀名改为bat。参数注解如下：-f 功能，选择了a,是最好用的，另外可以选择d，是最快的。 -m 模型类型，选择了GTRGAMMAI-s 序列所在文件名-n 后缀，自己设，这里设为5genes-q 基因分布所在文件名-# 分析10次Pause 停止运行Autorun.bat,开始跑数据。假如要分别对密码子的每一位进行独立的分析，如对编码蛋白的基因进行分析，则对基因分布文件进行修改，如改为DNA,ND_1=17

7、53-26713DNA,ND_2=1754-26713DNA,ND_3=1755-26713当程序跑完后，看最好的是哪次的结果。这次最好的结果是第一次run的，可以把结果拖进treeview里查看可以对Autorun进行修改而进行bootstrap，如改为：RAxML-7.0.3-WIN.exe -f a -x 12345 -p 12345 -m GTRGAMMAI -s 5genes.phy -n 5genes_boot50 -q 5genes.txt -# 50Pause程序产生了四个文件标尺0.1较为合适下面开始贝叶斯分析，用bioedit打开5genes.fas点击EXPORT-seq

8、uencealighment-nex/paup编辑此文件，查看missing=-;假如为missing=M gap=-;，则改为missing=-;将注解写在end;后begin mrbayes;charset 12S = 1 - 675;charset 16S = 676 - 1752;charset ND1 = 1753 - 2671;charset ND2 = 2672 - 3646;charset CYTB = 3647 - 4781;partition 5P = 5: 12S, 16S,ND1,ND2,CYTB;end;begin mrbayes; set partition=5P;

9、 Prset applyto=(all) statefreqpr=dirichlet(1,1,1,1) ratepr=variable; lset applyto=(all) nst=6 rates=invgamma; unlink shape=(all) pinvar=(all) statefreq=(all) revmat=(all);end;begin mrbayes;mcmc ngen=10000000 nruns=2 temp=0.2 samplefreq=1000 printfreq=1000;end;下面为注解Charset 基因片段Prset 先验，认为每种一样Lset likelyhood怎么设，all是所有partition分开，nst=6是GTR模型Unlike 所有模型非关联化Mcmc 加热，ngen为走多少步，nchains=4为默认值，可省略。在mcmc后加p可调出之前的数据，nrun为同时进行几个线程，最多为8个。Samplefreq为多少步确认一次Printfreq为多少步显示一次中括号为冷链，小括号为热链。 运行mrbayes，输入execute 5genes.nex当deviation为0，则两个结果一样，否则则不一样等待程序跑完，输入sumt burnin=50 这个值一般是所跑的四分之一。最后生成一致树，用treev