全基因組復(fù)制 WGD |（二）Ka/Ks及4Dtv值計(jì)算

轉(zhuǎn)自：http://www.itdecent.cn/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)上所有的核苷酸替換（無論是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)，說明基因組多樣性較多，或者說冗余基因較多，可能此刻發(fā)生了物種分化或者基因組復(fù)制。簡(jiǎn)單理解，這是一個(gè)變化巨大的過程。如果4倍簡(jiǎn)并位點(diǎn)替換率較低，或者說達(dá)到一個(gè)平穩(wěn)狀態(tài)，那么可能這個(gè)物種，并沒有重大事件發(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ù)處理，獲取最長轉(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è)問題，不少同學(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

來源：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è)問題，簡(jiǎn)書中也有其他作者指出來過，需要將腳本中密碼子表“my %codons=”中“U”改成“T”；這里RNA codon 和 DNA codon混在一起了，是有問題的，在此統(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來計(jì)算，PAML可以利用DNA或Protein數(shù)據(jù)庫使用最大似然法進(jìn)行系統(tǒng)發(fā)育分析，也可以用于評(píng)估系統(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。而且它還可對(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.itdecent.cn/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.itdecent.cn/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)過WGD后會(huì)消失或者變成假基因。
（2）ds則顯然寬松些，2簡(jiǎn)并位點(diǎn)和4簡(jiǎn)并位點(diǎn)，都經(jīng)過計(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.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

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

友情鏈接更多精彩內(nèi)容