一合敦、4Dtv（四倍簡并位點(diǎn)顛換率）的定義

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值可以反映在進(jìn)化過程中的物種相對分化事件以及全基因組復(fù)制事件。

Sanwen huang et al., 2009

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

? ? ? ? Ka/Ks表示的是非同義替換（Ka）和同義替換（Ks）之間的比例充岛，這一比值可以判斷編碼該蛋白的基因是否遭受了選擇壓力。同義突變表示耕蝉，突變并不影響氨基酸序列崔梗，進(jìn)而不會(huì)影響蛋白結(jié)構(gòu)與功能。而非同義突變則會(huì)影響氨基酸序列垒在，可能會(huì)使其結(jié)構(gòu)和功能發(fā)生改變蒜魄，可能會(huì)遭受自然選擇。一般我們認(rèn)為场躯，同義突變不受自然選擇权悟，而非同義突變會(huì)遭受自然選擇作用。在生物進(jìn)化分析中推盛，知曉物種的同義突變和非同義突變發(fā)生的速率是非常有意義的峦阁。同義突變頻率即為Ks值，非同義突變頻率即為Ka值耘成，非同義突變率與同義突變率的比值即為Ka/Ks值榔昔。若Ka/Ks > 1，則認(rèn)為存在正選擇效應(yīng)（positive selection）瘪菌；若Ka/Ks = 1撒会，則認(rèn)為存在中性選擇效應(yīng)；如若Ka/Ks < 1师妙，則認(rèn)為存在負(fù)選擇效應(yīng)诵肛，即純化效應(yīng)或凈化選擇（purifying selection）。

三默穴、4Dtv及Ka/Ks值的計(jì)算

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)錄本序列，去冗余赃额。

同一個(gè)基因可能存在多個(gè)轉(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）運(yùn)行MCScanX(.blast & .gff文件 )

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

運(yùn)行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值計(jì)算

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

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

? ? ? ? 其中需要注意到一個(gè)問題轰胁，不少同學(xué)在安裝KaKs_Calculator2.0時(shí)會(huì)報(bào)錯(cuò)，當(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

? ? ? ? 這個(gè)腳本其實(shí)存在一個(gè)問題，簡書中也有其他作者指出來過意述，需要將腳本中密碼子表“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個(gè)線程

echo "12" >proc

# -m參數(shù)指定多序列比對程序?yàn)閏lustalw2，-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三個(gè)文件的基因ID需保持一致阔逼，且file.cds及file.pep為fasta格式，即“>”后面只接基因ID號地沮，否則會(huì)報(bào)錯(cuò)颜价，如下：

Citrus_sinensis.homolog

Citrus_sinensis.cds

Citrus_sinensis.pep

（4）使用KaKs_Calculator計(jì)算ka、ks值, -m參數(shù)指定kaks值的計(jì)算方法為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腳本計(jì)算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é)果文件進(jìn)行排序和去冗余

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é)果還不錯(cuò)，圖形畫出來也不太丑乾蛤，還好每界，哈哈。

? ? ? 深知自身水平有限家卖，錯(cuò)誤之處眨层，歡迎指正。

? ? ? ? 以上上荡，各位晚安趴樱。

補(bǔ)充一波：

1、關(guān)于PAML

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

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

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

請參考官方文檔坯临，寫的比較詳細(xì)；KaKs_Calculator2.0_manual

3恋昼、關(guān)于wgd

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

參考：

http://www.reibang.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

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