R語言怎么構(gòu)建進(jìn)化樹

這篇文章主要介紹“R語言怎么構(gòu)建進(jìn)化樹”的相關(guān)知識，小編通過實(shí)際案例向大家展示操作過程，操作方法簡單快捷，實(shí)用性強(qiáng)，希望這篇“R語言怎么構(gòu)建進(jìn)化樹”文章能幫助大家解決問題。

樹的軟件構(gòu)建進(jìn)化樹的方法和軟件很多，前面我們講解構(gòu)建進(jìn)化樹的原理時提過，最準(zhǔn)確的方法為貝葉斯方法，但是貝葉斯的方法計(jì)算量太大太耗時，對于大量的數(shù)據(jù)不適用，其次就是最大似然法，這里我們解釋三種利用最大似然法構(gòu)建進(jìn)化樹的軟件：fasttree，iqtree2，RAxML-ng.輸入的數(shù)據(jù)為做群體遺傳進(jìn)化得到的SNP數(shù)據(jù)。

數(shù)據(jù)準(zhǔn)備做完重測序分析后，我們可以得到包含 SNP的變異結(jié)果vcf文件作為輸入文件。詳見課程：重測序數(shù)據(jù)分析課程，然后利用億速云/pop-evol-gwas:v1.3 鏡像進(jìn)行前期數(shù)據(jù)準(zhǔn)備與后續(xù)進(jìn)化樹構(gòu)建分析：

#運(yùn)行環(huán)境準(zhǔn)備：docker鏡像啟動
#鏡像下載：
docker pull 億速云/pop-evol-gwas:v1.3
#啟動遺傳進(jìn)化鏡像
docker run --rm -it -m 4G --cpus 1  -v D:\pop:/work 億速云/pop-evol-gwas:v1.3

#對vcf文件進(jìn)行數(shù)據(jù)過濾
vcftools --gzvcf all.varFilter.vcf.gz --recode --recode-INFO-all --stdout     --maf 0.05  --max-missing 0.7  --minDP 4  --maxDP 1000      \
    --minQ 30 --minGQ 0 --min-alleles 2  --max-alleles 2 --remove-indels |gzip - > clean.vcf.gz

#利用tassel軟件對文件進(jìn)行排序
run_pipeline.pl -Xmx30G  -SortGenotypeFilePlugin -inputFile clean.vcf.gz \
    -outputFile clean.sorted.vcf.gz -fileType VCF

#vcf文件格式轉(zhuǎn)換成Phylip格式，用于后續(xù)構(gòu)建進(jìn)化樹
run_pipeline.pl  -Xmx5G -importGuess  $workdir/00.filter/clean.sorted.vcf.gz  \
    -ExportPlugin -saveAs supergene.phy -format Phylip_Inter1. FastTree 構(gòu)建進(jìn)化樹

FastTree 是基于最大似然法構(gòu)建進(jìn)化樹的軟件，它最大的特點(diǎn)就是運(yùn)行速度快，支持幾百萬條序列的建樹任務(wù)。但是fasttree不支持bootstrap檢驗(yàn)以及支持的替換模型有限。官網(wǎng)如下：http://www.microbesonline.org/fasttree/

替換模型選擇

FastTree 支持核酸和蛋白的進(jìn)化樹構(gòu)建，對于核酸，可選的替換模型包括以下幾種：JC（Jukes-Cantor）、GTR（generalized time-reversible），默認(rèn)的模型為JC。對于蛋白質(zhì)，可選的替換模型包括以下幾種：JTT (Jones-Taylor-Thornton 1992)、LG(Le and Gascuel 2008)、WAG(Whelan & Goldman 2001) 默認(rèn)的模型為JTT。FastTree要求輸入的多序列比對結(jié)果為FASTA或者Phylip格式。
構(gòu)建進(jìn)化樹命令舉例：

fasttree -nt -gtr  supergene.fa   >  fasttree.nwk

FastTree 其他命令舉例：

#對于蛋白質(zhì)的進(jìn)化樹構(gòu)建，基本用法如下
fasttree protein.fasta > tree
#也可以選擇LG或者WAG替換模型，用法如下
fasttree -lg protein.fasta > tree
fasttree -wag protein.fasta > tree
#對于核酸序列，基本用法如下
fasttree -nt nucleotide.fasta > tree
#也可以選擇GTR替換模型，用法如下
fasttree -nt  -gtr nucleotide.fasta > tree

. IQ-tree 構(gòu)建進(jìn)化樹IQ-tree也是一款最大似然法構(gòu)建進(jìn)化樹的軟件，目前IQ-tree已經(jīng)更新到2.0版本功能和性能也有很大的提升，主要有四大功能，高效建樹（efficient tree reconstruction），模型選擇（modelfinder: fast and accurate model selection），超快bootstrap（ultrafast bootstrap approximation），大型數(shù)據(jù)（big data analysis），以上特點(diǎn)特別適用于高通量測序的大量SNP構(gòu)建進(jìn)化樹。

模型選擇

構(gòu)建進(jìn)化樹有很多模型對于初學(xué)者往往不知道那種模型最適合，iqtree提供自動模型選擇功能，使用的軟件為modelfinder。Modelfinder是一款速度超快的自動最佳模型選擇軟件。其在保證準(zhǔn)確性的情況下，速度上比jmodeltest（針對DNA）和prottest（針對蛋白）要快100倍（ModelFinder is up to 100 times faster than jModelTest/ProtTest.），使用命令舉例：

#自動選擇最佳模型并構(gòu)建進(jìn)化樹：-m MFP
iqtree -s supergene.phy -m MFP
#只是想找最佳模型不構(gòu)建進(jìn)化樹：

iqtree -s example.phy -m MF
#查找模型計(jì)算過程：


ModelFinder will test up to 546 protein models (sample size: 36415) ...
 No. Model         -LnL         df  AIC          AICc         BIC
  1  LG            10134094.366 349 20268886.731 20268893.505 20271854.186
  2  LG+I          10133927.677 350 20268555.354 20268562.167 20271531.312
  3  LG+G4         10043239.052 350 20087178.104 20087184.917 20090154.062
  4  LG+I+G4       10043175.024 351 20087052.048 20087058.900 20090036.508
  5  LG+R2         10063911.721 351 20128525.442 20128532.294 20131509.902
  6  LG+R3         10045448.117 353 20091602.235 20091609.165 20094603.701

MFP為ModelFinder Plus的縮寫，該參數(shù)使程序執(zhí)行ModelFinder選擇最適模型并完成建樹分析。ModelFinder為許多不同的模型計(jì)算初始簡約樹的邏輯概率，并產(chǎn)生Akaike information criterion (AIC),* corrected Akaike information criterion* (AICc), and* the Bayesian information criterion* (BIC)三個結(jié)果標(biāo)準(zhǔn)值，通常ModelFinder選擇BIC分?jǐn)?shù)最低的模型（當(dāng)然也可以指定AIC和AICc通過指定選項(xiàng)-AIC或者-AICc)。如果你想節(jié)約時間，可指定選擇的模型及編碼參數(shù)，例如:從WAG ,LG ,JTT 核苷酸替代模型里選一個: -mset WAG,LG,JTT；在+G和+I，以及+I+G三個里面選擇rate ：-mrate G,I,I+G；heterogeneity參數(shù)：-mfreq FU,F命令行如下：

iqtree -s example.phy -m MPF -mset WAG,LG,JTT  -mrate G,I,I+G -mfreq FU,F

指定模型參數(shù)設(shè)置格式：-m MODEL+FreqType+RateType，

MODEL：model name
+FreqType：(可選項(xiàng))frequency type
+RateType：(可選項(xiàng))rate heterogeneity type

替換模型MODEL包括：

• DNA模型:

JC/JC69, F81, K2P/K80, HKY/HKY85, TN/TrN/TN93, TNe,
K3P/K81, K81u, TPM2, TPM2u, TPM3, TPM3u, TIM, TIMe,
TIM2, TIM2e, TIM3, TIM3e, TVM, TVMe, SYM, GTR and 6-digit

• 蛋白質(zhì)模型：

BLOSUM62, cpREV, Dayhoff, DCMut, FLU, HIVb, HIVw, JTT,
JTTDCMut, LG, mtART, mtMAM, mtREV, mtZOA, mtMet, mtVer,
mtInv, Poisson, PMB, rtREV, VT, WAG

+FreqType堿基使用偏好：Base frequencies 可選設(shè)置：

每個核苷酸位點(diǎn)上的替代是隨機(jī)發(fā)生的，則A,T,C,G出現(xiàn)的頻率應(yīng)該大致相等。實(shí)際情況：DNA受到自然選擇的壓力，各個位點(diǎn)的堿基出現(xiàn)頻率并不相等。

+RateType：rate heterogeneity across sites可選設(shè)置：

指定替換模型構(gòu)建進(jìn)化樹命令舉例：

iqtree -s example.phy -m TIM2+I+G

Bootstrap法評估分支支持度

真實(shí)的進(jìn)化信息只有一個，而我們總是拿著有限的序列信息，希望去獲得他。能否獲得他，是一個問題。而我們使用的序列信息是否能真實(shí)且穩(wěn)定地反應(yīng)一個進(jìn)化信息，那么是另外一個事情。bootstrap法常用的，尤其是ML法構(gòu)建進(jìn)化樹上，分支可靠性檢驗(yàn)方法。但是這個計(jì)算邏輯最大的問題在于，抽樣重新跑，抽樣再重新跑，不斷重復(fù)，直到收斂或者是到指定的比如1000次。計(jì)算量大，耗時長。IQ-tree的作者團(tuán)隊(duì)在前述提出了一個快速的BS方法，最后整合到IQ-tree中。使用的方式是

iqtree -s example.phy -m TIM2+I+G -b 1000

超快（ultrafast bootstraping）

大概是IQTREE的精髓所在。顧名思義，ultrafast bootstrap approximation的特點(diǎn)就是超快。這里涉及到的細(xì)節(jié)，感興趣的讀者可以參見IQTREE的開發(fā)者寫的幾篇文章。作者認(rèn)為，UFBoot is 10 to 40 times faster than RAxML rapid bootstrap and obtains less biased support values。

iqtree -s example.phy -m TIM2+I+G -B 1000

除ultrafast bootstrap之外，IQTREE還提供了以下檢驗(yàn)樹的拓?fù)浣Y(jié)構(gòu)可信度的方法。

• -alrt：SH-aLRT檢驗(yàn)(4)，沒記錯的話FastTree2使用的就是這個吧？

• -abayes：approximate Bayes test，由瑞士蘇黎世應(yīng)用科學(xué)教授Maria Anisimova提出(5)

• -lbp：fast local bootstrap probability method，由Adachi and Hasegawa提出

iqtree -s example.phy -m TIM2+I+G -B 1000 -alrt 1000

如果你指定了多個檢驗(yàn)方法，那么其結(jié)果會在樹里（.treefile）呈現(xiàn)出來，不同檢驗(yàn)所得數(shù)值用斜線隔開，比如：((a,b)100/98:0.1,c:0.2)90/95最后，iqtree做群體遺傳進(jìn)化構(gòu)建進(jìn)化樹推薦命令：

iqtree2 -s supergene.phy -st DNA -T 2  -mem 8G \
    -m  GTR  -redo \
    -B 1000 -bnni \
    --prefix iqtree

3. RAxML 構(gòu)建進(jìn)化樹RAxML是最大似然法（maximum likelihood）建樹的經(jīng)典工具，其由來自德國海德堡理論科學(xué)研究所（Heidelberg Institute for Theoretical Studies）的Alexandros Stamatakis開發(fā),最新已經(jīng)更新了版本RAxML-NG，支持的替換模型更多，運(yùn)行速度更快。

RAxML 建樹原理

RAxML采用最大似然法建樹，即將系統(tǒng)樹的拓?fù)浣Y(jié)構(gòu)、分支長度及進(jìn)化模型等的全部或者部分作為需要估計(jì)的參數(shù)，在給定的數(shù)據(jù)集及進(jìn)化模型的基礎(chǔ)上，用最大似然法的標(biāo)準(zhǔn)-似然值最大化來估計(jì)這些參數(shù)。首先，要選擇進(jìn)化模型，以簡約樹或者聯(lián)接樹為基礎(chǔ)，采用似然法估計(jì)模型中各個參數(shù)。設(shè)置好參數(shù)后，以簡約樹或者聯(lián)接樹作為起始樹，進(jìn)行似然分析，最后用統(tǒng)計(jì)學(xué)方法從多個似然樹中尋找最佳得分樹。

RAxML-NG 使用

RAxML軟件支持輸入文件的格式可以是比對好的fasta格式或者phylip格式，例如DNA比對序列，核苷酸替代模型為GTR，rate heterogeneity設(shè)置為gamma分布，不做bootstraping，命令如下：

raxml-ng  --msa supergene.phy --model GTR+G  --prefix raxml_tree     --threads 2 --seed 123

如果是建樹和bootstrap一起做，可以加--all參數(shù)一步完成：

raxml-ng  --msa supergene.phy --model GTR+G  --prefix raxml_tree --threads 2 --seed 123

關(guān)于“R語言怎么構(gòu)建進(jìn)化樹”的內(nèi)容就介紹到這里了，感謝大家的閱讀。如果想了解更多行業(yè)相關(guān)的知識，可以關(guān)注億速云行業(yè)資訊頻道，小編每天都會為大家更新不同的知識點(diǎn)。