一、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.

? ? ? ? 其生物學(xué)意義為：共線性區(qū)域內(nèi)基因?qū)Φ?Dtv值可以反映在進化過程中的物種相對分化事件以及全基因組復(fù)制事件。

Sanwen huang et al., 2009

二、關(guān)于Ka/Ks

? ? ? ? Ka/Ks表示的是非同義替換（Ka）和同義替換（Ks）之間的比例，這一比值可以判斷編碼該蛋白的基因是否遭受了選擇壓力。同義突變表示，突變并不影響氨基酸序列，進而不會影響蛋白結(jié)構(gòu)與功能。而非同義突變則會影響氨基酸序列，可能會使其結(jié)構(gòu)和功能發(fā)生改變，可能會遭受自然選擇。一般我們認為，同義突變不受自然選擇，而非同義突變會遭受自然選擇作用。在生物進化分析中，知曉物種的同義突變和非同義突變發(fā)生的速率是非常有意義的。同義突變頻率即為Ks值，非同義突變頻率即為Ka值，非同義突變率與同義突變率的比值即為Ka/Ks值。若Ka/Ks > 1，則認為存在正選擇效應(yīng)（positive selection）；若Ka/Ks = 1，則認為存在中性選擇效應(yīng)；如若Ka/Ks < 1，則認為存在負選擇效應(yīng)，即純化效應(yīng)或凈化選擇（purifying selection）。

三、4Dtv及Ka/Ks值的計算

1、數(shù)據(jù)準(zhǔn)備

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、數(shù)據(jù)處理

取file.pep中最長轉(zhuǎn)錄本序列，去冗余。

同一個基因可能存在多個轉(zhuǎn)錄本

python3 extract_longest_iso.py -i Citrus_sinensis.fa -o L_Citrus_sinensis.fa

去冗余前44275條序列，去冗余后23394條序列；

去冗余后

python腳本如下

extract_longest_iso.py

3、獲取共線性基因?qū)?/h3>

（1）對蛋白序列構(gòu)建索引

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格式文件

僅保留染色體、轉(zhuǎn)錄本ID及坐標(biāo)信息

（4）運行MCScanX(.blast & .gff文件 )

? ? ? ? MCScanX是一款分析基因組內(nèi)或者基因組間的共線性區(qū)塊的軟件，它利用種內(nèi)或種間蛋白質(zhì)blastp比對結(jié)果，再結(jié)合編碼這些蛋白的基因在基因組中的位置坐標(biāo)，得到種內(nèi)或種間基因組的共線性區(qū)塊。MCScanX軟件安裝及詳細使用參見官網(wǎng)，安裝和使用都比較友好。http://chibba.pgml.uga.edu/mcscan2/#tm

運行MCScanX，獲取共線性基因?qū)?/p>

MCScanX Citrus_sinensis -g -3 -e 1e-10

? ? ? ? 得到 Citrus_sinensis.collinearity、Citrus_sinensis.tandem文件及Citrus_sinensis.html文件夾，其中我們需要的信息就在Citrus_sinensis.collinearity結(jié)果文件中。

Citrus_sinensis.collinearity結(jié)果文件

（5）提取共線性基因?qū)?/h4>

cat Citrus_sinensis.collinearity | grep "Cs" | awk '{print $3"\t"$4}' > Citrus_sinensis.homolog

Citrus_sinensis.homolog

4、Ka、Ks及4Dtv值計算

（1）準(zhǔn)備軟件及腳本

? ? ? ? 需要用到的軟件：KaKs_Calculator2.0及ParaAT；軟件安裝參考：http://cbb.big.ac.cn/Software

? ? ? ? 其中需要注意到一個問題，不少同學(xué)在安裝KaKs_Calculator2.0時會報錯，當(dāng)然我也碰到了，查了好久才解決，為了避免再次踩坑，特貼出解決方法；

###? make: *** [KaKs_Calculator] error.

解壓安裝包后，在“base.h”文件中最前面添加一行內(nèi)容：

#include <string.h>

然后保存、退出，再“make”命令安裝即可。

鏈接如下：https://www.researchgate.net/post/Can_someone_help_me_with_a_problem_about_KaKs_calculator20_used_in_linux_systems

? ? ? ? 需要用到的腳本：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混在一起了，是有問題的，在此統(tǒng)一用DNA codon。

calculate_4DTV_correction.pl

（2）數(shù)據(jù)準(zhǔn)備

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程序?qū)⒌鞍仔蛄斜葘Y(jié)果轉(zhuǎn)化為cds序列比對結(jié)果

# 新建proc文件, 設(shè)定12個線程

echo "12" >proc

# -m參數(shù)指定多序列比對程序為clustalw2，-p參數(shù)指定多線程文件，-f參數(shù)指定輸出文件格式為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參數(shù)指定kaks值的計算方法為YN模型

for i in `ls *.axt`;do KaKs_Calculator -i $i -o ${i}.kaks -m YN;done

# 將多行axt文件轉(zhuǎn)換成單行

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）合并所有同源基因?qū)Φ?dtv

for i in `ls *.4dtv`;do awk 'NR>1{print $1"\t"$3}' $i >>all-4dtv.txt;done

（7）合并所有同源基因?qū)Φ膋aks值

for i in `ls *.kaks`;do awk 'NR>1{print $1"\t"$3"\t"$4"\t"$5}' $i >>all-kaks.txt;done

（8）給結(jié)果文件進行排序和去冗余

sort all-4dtv.txt|uniq >all-4dtv.results

sort all-kaks.txt|uniq >all-kaks.results

（9）將kaks結(jié)果和4Dtv結(jié)果合并

join -1 1 -2 1 all-kaks.results all-4dtv.results >all-results.txt

（10）給結(jié)果文件添加標(biāo)題

sed -i '1i\Seq\tKa\tKs\tKa/Ks\t4dtv_corrected' all-results.txt

四、可視化作圖

1、數(shù)據(jù)準(zhǔn)備

Arabidopsis_thaliana-4dtv.txt、Citrus_grandis-4dtv.txt、Citrus_sinensis-4dtv.txt、Malus_domestic-4dtv.txt

2、數(shù)據(jù)處理

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 文件添加標(biāo)題

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

? ? ? ? 這篇寫了一晚上，總算寫完了，看著結(jié)果還不錯，圖形畫出來也不太丑，還好，哈哈。

? ? ? 深知自身水平有限，錯誤之處，歡迎指正。

? ? ? ? 以上，各位晚安。

補充一波：

1、關(guān)于PAML

? ? ? ? 其實計算Ka/Ks有很多種算法，KaKs_Calculator只是其中一種。目前在文章中見的較多的是用PAML中的codeml來計算，PAML可以利用DNA或Protein數(shù)據(jù)庫使用最大似然法進行系統(tǒng)發(fā)育分析，也可以用于評估系統(tǒng)進化過程中的參數(shù)和假設(shè)檢驗，還可以構(gòu)建系統(tǒng)進化樹，但效果不太好。也有不少做全基因組復(fù)制分析的軟件會調(diào)用PAML，其中wgd在做Ks分析時，用的就是PAML中的codeml來計算dN/dS。而且它還可對Ks分布的統(tǒng)計進行建模，例如核密度擬合(Kernel density estimate, KDE)或高斯混合模型(Gaussian mixture models)等。

PAML軟件的詳細用法請參考官方文檔及陳連福的生信博客：http://abacus.gene.ucl.ac.uk/software/pamlDOC.pdf；http://www.chenlianfu.com/?p=2948；http://m.itdecent.cn/p/3c26a5698f9c

2、關(guān)于KaKs_Calculator2.0運算模型選擇

請參考官方文檔，寫的比較詳細；KaKs_Calculator2.0_manual

3、關(guān)于wgd

用wgd做全基因組復(fù)制分析請參考：http://m.itdecent.cn/p/3cfeb49c821a

參考：

http://m.itdecent.cn/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

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