一.基本步骤
ParaAT.pl计算得到file.axt文件，如果加上-k参数，就会使用KaKs_Calculator的MA方法来计算（速度很慢），为了加快计算速度，可以使用YN的方法来计算。于是先用ParaAT.pl得到file.axt文件，然后使用KaKs_Calculator来计算。
ParaAT.pl -h test.homologs -n test.cds -a test.pep -p proc -m muscle -f axt -g -k -o result_dir #proc文件必须与输出位置在同一个目录下，不然会报错
# -h 指定同源基因对列表文件:
# -n 指定核酸序列文件
# -a 指定蛋白序列文件
# -p 指定多线程文件(必须是一个文件而不是数字)
# -m 指定比对工具：muscle这个序列比对软件有问题，可能会导致比对空文件的出现，mafft发现正常。主要是muscle不同版本的命令不同。
[clustalw2 | t_coffee | mafft | muscle , optional]
# -g 去除比对有gap的密码子
# -k 用 KaKs_Calculator 计算kaks值 ： 使用这个参数会使用使用多个方法进行计算，导致速度变慢
# -o 输出结果的目录，两次运行的输出目录不能相同
# -f 输出比对文件的格式：[axt | fasta | paml | codon | clustal, optional], default = fasta
# ParaAT.pl -help
KaKs_Calculator计算kaks值
find output/ -name '*.axt' -type f -print0 | parallel -0 -j 40 KaKs_Calculator -i {} -o ${kaks_outdir}/{/.}.kaks -m YN
上面的方法其实可以优化，只需要将
ParaAT.pl脚本进行适当修改即可。这样修改之后就可以一行命令执行了。
my $cmd = $KaKsCMD." -i ".$file.".cds_aln.axt -o ".$file.".cds_aln.axt.kaks >> msg.kaks";
# 在这行代码里面添加方法参数（-m YN）：
my $cmd = $KaKsCMD." -i ".$file.".cds_aln.axt -o ".$file.".cds_aln.axt.kaks -m YN >> msg.kaks";
run_kaks(){
cat ${qry}.cds ${ref}.cds > "${qry_basename}.cds"
cat ${qry}.pep ${ref}.pep > "${qry_basename}.pep"
echo "${thread}" > "threadNum.txt"
if ! [ -e "${qry_basename}-${ref_basename}.pair" ];then
ls /home/bqxiao/data5/01.pan/19.duplicate/2.analyse/data/alignment/${qry_basename}_${ref_basename}.blast | xargs -I {} awk -F "\t" "{ if (\$1 != \$2) print \$1\"\t\"\$2 }" {} > ${qry_basename}-${ref_basename}.pair
fi
nohup ParaAT.pl -h "${qry_basename}-${ref_basename}.pair" \
-n "${qry_basename}.cds" \
-a "${qry_basename}.pep" \
-p "threadNum.txt" \
-o "output" -m mafft -g -f axt -k > ParaAT.log &
}
ref="/home/bqxiao/data5/01.pan/19.duplicate/data/cacao/temp/cacao"
qrys="/home/bqxiao/data5/01.pan/01.newdata/cds-pep/Aralia_elata_Seem
/home/bqxiao/data5/01.pan/01.newdata/cds-pep/Coriandrum_sativum
"
ref_basename=$(basename ${ref})
thread=10
for qry in ${qrys};do
qry_basename=$(basename ${qry})
if ! [ -e "${qry_basename}" ];then
mkdir "${qry_basename}"
fi
cd ${qry_basename}
run_kaks
cd ..
done
二.使用这个脚本时可能遇到的问题
输出目录必须没有创建
这个脚本会自动修改基因名称(将基因名称后面的 点 去除)。目前已做修改，可以返回完整的基因名称，而不会因为名称里面有版本号而去除。主要是这两个函数，GenerateSequenceFiles和ReadSeqs。
-
-h 对应的参数必须是基因对列表，不能直接使用blast文件。
ls alignment/${qry_basename}_${ref_basename}.blast | xargs -I {} awk -F "\t" "{ if (\$1 != \$2) print \$1\"\t\"\$2 }" {} > output/${qry_basename}/${qry_basename}-${ref_basename}.pair
运行ParaAT.pl脚本的运行位置必须和输入文件及输出目录处在同一个目录下，不然运行到序列比对时，将会出现bug(无法找到多线程文件proc)
-
若是出现检查完cds文件、pep文件以及基因对文件之后，程序直接结束。
Generating homologous group files for alignments: 964 from 964 groups这两个数字不相等，说明cds和pep文件以及基因对文件中名称可能不一致，或者fasta文件中有注释信息，比如：>gene1 [sequence in Human]...
-
几个程序的比对速度以及准确性：
- 比对速度：Muscle>MAFFT>ClustalW>T-Coffee
- 比对准确性：MAFFT>Muscle>T-Coffee>ClustalW
这个脚本中的muscle命令如下：
elsif ($MSACMD =~ /muscle/) {
$cmd = "muscle -quiet -in ".$file.".pep -out ".$file.".pep_aln >> msg.msa";
}
# 显然这是muscle3的输入格式，而muscle5的输入格式如下：
# muscle -align input.fa -output aln.afa
- 所输入的cds文件和pep文件基因名称必须一致，而且都没有注释信息。可使用下面的命令过滤文件。
# Usage:
# ls *.pep | xargs -I {} filter.py {}
# ls *.cds | xargs -I {} filter.py {}
def process_file(input_filename, output_filename):
with open(input_filename, 'r', encoding='utf-8') as file:
lines = file.readlines()
processed_lines = []
for line in lines:
if line.startswith('>'):
# 去除>开头中中的注释部分内容，只保留基因名称
line = line.split()[0] + "\n"
processed_lines.append(line)
else:
processed_line = line.replace('.', '').replace('*', '')
if processed_line.strip(): # 检查行是否只剩换行符
processed_lines.append(processed_line)
with open(output_filename, 'w', encoding='utf-8') as file:
file.writelines(processed_lines)
import sys
if len(sys.argv) != 2:
print("Usage: python3 filter.py input_filename")
sys.exit(1)
input_filename = sys.argv[1]
output_filename = input_filename
process_file(input_filename, output_filename)
