使用KaKs_calculator计算Ka/Ks

安装kakscalculator2

chmod 757 ./KaKs_Calculator2.0/src/AXTConvertor
chmod 757 ./KaKs_Calculator2.0/bin/Linux/KaKs_Calculator

添加环境变量

echo 'export PATH=/share/software/KaKs_Calculator2.0/bin/Linux:$PATH' >>~/.bashrc
echo 'export PATH=/share/software/KaKs_Calculator2.0/src:$PATH' >>~/.bashrc
source ~/.bashrc

安装序列比对软件mafft 下载地址

wget https://mafft.cbrc.jp/alignment/software/mafft-7.471-with-extensions-src.tgz
tar -zxvf mafft-7.471-with-extensions-src.tgz
cd mafft-7.471-with-extensions
cd core
vim Makefile
#手动修改prefix和bindir

我的是修改为下面的路径
PREFIX = /share/softwares/mafft
BINDIR = /share/softwares/mafft/bin

make clean
make 
make install
echo `export PATH=/share/softwares/mafft/bin:$PATH` >>~/.bashrc
source ~/.bashrc

mafft官网
 mafft使用

下载安装ParaAT https://bigd.big.ac.cn/tools/paraat

ParaAT是整合工具，可以调用比对工具和KaKs。

准备输入文件

test.cds #每个基因最长的转录本的DNA序列
test.pep #每个基因最长的蛋白序列
cds和pep的获取方式，可以使用jcvi的工具包参考
proc 数字（小于CPU实际数量的数字）

生成 `test.homolog`是共线性分析的基因对参考

#构建数据库
makeblastdb -in test.pep -dbtype prot -parse_seqids -out test_sinensis
blastp -query test.pep -out test_sinensis.blast -db test_sinensis -outfmt 6 -evalue 1e-10 -num_threads 8 -qcov_hsp_perc 50.0 -num_alignments 5 2> blastp.log &
#MCScanX 需要2个输入文件test_sinensis.gff和test_sinensis.blast
MCScanX test_sinensis -g -3 -e 1e-10
#MCScanX输出的`test_sinensis.collinearity`即是共线性基因对
cat test_sinensis.collinearity|grep -v "^#"|cut -f 2,3 >test.homolog

开始分析

ParaAT.pl -h test.homolog -n test.cds -a test.pep -p proc -m mafft -f axt -g -k -o result_dir
#查看KaKs的详情  cat result_dir/*.kaks | cut -f 1,3,4,5 | grep -v 'Sequence' |less -S

#把KaKs输出到文件KaKs.txt
cat result_dir/*.axt.kaks | cut -f 1,3,4,5 | grep -v 'Sequence' >test.KaKs.txt
cat test.KaKs.txt| tr "\t" "," |sed '1i Sequence,Ka,Ks,Ka/Ks' >test.KaKs.csv

ParaAT.pl参数讲解：
-h, 同源基因名称文件
-n, 指定核酸序列文件
-a, 指定蛋白序列文件
-p, 指定多线程文件数字
-m, 指定比对工具(clustalw2 | t_coffee | mafft | muscle)，多选一
-g, 去除比对有gap的密码子
-k, 用KaKs_Calculator 计算kaks值
-o, 输出结果的目录
-f, 输出比对文件的格式
*** 也可通过-f参数得到其他软件分析ka/ks所需的格式

可视化（箱线图和小提琴图）

setwd("e:/****/GeneDuplication")
library("ggplot2")

data_kaks <- read.csv("test.KaKs.csv")
data_kaks$Item <- "normal" #如果有多组，可以在此列设置分组信息
ggplot(data=data_kaks, aes(x=Item, y=Ka.Ks, fill=Item)) + 
  geom_violin(alpha=0.8,width=1)+ guides(fill=F)+xlab(' ')+ylab('Ka/Ks')+
  geom_boxplot(alpha=0.5,varwidth = TRUE)+
  labs(title = "Distribution of Ka/Ks distances", size = 1.5)
ggsave('distribute of KaKs.pdf',dpi=300)
 ggplot(data=data_kaks, aes(x=Ka,group=Item)) + 
  geom_density(alpha=0.4,aes(color = Item))+ xlab('Ka value')+ylab('Density')+
 labs(title = "Distribution of Ka distances", size = 1.5)+guides()
ggsave('Distribution of Ka',dpi=300)
head(data_kaks)
#Sequence        Ka        Ks    Ka.Ks   Item
#1 Sl01G000030-Sl15G011500 0.0721061 0.0846677 0.851637 normal
#2 Sl01G000130-Sl15G011440 0.0162336 0.1356720 0.119654 normal
#3 Sl01G000170-Sl15G011380 0.0224351 0.1382140 0.162322 normal
#4 Sl01G000200-Sl15G011320 0.1947300 0.2629070 0.740682 normal
#5 Sl01G000330-Sl15G011250 0.0268044 0.1196890 0.223950 normal
#6 Sl01G000440-Sl15G011120 0.3277060 0.7499520 0.436969 normal

最后编辑于：2021.01.12 16:53:24