前期准备

• 不同物种的蛋白和cds序列：
  os.pep, os.cds, sb.pep, sb.cds
• 依赖程序：
  blast+, clustalw2, ParaAT, KaKs_Calculator
  注：所依赖程序需提前安装好，并添加到环境变量中
• 所用脚本：
  axt2one-line.py, calculate_4DTV_correction.pl

获得不同物种之间的同源序列

对参考物种的蛋白序列构建索引

makeblastdb -in sb.pep -dbtype prot

将目标物种的蛋白序列与参考物种进行比对，并保留最优匹配结果

blastp -query os.pep -db sb.pep -evalue 1e-5 -max_target_seqs 1 -num_threads 12 -out os2sb.blastp_out.m6 -outfmt 6

提取最优匹配的同源序列基因对

cut -f1-2 os2sb.blastp_out.m6|sort|uniq > os2sb.homolog

合并目标物种和参考物种的蛋白序列和cds序列

cat os.cds sb.cds >os_sb.cds
cat os.pep sb.pep >os_sb.pep

使用ParaAT程序将蛋白序列比对结果转化为cds序列比对结果

# 新建proc文件, 设定12个线程
echo "12" >proc
# -m参数指定多序列比对程序为clustalw2，-p参数指定多线程文件，-f参数指定输出文件格式为axt
ParaAT.pl -h os2sb.homolog -n os_sb.cds -a os_sb.pep -m clustalw2 -p proc -f axt -o os2sb_out

计算kaks值和4dtv值

cd os2sb_out
# 使用KaKs_Calculator计算kaks值, -m参数指定kaks值的计算方法为YN模型
for i in `ls *.axt`;do KaKs_Calculator -i $i -o ${i}.kaks -m YN;done
# 将多行axt文件转换成单行
for i in `ls *.axt`;do axt2one-line.py $i ${i}.one-line;done
# 使用calculate_4DTV_correction.pl脚本计算4dtv值
ls *.axt.one-line|while read id;do calculate_4DTV_correction.pl $id >${id%%one-line}4dtv;done
# 合并所有同源基因对的4dtv
for i in `ls *.4dtv`;do awk 'NR>1{print $1"\t"$2}' $i >>all-4dtv.txt;done
# 合并所有同源基因对的kaks值
for i in `ls *.kaks`;do awk 'NR>1{print $1"\t"$3"\t"$4"\t"$5}' $i >>all-kaks.txt;done
# 给结果文件排序并去冗余
sort all-4dtv.txt|uniq >all-4dtv.results
sort all-kaks.txt|uniq >all-kaks.results
# 清除中间文件
rm *one-line
rm all-4dtv.txt
rm all-kaks.txt
# 将kaks结果文件和4dtv结果文件进行合并
join -1 1 -2 1 all-kaks.results all-4dtv.results >all-results.txt
# 给结果文件添加标题
sed -i '1i\Seq\tKa\tKs\tKa/Ks\t4dtv_corrected' all-results.txt

使用calc_kaks_4dtv.sh脚本一步计算kaks值和4dtv值

nohup sh calc_kaks_4dtv.sh os sb &

查看脚本文件

#!/bin/bash
set -e
set -u

if [ $# -lt 2 ];then
    echo "Two parameters needed, but only $# given!"
    echo "Usage: sh calc_kaks_4dtv.sh <Species1> <Species2>"
exit 1;
fi

Species1=$1
Species2=$2

# 对参考物种的蛋白序列构建索引
makeblastdb -in ${Species2}.pep -dbtype prot
# 将目标物种的蛋白序列与参考物种进行比对，并保留最优匹配结果
blastp -query ${Species1}.pep -db ${Species2}.pep -evalue 1e-5 -max_target_seqs 1 -num_threads 10 -out ${Species1}2${Species2}.blastp_out.m6 -outfmt 6
# 提取最优匹配的同源序列基因对
cut -f1-2 ${Species1}2${Species2}.blastp_out.m6|sort|uniq > ${Species1}2${Species2}.homolog
# 合并目标物种和参考物种的蛋白序列和cds序列
cat ${Species1}.cds ${Species2}.cds >${Species1}_${Species2}.cds
cat ${Species1}.pep ${Species2}.pep >${Species1}_${Species2}.pep
# 使用ParaAT程序将蛋白序列比对结果转化为cds序列比对结果
ParaAT.pl -h ${Species1}2${Species2}.homolog -n ${Species1}_${Species2}.cds -a ${Species1}_${Species2}.pep -p proc -m clustalw2 -f axt -o ${Species1}2${Species2}_out
cd ${Species1}2${Species2}_out
# 使用KaKs_Calculator计算kaks值
for i in `ls *.axt`;do KaKs_Calculator -i $i -o ${i}.kaks -m YN;done
# 将多行axt文件转换成单行
for i in `ls *.axt`;do axt2one-line.py $i ${i}.one-line;done
# 使用calculate_4DTV_correction.pl脚本计算4dtv值
ls *.one-line|while read id;do calculate_4DTV_correction.pl $id >${id%%one-line}4dtv;done
# 合并所有同源基因对的4dtv值
for i in `ls *.4dtv`;do awk 'NR>1{print $1"\t"$2}' $i >>all-4dtv.txt;done
# 合并所有同源基因对的kaks值
for i in `ls *.kaks`;do awk 'NR>1{print $1"\t"$3"\t"$4"\t"$5}' $i >>all-kaks.txt;done
# 排序并去冗余
sort all-4dtv.txt|uniq >all-4dtv.results
sort all-kaks.txt|uniq >all-kaks.results
# 清除中间文件
rm *one-line
rm all-4dtv.txt
rm all-kaks.txt
# 将kaks结果文件和4dtv结果文件进行合并
join -1 1 -2 1 all-kaks.results all-4dtv.results >all-results.txt
# 给结果文件添加标题
sed -i '1i\Seq\tKa\tKs\tKa/Ks\t4dtv_corrected' all-results.txt

脚本运行完后会生成以下结果文件

image.png

最终得到的结果在os2sb_out文件夹中的all-results.txt文件中

image.png

绘制kaks值和4dtv值的密度分布图

setwd("/Users/Davey/Desktop")
library(ggplot2)
library(patchwork)
data <- read.table("all-results.txt",header=T,check.names=F)
head(data)
                       Seq        Ka       Ks     Ka/Ks 4dtv_corrected
1 Os01g0276600-Sb03g022840 0.2954040 2.146660 0.1376110      0.3852584
2 Os01g0276700-Sb03g022860 0.0372221 0.586840 0.0634280      0.2272330
3 Os01g0276800-Sb03g022870 0.1293250 0.560683 0.2306560      0.2456600
4 Os01g0721900-Sb03g158830 0.0769183 0.849214 0.0905759      0.1966413
5 Os01g0723100-Sb03g158910 0.1109170 1.198870 0.0925176      0.4060341
6 Os01g0723200-Sb03g158920 0.0398003 0.941168 0.0422882      0.1695289
p1 <- ggplot(data,aes(Ks))+ geom_line(stat="density",color="red")+theme_classic()+theme(axis.title = element_text(size=24),axis.text = element_text(size = 22,color = "black"))
p2 <- ggplot(data,aes(Ka/Ks))+ geom_line(stat="density",color="blue")+theme_classic()+theme(axis.title = element_text(size=24),axis.text = element_text(size = 22,color = "black"))
p3 <- ggplot(data,aes(`4dtv_corrected`))+ geom_line(stat="density",color="orange")+theme_classic()+theme(axis.title = element_text(size=24),axis.text = element_text(size = 22,color = "black"))
p1 + p2 + p3

image.png