安裝kakscalculator2

kakscalculator2下載地址
解壓后

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

添加環(huán)境變量

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

安裝序列比對(duì)軟件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
#手動(dòng)修改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官網(wǎng)
mafft使用

下載安裝ParaAT https://bigd.big.ac.cn/tools/paraat

ParaAT是整合工具，可以調(diào)用比對(duì)工具和KaKs。

準(zhǔn)備輸入文件

• test.cds #每個(gè)基因最長(zhǎng)的轉(zhuǎn)錄本的DNA序列
• test.pep #每個(gè)基因最長(zhǎng)的蛋白序列
  cds和pep的獲取方式，可以使用jcvi的工具包 參考
• proc 數(shù)字（小于CPU實(shí)際數(shù)量的數(shù)字）

生成 test.homolog是共線性分析的基因?qū)?參考

#構(gòu)建數(shù)據(jù)庫(kù)
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個(gè)輸入文件test_sinensis.gff和test_sinensis.blast
MCScanX test_sinensis -g -3 -e 1e-10
#MCScanX輸出的`test_sinensis.collinearity`即是共線性基因?qū)?cat test_sinensis.collinearity|grep -v "^#"|cut -f 2,3 >test.homolog

開(kāi)始分析

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參數(shù)講解：
-h, 同源基因名稱(chēng)文件
-n, 指定核酸序列文件
-a, 指定蛋白序列文件
-p, 指定多線程文件 數(shù)字
-m, 指定比對(duì)工具(clustalw2 | t_coffee | mafft | muscle)，多選一
-g, 去除比對(duì)有g(shù)ap的密碼子
-k, 用KaKs_Calculator 計(jì)算kaks值
-o, 輸出結(jié)果的目錄
-f, 輸出比對(duì)文件的格式
*** 也可通過(guò)-f參數(shù)得到其他軟件分析ka/ks所需的格式

可視化（箱線圖和小提琴圖）

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

data_kaks <- read.csv("test.KaKs.csv")
data_kaks$Item <- "normal" #如果有多組，可以在此列設(shè)置分組信息
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