全基因组复制 WGD |（二）Ka/Ks及4Dtv值计算
转自：https://www.jianshu.com/p/9d28de3d18e6
一、4Dtv（四倍简并位点颠换率）
1. 4倍简并位点（fourfold degenerate site）的定义
（1）如果密码子的某个位点上任何核苷酸都编码同样的氨基酸，则称这个位点为4倍简并位点。
（2）例如甘氨酸密码子（GGA, GGG, GGC, GGU）的第三个位点就是一个4倍简并位点，因为这个位点上所有的核苷酸替换（无论是A、G、U、C）都是同义的，即编码同一个氨基酸。
4倍是不是A, T, C, G的意思 ？
The 4DTv is calculated as the number of transversions at all four fold degenerate third codon positions divided by the number of fourfold degenerate third codon positions.
简单理解就是，4倍简并位点第三个核酸密码子的替换率。
image.png
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.
2. 4DTV的生物学意义
共线性区段所包含的基因对的4DTV值可反映物种在进化史中的物种相对分化事件以及全基因组复制事件。
image
image
其实，我的理解是，较多的基因对数存在4倍简并位点，说明基因组多样性较多，或者说冗余基因较多，可能此刻发生了物种分化或者基因组复制。简单理解，这是一个变化巨大的过程。如果4倍简并位点替换率较低，或者说达到一个平稳状态，那么可能这个物种，并没有重大事件发生，一直平衡发展。
3. 关于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、数据处理，获取最长转录本
去冗余前44275条序列，去冗余后23394条序列；
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
image
image
（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结果文件中。
image
（5）提取共线性基因对
cat Citrus_sinensis.collinearity | grep "Cs" | awk '{print $3"\t"$4}' > Citrus_sinensis.homolog
image
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/JinfengChen/Scripts/blob/master/FFgenome/03.evolution/distance_kaks_4dtv/bin/calculate_4DTV_correction.pl
https://github.com/scbgfengchao/4DTv/blob/master/axt2one-line.py
这个脚本其实存在一个问题，简书中也有其他作者指出来过，需要将脚本中密码子表“my %codons=”中“U”改成“T”；这里RNA codon 和 DNA codon混在一起了，是有问题的，在此统一用DNA codon。
image
（2）数据准备
Citrus_sinensis.homolog 、Citrus_sinensis.cd、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
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 &
-m参数指定多序列比对程序为clustalw2，
-p参数指定多线程文件，
-f参数指定输出文件格式为axt
此步需注意，file.homolog、file.cds、file.pep三个文件的基因ID需保持一致，且file.cds及file.pep为fasta格式，即“>”后面只接基因ID号，否则会报错，如下：
（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
image
image
合并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")
image
> ### 载入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)
image
补充一波：
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
4、4DTV计算上强调的是必须是4简并位点的同义替换, 对同义突变率的这种估计是一种相对保守的方法。而ds则不仅仅指4简并位点，还包括2简并位点。
（1）4DTV具有较高的阈值，第三个密码子是4简并位点的时候，才认为是冗余基因中的密码子，然后判断WGD事件。当然，这些冗余基因经过WGD后会消失或者变成假基因。
（2）ds则显然宽松些，2简并位点和4简并位点，都经过计算，然后进行冗余基因的判断，寻找峰值，判断WGD事件。
参考文献：
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
参考：
https://www.jianshu.com/p/c7cbb6fed1d7
https://github.com/JinfengChen/Scripts/blob/master/FFgenome/03.evolution/distance_kaks_4dtv/bin/calculate_4DTV_correction.pl
https://github.com/scbgfengchao/4DTv/blob/master/axt2one-line.py
http://chibba.pgml.uga.edu/mcscan2/#tm
http://cbb.big.ac.cn/Software
