用mcmctree做分歧时间估计的时候遇到报错:
MCMCTREE in paml version 4.9, March 2015
Reading options from mcmctree.ctl..
finetune is deprecated now.
Reading master tree.
(Bmor, ((Achi, Rpro), Nlug));
Reading sequence data.. 3 loci
*** Locus 1 ***
ns = 4 ls = 621319
Reading sequences, sequential format..
Reading seq # 1: Achi
Error in sequence data file: B at 621301 seq 1.
Make sure to separate the sequence from its name by 2 or more spaces.
问题是序列不符合phylip格式,每个名字没有区分开
1.使用以下命令行将orthofinder处理好的多序列比对结果,fasta格式改为phylip格式,失败
cat SpeciesTreeAlignment.fa |tr '\n' '\t'|sed 's/>/\n/g' |sed 's/\t/ /'|sed 's/\t//g'| awk 'NF > 0' >Aa.phy.tmp
awk '{print " "NR" "length($2)}' supergene.phy.tmp|tail -n 1 | cat - supergene.phy.tmp > Aa.phy
- 使用
cat SpeciesTreeAlignment.fa |tr '\n' '\t'|sed 's/>/\n>/g' |sed 's/\t/\n/'|sed 's/\t//g'| awk 'NF > 0' > Aa1.phy.tmp
转换为正常的fasta(去掉换行等)
使用R包将fasta转换为phylip格式
library(devtools)
library(ape)
library(phylotools)
data <- read.fasta("Aa1.phy.tmp")
dat2phylip(data, outfile = "out.phy")
还是得到相同报错
报错原因,查找可能是蛋白序列里带有U,U不在常见密码子当中,所以有些软件不识别会报错
sed 's/U/X/g' SpeciesTreeAlignments.fa > STA_delU.fa
再用得到的这个去除U fasta文件,执行以上命令行,得到phylip格式文件。序列个数没有问题。注意,千万不要把U替换为空,会影响文件序列长度。应该替换为X。
运行继续报错,但是距离成功不远了。
ns = 4 ls = 621319
Reading sequences, sequential format..
Reading seq # 4: Rpro
Sequences read..
Counting site patterns.. 0:00
56693 patterns at 621319 / 621319 sites (100.0%), 0:01
Counting frequencies..
56693 patterns, messy
*** Locus 2 ***
Error: seq err1: EOF.
seq file is not paml/phylip format. Trying nexus format.
原来问题是,软件自带测试数据是3组数据,而我们的数据只有一组!
所以非常简单,修改mcmctree.ctl。改为ndata=1
ndata = 1
seqtype = 2 * 0: nucleotides; 1:codons; 2:AAs
usedata = 3 * 0: no data; 1:seq like; 2:use in.BV; 3: out.BV
clock = 3 * 1: global clock; 2: independent rates; 3: correlated rates
RootAge = <1.0 * safe constraint on root age, used if no fossil for root.
运行成功。