轉(zhuǎn)自：http://www.reibang.com/p/9d28de3d18e6

一、4Dtv（四倍簡(jiǎn)并位點(diǎn)顛換率）

1. 4倍簡(jiǎn)并位點(diǎn)（fourfold degenerate site）的定義

（1）如果密碼子的某個(gè)位點(diǎn)上任何核苷酸都編碼同樣的氨基酸，則稱這個(gè)位點(diǎn)為4倍簡(jiǎn)并位點(diǎn)程癌。
（2）例如甘氨酸密碼子（GGA, GGG, GGC, GGU）的第三個(gè)位點(diǎn)就是一個(gè)4倍簡(jiǎn)并位點(diǎn)御蒲，因?yàn)檫@個(gè)位點(diǎn)上所有的核苷酸替換（無(wú)論是A笑诅、G愕秫、U好芭、C）都是同義的承疲，即編碼同一個(gè)氨基酸邻耕。

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.
簡(jiǎn)單理解就是燕鸽，4倍簡(jiǎn)并位點(diǎn)第三個(gè)核酸密碼子的替換率兄世。

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的生物學(xué)意義

共線性區(qū)段所包含的基因?qū)Φ?DTV值可反映物種在進(jìn)化史中的物種相對(duì)分化事件以及全基因組復(fù)制事件。

image

image

其實(shí)啊研，我的理解是御滩，較多的基因?qū)?shù)存在4倍簡(jiǎn)并位點(diǎn)，說(shuō)明基因組多樣性較多党远，或者說(shuō)冗余基因較多削解，可能此刻發(fā)生了物種分化或者基因組復(fù)制。簡(jiǎn)單理解沟娱，這是一個(gè)變化巨大的過(guò)程钠绍。如果4倍簡(jiǎn)并位點(diǎn)替換率較低，或者說(shuō)達(dá)到一個(gè)平穩(wěn)狀態(tài)花沉，那么可能這個(gè)物種柳爽，并沒(méi)有重大事件發(fā)生媳握，一直平衡發(fā)展。

3. 關(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ù)處理，獲取最長(zhǎng)轉(zhuǎn)錄本

去冗余前44275條序列拔妥，去冗余后23394條序列忿危；

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

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

makeblastdb -in Citrus_sinensis.pep -dbtype prot -parse_seqids -out Citrus_sinensis

（2）blastp比對(duì)

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）運(yùn)行MCScanX(.blast & .gff文件 )

MCScanX是一款分析基因組內(nèi)或者基因組間的共線性區(qū)塊的軟件没龙，它利用種內(nèi)或種間蛋白質(zhì)blastp比對(duì)結(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é)果文件中肛鹏。

image

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

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

image

4逸邦、Ka、Ks及4Dtv值計(jì)算

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

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

其中需要注意到一個(gè)問(wèn)題缕减，不少同學(xué)在安裝KaKs_Calculator2.0時(shí)會(huì)報(bào)錯(cuò)，當(dāng)然我也碰到了芒珠，查了好久才解決桥狡，為了避免再次踩坑，特貼出解決方法皱卓；

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

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

#include <string.h>

然后保存、退出娜汁，再make命令安裝即可嫂易。

需要用到的腳本：calculate_4DTV_correction.pl；axt2one-line.py

來(lái)源：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è)問(wèn)題掐禁，簡(jiǎn)書中也有其他作者指出來(lái)過(guò)怜械，需要將腳本中密碼子表“my %codons=”中“U”改成“T”；這里RNA codon 和 DNA codon混在一起了傅事，是有問(wèn)題的缕允，在此統(tǒng)一用DNA codon。

image

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

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

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

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參數(shù)指定多序列比對(duì)程序?yàn)閏lustalw2，
-p參數(shù)指定多線程文件洒擦，
-f參數(shù)指定輸出文件格式為axt

此步需注意椿争，file.homolog、file.cds熟嫩、file.pep三個(gè)文件的基因ID需保持一致秦踪，且file.cds及file.pep為fasta格式，即“>”后面只接基因ID號(hào)掸茅，否則會(huì)報(bào)錯(cuò)椅邓，如下：

（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

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 文件添加標(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")

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

補(bǔ)充一波：

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

PAML軟件的詳細(xì)用法請(qǐng)參考官方文檔及陳連福的生信博客：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)算模型選擇
請(qǐng)參考官方文檔，寫的比較詳細(xì)试和；KaKs_Calculator2.0_manual

3讯泣、關(guān)于wgd
用wgd做全基因組復(fù)制分析請(qǐng)參考：http://www.reibang.com/p/3cfeb49c821a

4、4DTV計(jì)算上強(qiáng)調(diào)的是必須是4簡(jiǎn)并位點(diǎn)的同義替換, 對(duì)同義突變率的這種估計(jì)是一種相對(duì)保守的方法阅悍。而ds則不僅僅指4簡(jiǎn)并位點(diǎn)好渠，還包括2簡(jiǎn)并位點(diǎn)。
（1）4DTV具有較高的閾值节视，第三個(gè)密碼子是4簡(jiǎn)并位點(diǎn)的時(shí)候拳锚，才認(rèn)為是冗余基因中的密碼子，然后判斷WGD事件寻行。當(dāng)然霍掺，這些冗余基因經(jīng)過(guò)WGD后會(huì)消失或者變成假基因。
（2）ds則顯然寬松些拌蜘，2簡(jiǎn)并位點(diǎn)和4簡(jiǎn)并位點(diǎn)杆烁，都經(jīng)過(guò)計(jì)算，然后進(jìn)行冗余基因的判斷简卧，尋找峰值兔魂，判斷WGD事件。

參考文獻(xiàn)：
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

參考：
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