因為全基因組復(fù)制(Whole genome duplications, WGD)是生物進(jìn)化的重要因素之一, 所以WGD分析也是基因組分析經(jīng)常用到的一種分析方法。舉個例子，我們之所以能在多條染色體之間發(fā)現(xiàn)一些古老基因連鎖現(xiàn)象，是因為被子植物在過去2億年時間里就出現(xiàn)了多次的全基因組復(fù)制和基因丟失事件（見下圖，Tang et al., 2008）

古老WGD檢測有兩種方法，一種是共線性分析，另一種則是根據(jù)Ks分布圖。其中Ks定義為平均每個同義位點上的同義置換數(shù)，與其對應(yīng)的還有一個Ka，指的是平均每個非同義位點上的非同義置換數(shù)。

如果沒有WGD或是大片段重復(fù)，那么基因組中的旁系同源基因的同義置換符合指數(shù)分布(exponential distribution), 反之，Ks分布圖中就會出現(xiàn)一個由于WGD導(dǎo)致的正態(tài)分布峰(normal distributed peak). 而古老WGD的年齡則可通過分析這些峰中的同源置換數(shù)目來預(yù)測(Tiley et al., 2018)。

以發(fā)表在Science上的_Papaver somniferum_ L. 基因組文章中的圖Fig 1C為例，文章分別分析了_P. somniferum_ 和其他幾個物種的Ks分布。從Ks分布圖可以看到_P. somniferum_經(jīng)歷了一次比較近的WGD事件(Guo et al., 2018)。

文章中關(guān)于WGD的分析流程參考附錄8.1 Whole genome duplication in opium poppy genome, 我總結(jié)了關(guān)鍵的幾點

• 使用megablast進(jìn)行自比對，尋找基因組中共線性的區(qū)塊
• 使用BLASTP基于RBH( reciprocal best hits )進(jìn)行蛋白之間的相互比對
• 使用KaKs_Calculator基于YN模型計算RBH基因?qū)χ械腒s(synonymous substitution rate)
• 為了區(qū)分Ks中peak是WGD事件還是小規(guī)模重復(fù)，作者使用MCScanX進(jìn)行共線性分析，發(fā)現(xiàn)93.9%的RBH基因都是基因組內(nèi)共線性

目前WGD的分析流程也有人發(fā)了文章，我通過關(guān)鍵字"wgd pipeline"搜索找到了如下幾個:

• GenoDup: https://github.com/MaoYafei/GenoDup-Pipeline
• WGDdetector: https://github.com/yongzhiyang2012/WGDdetector
• wgd: https://github.com/arzwa/wgd

這一篇介紹的是wgd的用法。

軟件安裝

wgd目前無法用bioconda直接安裝，所以安裝起來稍顯麻煩，這是因為它的依賴軟件有點多。wgd依賴于以下軟件

• BLAST
• MCL
• MUSCLE/MAFFT/PRANK
• PAML
• PhyML/FastTree
• i-ADHoRe

但是好消息是它依賴的軟件大部分都可以用bioconda進(jìn)行安裝

conda create -n wgd python=3.5 blast mcl muscle mafft prank paml  fasttree cmake libpng mpi=1.0=mpich
conda activate wgd

這里選擇了mpi=1.0=mpich, 原因是i-adhore依賴于mpich. 如果安裝了openmpi就會出現(xiàn)error while loading shared libraries: libmpi_cxx.so.40: cannot open shared object file: No such file or directory

之后的安裝就簡單多了

git clone https://github.com/arzwa/wgd.git
cd wgd
pip install .
# 或者一行命令
pip install git+https://github.com/arzwa/wgd.git

對于i-ADHoRe，需要先在http://bioinformatics.psb.ugent.be/webtools/i-adhore/licensing/同意許可，才能下載i-ADHoRe-3.0

由于我的miniconda3安裝在~/opt/下，所以安裝路徑為~/opt/miniconda3/envs/wgd/

tar -zxvf i-adhore-3.0.01.tar.gz
cd i-adhore-3.0.01
mkdir -p build && cd build
cmake .. -DCMAKE_INSTALL_PREFIX=~/opt/miniconda3/envs/wgd/
make -j 4 
make insatall

軟件介紹

WGD一共有9個子模塊

• kde: 對Ks分布進(jìn)行KDE擬合
• ksd : Ks分布構(gòu)建
• mcl:All-vs-ALl的BLASP比對 + MCL分類分析.
• mix: Ks分布的混合建模.
• pre: 對CDS文件進(jìn)行預(yù)處理
• syn: 調(diào)用I-ADHoRe 3.0利用GFF文件進(jìn)行共線性分析
• viz: 繪制柱狀圖和密度圖
• wf1: 全基因組旁系同源基因組(paranome)的Ks標(biāo)準(zhǔn)分析流程，調(diào)用mcl, ksd和syn
• wf2: one-vs-one 同源基因(ortholog)的Ks標(biāo)準(zhǔn)分析流程，調(diào)用wcl和ksd

流程示意圖如下:

使用方法

以甘蔗的基因組 Saccharum spontaneum L 為例，基因組為8倍體，共32條染色體(2n = 4x8 = 32)

下載CDS和GFF注釋文件 tutorial

mkdir -p wgd_tutorial && cd wgd_tutorial
wget http://www.life.illinois.edu/ming/downloads/Spontaneum_genome/Sspon.v20190103.cds.fasta.gz
wget http://www.life.illinois.edu/ming/downloads/Spontaneum_genome/Sspon.v20190103.gff3.gz
gunzip *.gz

先用conda activate wgd啟動我們的分析環(huán)境，然后就開始分析了

第一步: 使用wgd mcl鑒定基因組內(nèi)的同源基因

wgd mcl -n 20 --cds --mcl -s Sspon.v20190103.cds.fasta -o Sspon_cds.out
# -n: 線程
# --cds: 輸入為cds
# --mcl: mcl聚類

這一步運行過程中，wgd會先檢查cds序列是否有效，也就是是否以ATG(起始密碼子)開頭，且以TAA/TAG/TGA(終止密碼子)結(jié)尾，然后將cds轉(zhuǎn)成氨基酸序列，之后用Blastp進(jìn)行相互比對，然后根據(jù)blastp結(jié)果用mcl聚類的方式尋找旁系同源基因。

輸出結(jié)果在Sspon_cds.out，有兩個結(jié)果輸出

• blast.tsv: BLASTP的outfmt6輸出結(jié)果
• blast.tsv.mcl: MCL聚類結(jié)果，每一行可以認(rèn)為是一個基因家族(gene family)

第二步: 使用wgd ksd構(gòu)建Ks分布

wgd ksd --n_threads 80 Sspon_cds.out/Sspon.v20190103.cds.fasta.blast.tsv.mcl Sspon.v20190103.cds.fasta

這一步也是先過濾cds中的無效數(shù)據(jù)，然后用mafft(默認(rèn))/muscle/prank對每個基因家族進(jìn)行多重序列聯(lián)配，用codeml計算dN/dS, 用alc/fasttree(默認(rèn))/phyml建樹

輸出結(jié)果在wgd_ksd目錄下

• ks.tsv: 每個基因家族中基因?qū)Φ母黜椊y(tǒng)計，其中包括Ka和Ks
• ks.svg: Ks分布，見下圖

第三步: 如果基因組質(zhì)量不錯，那么可以使用wgd syn進(jìn)行共線性分析。它能幫我們找到基因組內(nèi)的共線性區(qū)塊，以及相應(yīng)的錨定點

wgd syn --feature gene --gene_attribute ID -ks wgd_ksd/Sspon.v20190103.cds.fasta.ks.tsv Sspon.v20190103.gff3 Sspon_cds.out/Sspon.v20190103.cds.fasta.blast.tsv.mcl
#--feature: 從第三列選擇特征
#--gene_attribute: 從第九列提取編號

輸出圖片以.svg結(jié)尾，如下所示，圖中顏色紅藍(lán)代表Ks得分。

Ks的下游分析通常還包括對Ks分布的統(tǒng)計建模，這可以使用wgd kde進(jìn)行核密度擬合(Kernel density estimate, KDE)或用wgd mix建立高斯混合模型(Gaussian mixture models)

# KDE
wgd kde wgd_ksd/Sspon.v20190103.cds.fasta.ks.tsv
# Gaussian
wgd mix wgd_ksd/Sspon.v20190103.cds.fasta.ks.tsv

wgd kde輸出kde.svg, 而wdg mix則生成一個 wgd_mix文件夾。

混合模型常用于基于Ks分布研究WGD。基于一些基本的分子進(jìn)化假設(shè)，我們預(yù)期Ks分布中由WGD導(dǎo)致的peak應(yīng)該符合高斯分布，能夠用log-normal分布近似。但是考慮到混合模型容易過擬合，特別是Ks分布，因此作者不建議使用混合模型作為多次WGD假設(shè)的正式統(tǒng)計測試。

參考文獻(xiàn)

Tang, H., Bowers, J.E., Wang, X., Ming, R., Alam, M., and Paterson, A.H. (2008). Synteny and Collinearity in Plant Genomes. Science 320, 486–488.

Tiley, G.P., Barker, M.S., and Burleigh, J.G. (2018). Assessing the Performance of Ks Plots for Detecting Ancient Whole Genome Duplications. Genome Biol Evol 10, 2882–2898.

Guo, L., Winzer, T., Yang, X., Li, Y., Ning, Z., He, Z., Teodor, R., Lu, Y., Bowser, T.A., Graham, I.A., et al. (2018). The opium poppy genome and morphinan production. Science 362, 343–347.

Zwaenepoel, A., and Van de Peer, Y. (2019). wgd—simple command line tools for the analysis of ancient whole-genome duplications. Bioinformatics 35, 2153–2155.

總結(jié)

以上是生活随笔為你收集整理的python 密度聚类 使用_使用wgd进行全基因组复制分析的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。