06进化树流程-常规
Tips：
1.raxmlHPC-PTHREADS-SSE3 -o参数为root,与外层程序-r都是根
2.外群名对应fa文件里的名字
3.模型预测时,要输入jmodel能识别的格式(如.fas) 最好先用mega转换一下
〇.矫正
0.1 下载ref
## 复制表格中的物种 至 tmp_tre文件(以+号为分割)
编号目录]
awk -F "+" '{ print $2 }' tmp_tre | awk '$1=$1' > list_tre && cat list_tre && down- tre && cd ref_tre/ && cat fasta/* > allref.fa && cd ../
$ python3 /share/nas1/yuj/script/chloroplast/download_genome_from_ncbi.py -f gb -a list_tre # 备用下载
## 复制组装结果至ref_tre文件夹
1)叶绿体
for i in $( ls analysis/assembly) ;do echo $i;cp analysis/assembly/$i/finish/*.fsa ref_tre/$i.fsa;done
for i in $( ls analysis/assembly) ;do echo $i;cp analysis/annotation/$i/gene_anno/$i.gbk ref_tre/gbk/;done
2)线粒体
for i in $( ls analysis/assembly) ;do echo $i;cp analysis/assembly/$i/finish/*.fsa ref_tre/$i.fsa;done
for i in $( ls analysis/assembly) ;do echo $i;cp analysis/annotation/$i/$i.gbk ref_tre/gbk/;done
0.2 全基因组矫正
0.2.1 线粒体
(1)全部修改起点
cd ref_tre/fasta
for i in *.fasta*;do nuc ../M23_FULLMT.fsa $i -l 1 && mt_move_pos.py -fa $i -s `awk 'NR==4 {print $3}' out.delta` ;done
cd ../ && mkdir fasta00 && mv fasta/*.fasta fasta00 && cd fasta && rename fsa2 fasta *.fsa2 && cd ../ && cat fasta/*.fasta *.fsa > allinone.fa
nuc *.fsa allinone.fa && mum
0.2.2 叶绿体
0.看一眼基因数量
cd ref_tre && cp_get_cds.py -i gbk/ -o out/
# 提取完整序列，考虑到一些有重复的项目,默认去重,-d不去重
1.查看长度 去掉长度明显异常的(如4204项目的高级分析物种掺入了线粒体)
ref_tre]$ cd fasta && for i in * ;do ir.py -i $i -l;done && cd ../../
2.看ir程序结果 判断各物种序列情况(是否起点有错)
GP-XXX]$ cd ref_tre/fasta/ && for i in *.fasta;do ir $i ;done > ir.log
(1)(批量修改)起点不对
fasta]
python3 /share/nas1/yuj/script/chloroplast/phytree/cp_batch_adjust_genome_start.py -i1 ./
# 等待几秒钟,生成ref_tre/complete2 文件夹
(2)执行下面命令再检查一遍 可以省略该步骤
fasta]$ cd ../complete2 && for i in *.fasta;do ir $i ;done > ir.log && cd ../
3.看nuc程序结果 判断各物种序列情况(ssc反向/整体反向)
ref_tre]$ cat complete2/*.fasta *.fsa > allinone.fa
ref_tre]$ nuc *.fsa allinone.fa && mum
# ssc反向调整(叶绿体)
$ cp_ssc.pl -i -o
一.Full Genome
1.1 最大似然法 FULLML
1.1.1 叶绿体/线粒体
ref_tre]$
perl /share/nas6/pub/pipline/genome-assembly-seq/chloroplast-genome-seq/v1.2/phytree/raxml_phytree.pl -i allinone.fa -p fullml -o full_ml_phytree &
# perl /share/nas6/pub/pipline/genome-assembly-seq/chloroplast-genome-seq/v1.2/phytree/raxml_phytree.pl -i allinone.fa -p 物种名(可以随便起) -r 外群名(多个外群用,号隔开) -o full_ml_phytree &
1.1.2 分步拆解
0.创建文件夹
ref_tre]$
mkdir full_ml_phytree && cd full_ml_phytree/
1.mafft比对
full_ml_phytree]$
mafft --auto --quiet --thread 30 ../allinone.fa > fullml.aln
# 去掉换行 perl /share/nas6/xul/program/mt2/phytree/gene_tree/src/fasta2line.pl -i fullml.aln -o fullml.aln.fasta
2.trimal修剪
full_ml_phytree]$
/share/nas6/xul/soft/trimal/trimal-1.4.1/source/trimal -in fullml.aln -out fullml.trim.aln -automated1
# perl /share/nas6/xul/program/mt2/phytree/gene_tree/src/fasta2line.pl -i allinone.aln -o all.fasta.line # 转换成一行序列
3.raxml建树
full_ml_phytree]$
mkdir raxml && touch raxml/RAxML.log
raxmlHPC-PTHREADS-SSE3 -f a -s fullml.trim.aln -T 8 -# 1000 -m GTRGAMMA -x 12345 -p 12345 -n genome -w `realpath raxml` 1> raxml/RAxML.log
# raxmlHPC-PTHREADS-SSE3 -f a -s fullml.trim.aln -T 8 -# 1000 -m GTRGAMMA -o [外群名] -x 12345 -p 12345 -n genome -w raxml的绝对路径 1> raxml/RAxML.log
4.复制改名进一步处理
full_ml_phytree]$
cp raxml/RAxML_bipartitionsBranchLabels.genome fullml.genome.nwk
5.第3步、第4步合并
mkdir raxml && touch raxml/RAxML.log && raxmlHPC-PTHREADS-SSE3 -f a -s fullml.trim.aln -T 8 -# 1000 -m GTRGAMMA -x 12345 -p 12345 -n genome -w `realpath raxml` 1> raxml/RAxML.log && cp raxml/RAxML_bipartitionsBranchLabels.genome fullml.genome.nwk &
## 模型预测
java -jar /share/nas6/zhouxy/biosoft/jmodeltest/jmodeltest-2.1.10/jModelTest.jar -i -f -g 4 -BIC -AIC -AICc -DT -tr 15 -o jmodeltest -d fullml.trim.aln &
1.2 贝叶斯法 FULLBI
## 前面ml有结果，直接将“aln.fa”和“模型预测结果”输入以下程序
full_ml_phytree]$
python3 /share/nas1/yuj/script/chloroplast/phytree/or_convert_to_nex.py -i fullml.trim.aln -m jmodeltest -o fullbi.nex
cd ../ && mkdir full_bi && mv full_ml_phytree/fullbi.nex full_bi/ && cd full_bi
mb -i fullbi.nex & # s1节点运行报错
# 完整步骤
1.采用全基因组做进化树分析，将环形序列设置相同起点，物种间序列用MAFFT软件（v7.427，--auto模式）进行多序列比对（.fasta）
2.将比对好的数据（.fasta）用jModelTest计算最优核苷酸替代模型，在贝叶斯文件夹里，详见 https://www.jianshu.com/p/e5cfad89a1a4
$ java -jar /share/nas6/zhouxy/biosoft/jmodeltest/jmodeltest-2.1.10/jModelTest.jar -i -f -g 4 -BIC -AIC -AICc -DT -tr 15 -o jmodeltest -d Single_gene_seq.trim.aln
3.第1步比对好的数据（.fasta）导入Phylosuite,并根据第2步的结果准备贝叶斯法建树的.nex文件,修改详见https://www.jianshu.com/p/8b10ef5c26e1
4.将.nex 文件导入MrBayes软件进行贝叶斯建树(参数已包含在.nex文件里)
$ mb -i fullbi.nex &
1.3 邻接法 FULLNJ
/share/nas1/yuj/software/rapidNJ_2.3.3/bin/rapidnj -i fa fullml.trim.aln -b 1000 > fullnj.nwk
二.Shared CDS
2.0 Tips
# 叶绿体
# 取目标文件CDS
$ perl /share/nas6/xul/program/chloroplast/phytree/gene_tree/find_total_gene.pl -i gbk/ -p KX452726物种名/文件名 -o gene # gbk放目标文件及随便一个参考
# gene/feature/ref.gene.seq为参考的所有cds
# gene/feature/$prefix.gene.seq为目标文件的所有cds
$ cp_annotation_one_gene_by_ref_gbk2.pl -i1 fasta/NC_026867.1.fasta -i2 ref/gbk/Tkuehnei.gbk -g 基因名 # 根据参考gbk查找基因
$ cp_annotation_one_gene_by_ref_gbk2.pl -i1 Thainanus_FULLMT.fsa -i2 ref_ann/gbk/NC_050664.1.gbk -g COX1 # 根据参考gbk查找基因
$ perl /share/nas6/xul/program/mt2/phytree/gene_tree/src/fasta2line.pl -i mafft/gene1.COX1.fasta.aln -o mafft/gene1.COX1.fasta # 比对去掉分行
# 线粒体
$ perl /share/nas6/pub/pipline/genome-assembly-seq/mitochondrial-genome-seq/v2.0/advance/phytree/raxml_phytree.pl # 建树
PS:
D:\ProgramData\Anaconda3\lib\site-packages\Bio\GenBank\__init__.py 352 358 1236
seq_type = 'dna circular' # 20220519 自己添加的 1229
改了biopython库,以后要确保这些gbk都是环形的才行,不然就要改回去
2.1 最大似然法 CDSML
2.1.0 真菌
1.ref_tre]
python3 /share/nas1/yuj/script/fungi/phytree/fungi_from_gbk_get_cds_V1.0.py -i gbk/ -o out_all -d && mkdir cds_all && cp out_all/cds/* cds_all # - -d 保留重复
python3 /share/nas1/yuj/script/fungi/phytree/fungi_extract2mafft_V1.5.py -i cds_all/ -o1 gene/extract -o2 gene/mafft -c2
2.1简单项目----删除非公有cds
ref_tre]$
cd gene/mafft && for i in *.fasta;do echo $i `grep ">" $i | wc -l` ;done | sort -V
2.2有些问题----先从上一个文件夹开始修改
2.2.1删除非公有cds
cd gene/extract
for i in *.fasta;do echo $i `grep ">" $i | wc -l` ;done | sort -k 2 -r -V
rm 非公有cds
2.2.2找出未注释的cds或去掉多余的一个
!!!注意!!!>id相同的话,只会显示前一条序列的结果
3. mafft]$
perl /share/nas6/xul/program/mt2/phytree/connection_head_tail.pl -i *.fasta -o ../
# 后缀有2的脚本是另一种排序 perl /share/nas6/xul/program/mt2/phytree/connection_head_tail2.pl -i *.fasta -o ../
# 若都不好用,使用叶绿体的脚本 perl /share/nas6/xul/program/chloroplast/phytree/connection_head_tail.pl -i *.fasta -o ../
4. mafft]$
cd ../ && /share/nas6/xul/soft/trimal/trimal-1.4.1/source/trimal -in gene_seq.aln -out gene_seq.trim.aln -automated1
# 可能需要把gene_seq.aln改为.fasta文件
5. gene]$
nohup perl /share/nas6/pub/pipline/genome-assembly-seq/chloroplast-genome-seq/v1.2/phytree/raxml_phytree.pl -i gene_seq.trim.aln -aln -o cds_ml_phytree &
# 建议采用叶绿体程序,为 后续 cds贝叶斯 做准备
# 线粒体自己的程序,不会进行模型预测
# nohup perl /share/nas6/pub/pipline/genome-assembly-seq/mitochondrial-genome-seq/v2.0/advance/phytree/raxml_phytree.pl -i gene_seq.trim.aln -aln -o cds_ml_phytree &
2.1.1 线粒体
程序路径：
线粒体
python3 /share/nas1/yuj/script/mitochondrion/phytree/mt_from_gbk_get_cds_V1.0.py -h
python3 /share/nas1/yuj/script/mitochondrion/phytree/mt_extract2mafft_V1.5.py -h
叶绿体
python3 /share/nas1/yuj/script/chloroplast/phytree/cp_from_gbk_get_cds_V2.5.py -h
python3 /share/nas1/yuj/script/chloroplast/phytree/cp_extract_and_mafft_V1.0.py -h
1. ref_tre]
mt_get_cds.py -i gbk/ -o out_all -d && mkdir cds_all && cp out_all/cds/* cds_all # -d 保留重复
# python3 /share/nas1/yuj/script/mitochondrion/phytree/mt_from_gbk_get_cds_V2.0.py -i gbk/ -o out_all && mkdir cds_all && cp out_all/cds/* cds_all
for i in cds_all/*.fasta;do if [ ! -s $i ]; then echo $i && rm $i; fi; done # 看看有没有空文件,有就删除
python3 /share/nas1/yuj/script/mitochondrion/phytree/mt_extract2mafft_V1.5.py -i cds_all/ -o1 gene/extract -o2 gene/mafft -c2
2. ref_tre]$
cd gene/mafft
3. mafft]$
perl /share/nas6/xul/program/mt2/phytree/connection_head_tail.pl -i *.fasta -o ../
# 后缀有2的脚本是另一种排序 perl /share/nas6/xul/program/mt2/phytree/connection_head_tail2.pl -i *.fasta -o ../
# 若都不好用,使用叶绿体的脚本 perl /share/nas6/xul/program/chloroplast/phytree/connection_head_tail.pl -i *.fasta -o ../
4. mafft]$
cd ../ && /share/nas6/xul/soft/trimal/trimal-1.4.1/source/trimal -in gene_seq.aln -out gene_seq.trim.aln -automated1
5. gene]$
nohup perl /share/nas6/pub/pipline/genome-assembly-seq/chloroplast-genome-seq/v1.2/phytree/raxml_phytree.pl -i gene_seq.trim.aln -aln -o cds_ml_phytree &
# 输入的是比对好的文件 -aln
# cannot be read as an alignment错误是因为序列🆔包含了“.”，需要去掉“.”
# 建议采用叶绿体程序,为 后续 cds贝叶斯 做准备
## 线粒体自己的程序,不会进行模型预测
# nohup perl /share/nas6/pub/pipline/genome-assembly-seq/mitochondrial-genome-seq/v2.0/advance/phytree/raxml_phytree.pl -i gene_seq.trim.aln -aln -o cds_ml_phytree &
2.1.2 叶绿体
1. ref_tre]$
cp_get_cds.py -i gbk/ -o out_all -d && mkdir cds_all && cp out_all/cds/* cds_all
# -d 保留重复
for i in cds_all/*.fasta;do if [ ! -s $i ]; then echo $i && rm $i; fi; done # 看看有没有空文件,有就删除
cd cds_all
sed -i 's/pafII/ycf4/g' *
sed -i 's/pafI/ycf3/g' *
sed -i 's/pbf1/psbN/g' *
sed -i 's/clpP1/clpP/g' *
cd ..
python3 /share/nas1/yuj/script/chloroplast/phytree/cp_extract_and_mafft_V2.0.py -i cds_all/ -o1 gene/extract -o2 gene/mafft -c2
2. ref_tre]$
cd gene/mafft
/share/nas6/xul/program/chloroplast/phytree/connection_head_tail.pl -i *.fasta -o ../ # 可能会有bug,目前只出现过一次,可以换用 perl /share/nas6/xul/program/mt2/phytree/connection_head_tail.pl -i *.fasta -o ../
cd ../ && /share/nas6/xul/soft/trimal/trimal-1.4.1/source/trimal -in Single_gene_seq.aln -out gene_seq.trim.aln -automated1
# 可能需要把gene_seq.aln改为.fasta文件
3. gene]$
$ nohup perl /share/nas6/pub/pipline/genome-assembly-seq/chloroplast-genome-seq/v1.2/phytree/raxml_phytree.pl -i gene_seq.trim.aln -aln -o cds_ml_phytree &
# 建议采用叶绿体程序,为 后续 cds贝叶斯 做准备
nohup perl /share/nas6/pub/pipline/genome-assembly-seq/chloroplast-genome-seq/v1.2/phytree/raxml_phytree.pl -i gene_seq.trim.aln -aln -r Glycine_max_cultivar_PI_548384_MW357267,Glycine_soja_KY241814 -o cds_ml_phytree2 &
2.1.3 叶绿体(xul版)
1.
$ perl /share/nas6/xul/program/chloroplast/phytree/gene_tree/find_total_gene.pl -i gbk/(包含所有gbk) -p 物种名/文件名(必须是测序物种名) -o gene
# 以下为分步子程序,一般不用
$ perl /share/nas1/xul/program/chloroplast/cds_analysis/src/extract.dir.cds.pl -i gbk/ -p MN649876 -o gene/feature # 提取cds,文件ref.gene.seq
$ perl /share/nas1/xul/program/chloroplast/cds_analysis/src/blast2match.no.ir.pl -i gene/feature/MN649876.gene.seq -r gene/feature/ref.gene.seq -o gene/blast
$ perl /share/nas1/xul/program/chloroplast/cds_analysis/src/mafft.dir.fasta.pl -i gene/blast/fasta -o gene/mafft
$ perl /share/nas2/xul/program/chloroplast/cds_analysis/src/mafft.dir.fasta.pl -i gene/blast/fasta -o gene/mafft
2.cd gene/mafft
3.nu=`ls ../../gbk/*.gbk |wc |awk '{print $1}'` && perl /share/nas6/xul/program/chloroplast/phytree/connection_head_tail.pl -i ` ls |grep -v trn |grep -v rrn |xargs grep -c ">" |grep ":$nu" |cut -d: -f 1 |tr "\n" " " ` -o ../../gene
# 若上述不好用,简化一下输入
$ perl /share/nas6/xul/program/chloroplast/phytree/connection_head_tail.pl -i *.fasta -o ../
4.cd ../ && /share/nas6/xul/soft/trimal/trimal-1.4.1/source/trimal -in Single_gene_seq.aln -out Single_gene_seq.trim.aln -automated1
5.nohup perl /share/nas6/pub/pipline/genome-assembly-seq/chloroplast-genome-seq/v1.2/phytree/raxml_phytree.pl -i Single_gene_seq.trim.aln -aln -o cds_ml_phytree & # 这个可以运行,但是可能没比对好
2.1.4 线粒体(xul版)
1.
$ perl /share/nas6/xul/program/mt2/phytree/gene_tree/find_total_gene.pl -i gbk/ -p 物种名/文件名(必须是测序物种名) -o gene
# 以下为分步子程序,一般不用,用的话需要判断
$ perl /share/nas6/xul/program/mt2/phytree/gene_tree/src/extract.dir.cds.pl -i $indir -p $prefix -o $outdir/feature
# 判断gene/feature/ref.gene.seq数量,基因名字保持一致(以大写为准),填上没有显示的基因名
$ perl /share/nas6/xul/program/mt2/phytree/gene_tree/src/blast2match.no.ir.pl -i gene/feature/$prefix.gene.seq -r gene/feature/ref.gene.seq -o gene/blast
# 判断gene/blast/fasta/数量及方向(方向应该没错,数量前面判断)
$ perl /share/nas6/xul/program/mt2/phytree/gene_tree/src/mafft.dir.fasta.pl -i $outdir/blast/fasta -o $outdir/mafft
$ mafft --auto --quiet --thread 30 gene/mafft/gene5.atp8.fasta1 > gene/mafft/gene5.atp8.fasta # 手动比较的要自己去掉换行,单个单个运行
2.2 贝叶斯法Bayes CDSBI
2.2.1 以线粒体为例
# 之前有cds ml结果
5.gene]$
$ python3 /share/nas1/yuj/script/chloroplast/phytree/or_convert_to_nex.py -i gene_seq.trim.aln -m cds_ml_phytree/jmodeltest -o cdsbi.nex && mkdir cds_bi && mv cdsbi.nex cds_bi/ && cd cds_bi
6.将.nex 文件导入MrBayes软件进行贝叶斯建树(参数已包含在.nex文件里)
$ mb -i cdsbi.nex
# 完整步骤，前4步同cdsml
1.ref_tre]$
$ python3 /share/nas1/yuj/script/mitochondrion/phytree/mt_from_gbk_get_cds_V1.0.py -h
$ python3 /share/nas1/yuj/script/mitochondrion/phytree/mt_extract2mafft_V1.0.py -h
2.cd gene/mafft
3.perl /share/nas6/xul/program/mt2/phytree/connection_head_tail.pl -i *.fasta -o ../ # 后缀有2的脚本是另一种排序 perl /share/nas6/xul/program/mt2/phytree/connection_head_tail2.pl -i *.fasta -o ../
# 若都不好用,使用叶绿体的脚本 perl /share/nas6/xul/program/chloroplast/phytree/connection_head_tail.pl -i *.fasta -o ../
4.cd ../ && /share/nas6/xul/soft/trimal/trimal-1.4.1/source/trimal -in gene_seq.aln -out gene_seq.trim.aln -automated1
# 把gene_seq.aln改为.fasta文件
5.第4步比对好的数据（改名为.fasta）导入Phylosuite,并根据上述提到ML法的模型预测结果准备贝叶斯法建树的.nex文件(看BIC那一行)
$ java -jar /share/nas6/zhouxy/biosoft/jmodeltest/jmodeltest-2.1.10/jModelTest.jar -i -f -g 4 -BIC -AIC -AICc -DT -tr 15 -o jmodeltest -d Single_gene_seq.trim.aln # 模型预测
6.将.nex 文件导入MrBayes软件进行贝叶斯建树(参数已包含在.nex文件里)
$ mb -i cdsbi.nex
2.2.2 贝叶斯参数
Note
nst
JC, nst=1
F81, nst=1
K80, nst=2
HKY, nst=2
GTR, nst=6
rates 定义位点之间的替换率，有以下几种选择：
equal：位点的替换率无差异
gamma：位点的替换率呈 gamma 分布，对应+G
adgamma：位点的替换率自相关，边缘位点替换率呈 gamma 分布，相邻位点有相关的替换率
propinv：一定比例位点的替换率是恒定的，对应+I
invgamma：一定比例位点的替换率是恒定的，剩下位点的替换率呈 gamma 分布，对应+I+G。
Ngammacat：配合上面的参数，如果替换速率设置为Gamma、Invgamma、Adgamma，则需要设置此选项。一般为4或5,增加会更准确但是会变慢,4折中
shapepr = fixed(1.1630)，括号内的为gamma值；
----示例 GTR+I+G
lset nst=6 rates=invgamma Ngammacat=5;
prset statefreqpr = fixed(0.3928,0.0901,0.0754,0.4416) revmat = fixed(0.7601,8.4971,0.7601,1.0000,8.4971,1.0000) pinvar=fixed(0.3690) shapepr = fixed(1.1630);
2.3 邻接法 CDSNJ
前面步骤都一样
/share/nas1/yuj/software/rapidNJ_2.3.3/bin/rapidnj -i fa gene_seq.trim.aln -b 1000 > phytree.nwk
三.修改名字
程序路径：
python E:\OneDrive\jshy信息部\Script\chloroplast\phytree\phytree_trans_nwk_name_V2.0.py -h
python3 /share/nas1/yuj/script/chloroplast/phytree/phytree_trans_nwk_name_V3.0.py -h
Tips：
Bayes树处理的时候,对最开始的文件进行修改,不要先转格式,改完名再转格式
使用前先查看待修改的树文件,通常来说运行-f1
如出现f2中的两种情况,则使用-f2
都用不了再-f3,输入准备好的文件
0.xul版,先不用
ref_tre]$ perl /share/nas6/xul/project/xianliti2/GP-20200617-2017_xizangnongmuxueyuan_1sample_lingdangzi_xianliti2/analysis/phytree/ref/get_id_genbank.pl gbk/*.gbk > id.list
# perl /share/nas6/xul/project/xianliti2/GP-20200617-2017_xizangnongmuxueyuan_1sample_lingdangzi_xianliti2/analysis/phytree/ref/trans_name.pl phytree/sample.cds.nwk > cds.nwk 可以不用
1.提取ID
ref_tre]$
# 适用于-f3(暂时用着)
cd fasta/ && for i in *.f*;do awk '{if(NR==1) print $0}' $i;done | awk -F ">" '{ print $2 }' > ../ori.list && cd ../ && cd fasta/ && for i in *.f*;do awk '{if(NR==1) print $0}' $i;done | awk -F ">" '{ print $2 }' | awk -F ".1_" '{ print $2"""_"""$1""".1" }' > ../new.list && cd ../ && for i in `cat ori.list`;do echo "${i//.1_/_1_}"; done > orifullbi.list
cd complete2/ && for i in *.f*;do awk '{if(NR==1) print $0}' $i;done | awk -F ">" '{ print $2 }' > ../ori.list && cd ../ && cd complete2/ && for i in *.f*;do awk '{if(NR==1) print $0}' $i;done | awk -F ">" '{ print $2 }' | awk -F ".1_" '{ print $2"""_"""$1""".1" }' > ../new.list && cd ../ && for i in `cat ori.list`;do echo "${i//.1_/_1_}"; done > orifullbi.list
$ grep ">" allinone.fa | awk -F ">" '{print $2}' > ori.list && grep ">" allinone.fa | awk -F ">" '{print $2}' | awk -F "_" '{print $2"""_"""$3"""_"""$4"""_"""$5"""_"""$6"""_"""$1}' > new.list # 不太好用
# 适用于-f1(程序修改后不太好用)
# perl /share/nas6/xul/project/xianliti2/GP-20200617-2017_xizangnongmuxueyuan_1sample_lingdangzi_xianliti2/analysis/phytree/ref/get_id_genbank.pl gbk/*.gbk > id.list
2.改名
ref_tre]$
cd full_ml_phytree && python3 /share/nas1/yuj/script/chloroplast/phytree/phytree_trans_nwk_name_V3.0.py -f3 -id1 ../ori.list -id2 ../new.list -2 *.genome.nwk -3 fullml.nwk && cd ../
cd full_bi && python3 /share/nas1/yuj/script/chloroplast/phytree/phytree_trans_nwk_name_V3.0.py -f3 -id1 ../orifullbi.list -id2 ../new.list -2 *.con.tre -3 fullbi.tre && cd ../
cd gene/cds_ml_phytree && python3 /share/nas1/yuj/script/chloroplast/phytree/phytree_trans_nwk_name_V3.0.py -f3 -id1 ../../ori.list -id2 ../../new.list -2 *.genome.nwk -3 cdsml.nwk && cd - # 不需要
cd gene/cds_ml_phytree && cp sample.genome.nwk cdsml.nwk && cd -
cd gene/cds_bi && python3 /share/nas1/yuj/script/chloroplast/phytree/phytree_trans_nwk_name_V3.0.py -f3 -id1 ../../ori.list -id2 ../../new.list -2 *.con.tre -3 cdsbi.tre && cd -
3.
ref_tre]$
mkdir phytree && cp full_ml_phytree/fullml.nwk phytree && cp full_bi/fullbi.tre phytree && cp gene/cds_ml_phytree/cdsml.nwk phytree && cp gene/cds_bi/cdsbi.tre phytree
四.方法描述
2.6 系统发生分析
选取木犀科，特别是女贞属已经发表的叶绿体基因组序列进行系统发生关系重建。 数据集包括 6 个女贞属的物种和其他 21 个木犀科的物种。 叶绿体全基因组使用 MAFFT 软件进行比对，利用 Gblocks 0.91b 删除序列比对不可靠的区域。 利用最大似然法和贝叶斯的方法分别构建系统发生关系。 在贝叶斯信息准则下，使用 ModelFinder 寻找最优的核苷酸替代模型。 最大似然树使用 RAxML v.8.1.24 软件进行构建，模型使用 ModelFinder 筛选出的最优模型 GTR+G，分支的拓扑结构的支持率（BS，bootstrap）使用 500 次重复抽样进行计算。 贝叶斯树使用 MrBayes 3.0 进行计算，运行 2 个重复，每个重复包含 4 条马可夫链（1 条冷链和 3 条热链），总运行 20,000,000 代，抽样频率为 1,000，运行结束后使用 Tracer 1.6 查看，当达到收敛并且所有参数均取样足够时（ESS > 200），舍弃前 25% 的抽样（burn in）。 合并剩余的抽样，生成 50% 的多数一致树（50% majority-rule consensus tree），并计算各分支的后验概率（PP，posterior probability）。
采用全基因组做进化树分析，将环形序列设置相同起点，
使用MAFFT软件（v7.427，--auto模式）进行多序列比对，利用trimAl(v1.4.rev15) 删除序列比对不可靠的区域。
利用最大似然法和贝叶斯的方法分别构建系统发生关系。
在贝叶斯信息准则下，使用jModelTest v2.1.10寻找最优的核苷酸替代模型。
最大似然树使用 RAxML v.8.2.10 软件进行构建，模型使用jModelTest v2.1.10 筛选出的最优模型 GTR+I+G，
分支的拓扑结构的支持率（BS，bootstrap）使用 1000 次重复抽样进行计算。
贝叶斯树使用 MrBayes v3.2.7进行构建，马尔科夫链蒙特卡洛方法（Markov Chain Monte Carlo, MCMC），
运行 2 个重复，每个重复包含 4 条马可夫链（1 条冷链和 3 条热链），
总运行 20,000,000 代，抽样频率为 1,000，舍弃前 25% 的抽样（burn in），
最终得到多数一致树（majority-rule consensus tree）。
全基因组贝叶斯法
1.
采用全基因组做进化树分析，将环形序列设置相同起点，物种间序列用MAFFT v7.427软件（--auto模式）进行多序列比对，将比对好的数据用MrBayes v3.2.7a（http://nbisweden.github.io/MrBayes/）软件, 选用GTR+I+G模型，Ngammacat设置为5，statefreqpr、revmat、pinvar、shapepr根据jModelTest软件寻找的最佳模型来设置，其余参数默认，构建贝叶斯进化树，分析结果如下图所示：
采用全基因组做进化树分析，将环形序列设置相同起点，物种间序列用MAFFT v7.427软件（--auto模式）进行多序列比对，
将比对好的数据用MrBayes v3.2.7a（http://nbisweden.github.io/MrBayes/）软件,
选用GTR+I+G模型，Ngammacat设置为5，statefreqpr、revmat、pinvar、shapepr根据jModelTest软件寻找的最佳模型来设置，
其余参数默认，构建贝叶斯进化树，分析结果如下图所示
2.
采用全基因组做进化树分析，将环形序列设置相同起点，物种间序列用MAFFT软件（v7.427，--auto模式）进行多序列比对，
将比对好的序列首尾连接，用trimAl(v1.4.rev15)修剪，然后使用MrBayes v3.2.7，马尔科夫链蒙特卡洛方法（Markov Chain Monte Carlo, MCMC），
迭代运算1百万代，每 100 代取样一次。运行结果最初的 25%系统树被删除（burn-in），最终得到多数一致树（majority-rule consensus tree）。
分析结果如下图所示：
全基因组最大似然法
采用全基因组做进化树分析，将环形序列设置相同起点，物种间序列用MAFFT v7.427软件（--auto模式）进行多序列比对，将比对好的数据用trimAl v1.4.rev15修剪，然后使用RAxML v8.2.10（https://cme.h-its.org/exelixis/software.html）软件，选用GTRGAMMA模型，rapid Bootstrap分析，bootstrap=1000，构建最大似然进化树，分析结果如下图所示：
采用全基因组做进化树分析，将环形序列设置相同起点，物种间序列用MAFFT v7.427软件（--auto模式）进行多序列比对，
将比对好的数据用trimAl v1.4.rev15修剪，
然后使用RAxML v8.2.10（https://cme.h-its.org/exelixis/software.html）软件，
选用GTRGAMMA模型，rapid Bootstrap分析，bootstrap=1000，
构建最大似然进化树，分析结果如下图所示：
全基因组邻接法
采用全基因组做进化树分析，将环形序列设置相同起点，
物种间序列用MAFFT v7.427软件（--auto模式）进行多序列比对，
将比对好的数据用trimAl v1.4.rev15修剪，
然后用rapidnj v2.3.3（https://github.com/somme89/rapidNJ）软件，
进行1000次bootstrap重采样，构建邻接法进化树，分析结果如下图所示：
CDS贝叶斯法
采用共有CDS做进化树分析，将每组CDS序列用MAFFT v7.427软件（--auto模式）进行多序列比对，将比对好的CDS序列首尾相接（每个物种的CDS序列串联在一起）后
用MrBayes v3.2.7a（http://nbisweden.github.io/MrBayes/）软件，选用GTR+I+G模型，Ngammacat设置为5，statefreqpr、revmat、pinvar、shapepr根据jModelTest软件寻找的最佳模型来设置，其余参数默认，其余参数默认，构建贝叶斯进化树，分析结果如下图所示：
采用共有CDS做进化树分析，将每组CDS序列用MAFFT v7.427软件（--auto模式）进行多序列比对，
将比对好的CDS序列首尾相接（每个物种的CDS序列串联在一起）后
用MrBayes v3.2.7a（http://nbisweden.github.io/MrBayes/）软件，
选用GTR+I+G模型，Ngammacat设置为5，statefreqpr、revmat、pinvar、shapepr根据jModelTest软件寻找的最佳模型来设置，
其余参数默认，构建贝叶斯进化树
CDS最大似然法
采用共有CDS做进化树分析，将每组CDS序列用MAFFT v7.427软件（--auto模式）进行多序列比对，将比对好的CDS序列首尾相接（每个物种的CDS序列串联在一起）后
用RAxML v8.2.10（https://cme.h-its.org/exelixis/software.html）软件，选用GTRGAMMA模型，rapid Bootstrap分析，bootstrap=1000，构建最大似然进化树，分析结果如下图所示：
采用共有CDS做进化树分析，将每组CDS序列用MAFFT v7.427软件（--auto模式）进行多序列比对，
将比对好的CDS序列首尾相接（每个物种的CDS序列串联在一起）后
用RAxML v8.2.10（https://cme.h-its.org/exelixis/software.html）软件，
选用GTRGAMMA模型，rapid Bootstrap分析，bootstrap=1000，
构建最大似然进化树
CDS邻接法
采用共有CDS做进化树分析，将每组CDS序列用MAFFT v7.427软件（--auto模式）进行多序列比对，
将比对好的CDS序列首尾相接（每个物种的CDS序列串联在一起）后
用rapidnj v2.3.3（https://github.com/somme89/rapidNJ）软件，
进行1000次bootstrap重采样，构建邻接法进化树
最大似然法参数
Note
raxmlHPC-PTHREADS-SSE3 -f a -s Pecans.trim.aln -T 8 -# 1000 -m GTRGAMMA -x 12345 -p 12345 -n genome -w ../raxml 1>../raxml/RAxML.log
-m GTRCAT: GTR近似法，对每个站点的替代率进行优化，并将这些单独的替代率分为-c指定的比率类别。
GTRCAT：GTR近似法，对每个站点的替代率进行优化，并将这些单独的替代率分为由-c指定的比率类别。这只是一个
这只是GTRGAMMA的一种变通方法，所以要确保不要根据GTRCAT的似然值来比较其他拓扑结构。
似然值进行比较。因此，你不能将GTRCAT与-f e（树评估）结合使用，也不能将GTRCAT与多个分析结合使用。
不能与原始排列的多重分析（-# |-N）选项结合使用。这是由于
这是因为作者假设你想根据可能性来比较树，如果你在原始排列上进行多
运行时，作者假定你想根据可能性来比较树木。如果你指定了例如-m GTRCAT和-# 10，程序会自动
使用GTRMIX（见下文）。
-m GTRMIX: 这个选项将使RAxML在GTRCAT下进行树状推断（搜索一个好的拓扑结构）。
GTRCAT。当分析完成后，RAxML将把它的模型切换到GTRGAMMA，并评估最终的
在GTRGAMMA下评估最终的树形拓扑结构，使其产生稳定的似然值。
-m GTRGAMMA: GTR(General Time Reversible)核苷酸替换模型[15]
GTRGAMMA 的实现采用了 4 个离散的速率类别，在速度和准确度之间进行了可接受的权衡。
准确度之间可以接受的权衡。请注意，由于性能的原因，这是硬编码的，也就是说，离散速率类别的数量不能被改变。
即用户不能改变离散率类别的数量。
-m GTRCAT GAMMA。用特定地点的进化率进行树的推断。然而，这里的速率是
使用4个离散的GAMMA速率进行分类，遵循Yang[17]
但我还是不喜欢这个主意（见第6节的讨论）。
-m GTRMIXI: 与GTRMIX相同，但要估计不变位点的比例。
-m GTRCAT GAMMAI: 与GTRCAT_GAMMA相同，但要估计出可变位点的比例。
公司原有教程
贝叶斯进化树构建
1.对序列进行多序列比对生成
软件mafft,clustalo,muscle等
合并之后用
mafft
2.多序列比对过滤
软件trimal
trimal -in <inputfile> -out <outputfile> -(other options)
3.格式转换
这之前都是fasta
perl /share/nas6/xul/program/chloroplast/phytree/aln2nexus.pl
4.作树
/share/nas6/xul/soft/mrbayes/MrBayes-3.2.7/src/mb
execute example.nex # 导入nex文件 /share/nas6/zhangxq/project/xianliti/GP-20190626-1437_fujiannonglin_hanxiaohong_1sample_hongtouwujing_xianliti/analysis/MrBayes/tree.nexus
lset nst=6 rates=invgamma # 设置进化模型参数.本例中设定数据为DNA数据. 蛋白的需要将DNA改成protein lset nst=6 Rates=invgamma Nucmodel=Protein
mcmc ngen=1000000 samplefreq=100 printfreq=100 diagnfreq=5000 // mcmc ngen=10000 samplefreq=10 printfreq=10 diagnfreq=5000 # Average standard deviation of split frequencies < 0.01 n
# ngen则是运行的长度，默认1,000,000次；samplefreq则是取样频率，每隔多少次运行次数取一次样；printfreq是打印频率，即每运行多少次将打印一行结果到屏幕上，默认为500；
diagnfreq则代表每运行多少次分析一次结果，得出 Average standard deviation of split frequencies，默认是5,000.
如果在设定的代数运行完毕后，给出的 Average standard deviation of split frequencies的值小于0.01，则根据提示输入‘no'来停止运行，反之则输入'yes'继续运行直到满足其值小于0.01为止。
If you are intersted mainly in the well-supported parts of the tree, a standard deviation below 0.05 may be adequate.
no
sump burnin=2500 // 250 # (ngen/samplefreq)*0.25 # 使用sump来对参数值进行归纳。设置的burnin值为 (ngen / samplefreq) * 0.25 。程序给出一个概括的表，要确保PSRF一列中的值接近 1.0，否则需要运行该多的代数。
sumt burnin=2500 // 250 # 使用sumt来构树。burnin值和前一个相同
5.MEGA展示
针对打不开的情况可以使用figtree
显示置信度的方法，点击branchlabels标签，display调成prob（percent），下面也选为percent然后保存为nwk
tree标签下的transform勾选之后可以对齐。
最好能保存一个nwk文件，但好像有一些小问题。保存一个pdf文件。
过长的名称会导致遮盖，右击选择轮廓，扩大一下，然后右击GPU预览回到正常，调整位置，右侧菜单释放蒙版