WGD（全基因组复制）分析——Ka/Ks及4Dtv值计算
一、4Dtv（四倍简并位点颠换率）的定义
4DTV stands for four-fold synonymous (degenerative) third-codon transversion. It represents a transversion in the third nucleotide position within four codons that does not result in a change in corresponding amino acid identity within the protein it codes for. Such an estimate of synonymous mutation rate within a transcribed region of a gene but not in region that experiences selection is a conserved means of estimating divergence within the more recent evolutionary past. Distances corresponding to the 'salicoid' whole-genome duplication events were delineated based on discrete peaks in 4DTV distributions. 4Dtv (transversion of four-fold degenerate site) values of each block were calculated using the sum of transversion of four-fold degenerate sites divided by the sum of four-fold degenerate sites.
其生物学意义为：共线性区域内基因对的4Dtv值可以反映在进化过程中的物种相对分化事件以及全基因组复制事件。
Sanwen huang et al., 2009
二、关于Ka/Ks
Ka/Ks表示的是非同义替换（Ka）和同义替换（Ks）之间的比例，这一比值可以判断编码该蛋白的基因是否遭受了选择压力。同义突变表示，突变并不影响氨基酸序列，进而不会影响蛋白结构与功能。而非同义突变则会影响氨基酸序列，可能会使其结构和功能发生改变，可能会遭受自然选择。一般我们认为，同义突变不受自然选择，而非同义突变会遭受自然选择作用。在生物进化分析中，知晓物种的同义突变和非同义突变发生的速率是非常有意义的。同义突变频率即为Ks值，非同义突变频率即为Ka值，非同义突变率与同义突变率的比值即为Ka/Ks值。若Ka/Ks > 1，则认为存在正选择效应（positive selection）；若Ka/Ks = 1，则认为存在中性选择效应；如若Ka/Ks < 1，则认为存在负选择效应，即纯化效应或净化选择（purifying selection）。
三、4Dtv及Ka/Ks值的计算
1、数据准备
Arabidopsis_thaliana.cds、Arabidopsis_thaliana.pep、Arabidopsis_thaliana.gff3；Citrus_grandis.cds、Citrus_grandis.pep、Citrus_grandis.gff3；Citrus_sinensis.cds、Citrus_sinensis.pep、Citrus_sinensis.gff3；Malus_domestica.cds、Malus_domestica.pep、Malus_domestica.gff3；
2、数据处理
取file.pep中最长转录本序列，去冗余。
同一个基因可能存在多个转录本
python3 extract_longest_iso.py -i Citrus_sinensis.fa -o L_Citrus_sinensis.fa
去冗余前44275条序列，去冗余后23394条序列；
去冗余后
python脚本如下
extract_longest_iso.py
3、获取共线性基因对
（1）对蛋白序列构建索引
makeblastdb -in Citrus_sinensis.pep -dbtype prot -parse_seqids -out Citrus_sinensis
（2）blastp比对
nohup blastp -query Citrus_sinensis.pep -out Citrus_sinensis.blast -db Citrus_sinensis -outfmt 6 -evalue 1e-10 -num_threads 8 -qcov_hsp_perc 50.0 -num_alignments 5 2> blastp.log &
（3）gff3文件处理
grep '\smRNA\s' Citrus_sinensis.gff3 | awk '{print $1"\t"$4"\t"$5"\t"$9}' | awk -F 'ID=' '{print $1$2}' | awk -F 'Parent=' '{print $1}' | awk '{print $1"\t"$4"\t"$2"\t"$3}' | awk -F ';' '{print $1"\t"$3}' | awk '{print $1"\t"$2"\t"$3"\t"$4}' > Citrus_sinensis.gff
gff3格式文件
仅保留染色体、转录本ID及坐标信息
（4）运行MCScanX(.blast & .gff文件 )
MCScanX是一款分析基因组内或者基因组间的共线性区块的软件，它利用种内或种间蛋白质blastp比对结果，再结合编码这些蛋白的基因在基因组中的位置坐标，得到种内或种间基因组的共线性区块。MCScanX软件安装及详细使用参见官网，安装和使用都比较友好。http://chibba.pgml.uga.edu/mcscan2/#tm
运行MCScanX，获取共线性基因对
MCScanX Citrus_sinensis -g -3 -e 1e-10
得到 Citrus_sinensis.collinearity、Citrus_sinensis.tandem文件及Citrus_sinensis.html文件夹，其中我们需要的信息就在Citrus_sinensis.collinearity结果文件中。
Citrus_sinensis.collinearity结果文件
（5）提取共线性基因对
cat Citrus_sinensis.collinearity | grep "Cs" | awk '{print $3"\t"$4}' > Citrus_sinensis.homolog
Citrus_sinensis.homolog
4、Ka、Ks及4Dtv值计算
（1）准备软件及脚本
需要用到的软件：KaKs_Calculator2.0及ParaAT；软件安装参考：http://cbb.big.ac.cn/Software
其中需要注意到一个问题，不少同学在安装KaKs_Calculator2.0时会报错，当然我也碰到了，查了好久才解决，为了避免再次踩坑，特贴出解决方法；
### make: *** [KaKs_Calculator] error.
解压安装包后，在“base.h”文件中最前面添加一行内容：
#include <string.h>
然后保存、退出，再“make”命令安装即可。
需要用到的脚本：calculate_4DTV_correction.pl；axt2one-line.py
https://github.com/scbgfengchao/4DTv/blob/master/axt2one-line.py
这个脚本其实存在一个问题，简书中也有其他作者指出来过，需要将脚本中密码子表“my %codons=”中“U”改成“T”；这里RNA codon 和 DNA codon混在一起了，是有问题的，在此统一用DNA codon。
calculate_4DTV_correction.pl
（2）数据准备
Citrus_sinensis.homolog、Citrus_sinensis.cds、Citrus_sinensis.pep；Citrus_grandis.homolog、Citrus_grandis.cds、Citrus_grandis.pep；Arabidopsis_thaliana.homolog、Arabidopsis_thaliana.cds、Arabidopsis_thaliana.pep；Malus_domestica.homolog、Malus_domestica.cds、Malus_domestica.pep；
（3）使用ParaAT程序将蛋白序列比对结果转化为cds序列比对结果
# 新建proc文件, 设定12个线程
echo "12" >proc
# -m参数指定多序列比对程序为clustalw2，-p参数指定多线程文件，-f参数指定输出文件格式为axt
nohup ParaAT.pl -h Citrus_sinensis.homolog -n Citrus_sinensis.cds -a Citrus_sinensis.pep -m clustalw2 -p proc -f axt -o Citrus_sinensis_out 2> ParaAT.log &
此步需注意，file.homolog、file.cds、file.pep三个文件的基因ID需保持一致，且file.cds及file.pep为fasta格式，即“>”后面只接基因ID号，否则会报错，如下：
Citrus_sinensis.homolog
Citrus_sinensis.cds
Citrus_sinensis.pep
（4）使用KaKs_Calculator计算ka、ks值, -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
（5）使用calculate_4DTV_correction.pl脚本计算4dtv值
ls *.axt.one-line|while read id;do calculate_4DTV_correction.pl $id >${id%%one-line}4dtv;done
（6）合并所有同源基因对的4dtv
for i in `ls *.4dtv`;do awk 'NR>1{print $1"\t"$3}' $i >>all-4dtv.txt;done
（7）合并所有同源基因对的kaks值
for i in `ls *.kaks`;do awk 'NR>1{print $1"\t"$3"\t"$4"\t"$5}' $i >>all-kaks.txt;done
（8）给结果文件进行排序和去冗余
sort all-4dtv.txt|uniq >all-4dtv.results
sort all-kaks.txt|uniq >all-kaks.results
（9）将kaks结果和4Dtv结果合并
join -1 1 -2 1 all-kaks.results all-4dtv.results >all-results.txt
（10）给结果文件添加标题
sed -i '1i\Seq\tKa\tKs\tKa/Ks\t4dtv_corrected' all-results.txt
四、可视化作图
1、数据准备
Arabidopsis_thaliana-4dtv.txt、Citrus_grandis-4dtv.txt、Citrus_sinensis-4dtv.txt、Malus_domestic-4dtv.txt
2、数据处理
cat citrus_sinensisl-4dtv.txt | awk '{print "C.sinensis""\t"$2}' > C.sinensis-4dtv.txt
Citrus_sinensis-4dtv.txt
C.sinensis-4dtv.txt
合并A.thaliana-4dtv.txt、C.grandis-4dtv.txt、C.sinensis-4dtv.txt、M.domestic-4dtv.txt文件
cat A.thaliana-4dtv.txt C.grandis-4dtv.txt C.sinensis-4dtv.txt M.domestic-4dtv.txt > 4species-4dtv.txt
### 给4species-4dtv.txt 文件添加标题
sed -i '1i\Species\tValue' 4species-4dtv.txt
3、RStudio作图
setwd("C:/Users/***/Desktop/4Dtv/")
### 读入4species-4dtv.txt文件
dtv <- read.table("4species-4dtv.txt", header = T, check.names =F, sep = "\t")
4species-4dtv.txt
### 载入R包
library(ggplot2)
library(hrbrthemes)
library(dplyr)
library(tidyr)
library(viridis)
### 绘图
p <- ggplot(data=dtv, aes(x=Value, group=Species, fill=Species)) + geom_density(adjust=1.5, alpha=.4) + theme_ipsum() + labs(title = "Distribution of 4Dtv distances", x = "4Dtv Value", y= "Density", size = 1.5)
Distribution of 4Dtv
这篇写了一晚上，总算写完了，看着结果还不错，图形画出来也不太丑，还好，哈哈。
深知自身水平有限，错误之处，欢迎指正。
以上，各位晚安。
补充一波：
1、关于PAML
其实计算Ka/Ks有很多种算法，KaKs_Calculator只是其中一种。目前在文章中见的较多的是用PAML中的codeml来计算，PAML可以利用DNA或Protein数据库使用最大似然法进行系统发育分析，也可以用于评估系统进化过程中的参数和假设检验，还可以构建系统进化树，但效果不太好。也有不少做全基因组复制分析的软件会调用PAML，其中wgd在做Ks分析时，用的就是PAML中的codeml来计算dN/dS。而且它还可对Ks分布的统计进行建模，例如核密度拟合(Kernel density estimate, KDE)或高斯混合模型(Gaussian mixture models)等。
PAML软件的详细用法请参考官方文档及陈连福的生信博客：http://abacus.gene.ucl.ac.uk/software/pamlDOC.pdf；http://www.chenlianfu.com/?p=2948；https://www.jianshu.com/p/3c26a5698f9c
2、关于KaKs_Calculator2.0运算模型选择
请参考官方文档，写的比较详细；KaKs_Calculator2.0_manual
3、关于wgd
用wgd做全基因组复制分析请参考：https://www.jianshu.com/p/3cfeb49c821a
参考：
https://www.jianshu.com/p/c7cbb6fed1d7
https://github.com/scbgfengchao/4DTv/blob/master/axt2one-line.py
http://chibba.pgml.uga.edu/mcscan2/#tm
Huang, S., Li, R., Zhang, Z. et al. The genome of the cucumber, Cucumis sativus L.. Nat Genet 41, 1275–1281 (2009). https://doi.org/10.1038/ng.475
Kumar, S., Stecher, G., Li, M., Knyaz, C. & Tamura, K. MEGA X: molecular evolutionary genetics analysis across computing platforms. Mol. Biol. Evol. 35,1547–1549 (2018).
Altschul, S.F., Gish, W., Miller, W., Myers, E.W. & Lipman, D.J. Basic local alignment search tool (1990) J. Mol. Biol. 215:403-410
Masami Hasegawa, Hirohisa Kishino and Taka-aki Yano. (1985) Dating of the Human-Ape Splitting by a Molecular Clock of Mitochondrial DNA. Journal of Molecular Evolution. 22:160-174
Wang, D., Zhang, Y., Zhang, Z., Zhu, J. & Yu, J. KaKs_Calculator 2.0: a toolkit incorporating gamma-series methods and sliding window strategies. Genomics, Proteom. Bioinforma. 8, 77–80 (2010).
Zhang, Z. et al. (2012) ParaAT: A parallel tool for constructing multiple protein-coding DNA alignments, Biochem Biophys Res Commun, 419, 779-781
