用于beast并行化運行的beagle的安裝

此處的beagle和用于填充vcf文件缺失基因型的beagle不是同一個軟件。

# 下載源碼
git clone https://github.com/beagle-dev/beagle-lib.git
cd beagle-lib
mkdir build && cd build

# 配置編譯（安裝到路徑/share/home/XXX/soft/beast2/beagle）
# 清理之前的配置
cd ~/soft/beast2/beagle-lib/build
rm -rf *
##編譯使用cpu版本的beagle
cmake .. \
  -DCMAKE_INSTALL_PREFIX=/share/home/XXX/soft/beast2/beagle \
  -DBUILD_SHARED_LIBS=ON \
  -DBUILD_CUDA=OFF \
  -DBUILD_JNI=ON \
  -DBUILD_OPENCL=OFF \
  -DCMAKE_CXX_COMPILER=/share/home/XXX/soft/gcc/8.3.0/bin/g++ \
  -DCMAKE_C_COMPILER=/share/home/XXX/soft/gcc/8.3.0/bin/gcc 
# 編譯和安裝
make -j8
make install
# 把下面的一行內(nèi)容添加到~/.bashrc
export LD_LIBRARY_PATH=/share/home/XXX/soft/beast2/beagle/lib:$LD_LIBRARY_PATH
source ~/.bashrc

上面安裝的是cpu版本的beagle,我在此處指定了編譯的gcc和g++的絕對路徑。
如果你的服務(wù)器有GPU的話，可以使用下面的命令編譯支持顯卡的beagle

#編譯使用GPU版本的beagle 
cmake .. \
  -DCMAKE_INSTALL_PREFIX=/share/home/XXX/soft/beast2/beagle \
  -DBUILD_SHARED_LIBS=ON \
  -DBUILD_CUDA=ON \
  -DBUILD_JNI=ON \
  -DCUDA_TOOLKIT_ROOT_DIR=/share/home/XXX/CUDA11.8 \
  -DBUILD_OPENCL=OFF \
  -DCMAKE_CXX_COMPILER=/share/home/XXX/soft/gcc/8.3.0/bin/g++ \
  -DCMAKE_C_COMPILER=/share/home/XXX/soft/gcc/8.3.0/bin/gcc

編譯cuda版本的beagle報錯fatal error: libhmsbeagle/GPU/kernels/BeagleCUDA_kernels.h: No such file or directory，官方github說是gcc的版本不對，我此處修改為8.3.0版本可以成功編譯。
注意cmake指定參數(shù)后的所有路徑都要使用絕對路徑，不要使用相對路徑，例如~,這個會報錯。

準(zhǔn)備xml文件(在windows上使用BEAUti)

現(xiàn)在假設(shè)你已經(jīng)擁有比對好的多個物種的單拷貝序列的fa文件。
如果你的序列很長，例如超過10K，而且物種很多，則使用BEAUti會很慢。
所有我使用的魔法工具來完成。

截取前1K的序列

#提取序列的前1-1000個堿基用于測試數(shù)據(jù)
seqkit subseq -r 1:1000 All.min4.fasta|seqkit seq -w 0  >test.fa

test.fa導(dǎo)入BEAUti，進行設(shè)置，獲得xml文件。
打開 BEAUti2（通常位于 BEAST2 安裝目錄）。
導(dǎo)入序列文件：
File → Import Alignment → 選擇 test.fasta。
設(shè)置分區(qū)（Partitions）：
如果數(shù)據(jù)有多個基因/位點，可定義分區(qū)模型，此處選擇 nucleotide（DNA序列）或aminoacid （氨基酸序列）

選擇進化模型(Model site)：
Site Model → 選擇 Substitution Model（如 HKY、GTR）。
設(shè)置 Gamma 或 Invariant Sites（如 +G 或 +I）。
此處可以使用iqtree估計的最佳模型也可以使用beast的拓展包bModelTest做核酸模型的自動選擇。

設(shè)置分子鐘（Clock Model）：
Strict Clock（嚴格分子鐘）或 Relaxed Clock（松弛分子鐘，如 Log Normal）。實際使用Log Normal對數(shù)松弛分子鐘

設(shè)置樹先驗（Tree Prior）：
Coalescent（種群歷史）或 Yule（物種分化）。

設(shè)置 MCMC 參數(shù)：
一個粗略的計算方法是總代數(shù)=3000*(序列條數(shù)的平方)。
這里有100條序列，就是30,000,000代
Chain Length（鏈長，通常 10,000,000 次迭代）。

樣本數(shù)量=Chain Length /采樣頻率，要求樣本數(shù)量至少有10,000個。
所以此處采樣頻率最大為30,000,000/10,000=3000，此處實際設(shè)置為1000
Log Parameters（采樣頻率，如每 1000 代記錄一次）。
保存 .xml 文件：
File → Save → 生成 test1K.xml。

使用我自定義的腳本fa_to_xml.py
對于大的fa文件，可以使用我自己編寫的腳本fa_to_xml.py,把fa序列轉(zhuǎn)為需要的xml格式，使用這個腳本后可以節(jié)省幾天的時間。

python3 fa_to_xml.py All.min4.fasta test1K.xml test.xml

fa_to_xml.py腳本需要3個參數(shù)，

參數(shù)1是完整的輸入的fa文件
參數(shù)2是輸入上面使用BEAUti生成的xml的模板文件
參數(shù)3是輸出xml的文件名
經(jīng)過此處的修改后，輸出的test.xml文件即是下一步真正分析的xml文件。

使用beast2進行分析

如果你只安裝了cpu版本的beagle

beast -beagle -threads 24 test.xml

強制使用cpu的beagle

beast -beagle_CPU -threads 24 test.xml

強制使用GPU的beagle，放心增大使用的核心數(shù)，GPU的核心數(shù)比cpu多的多

beast -beagle_GPU -beagle_order 1,2 -threads 400 -overwrite  All.xml

beast常用的參數(shù)

-overwrite參數(shù)用于覆蓋之前斷掉的輸出文件
-threads 40 指定使用的cpu數(shù)量
-beagle/-beagle_CPU/-beagle_GPU 指定使用beagle進行加速，如果不指定，則只能使用1個cpu，會非常慢
-D chainLength=1000000 指定新的chainlength,替換xml中的值
-prefix 字符串指定輸出文件的前綴
-seed 9527 指定隨機數(shù)，不同批次的分析指定不同的隨機數(shù)種子和輸出前綴，這樣后面可以把多次運行的log和tree合并到一起。
注意：我這里使用了參數(shù)-beagle_order 1,2來指定使用編號為1和2的GPU顯卡。
beast -beagle_info命令會顯示，你當(dāng)前環(huán)境中可用的beagle的資源，我的是顯示如下：說明有2個GPU,所以可以指定使用order 1,2.

0 : CPU (x86_64)
    Flags: PRECISION_SINGLE PRECISION_DOUBLE COMPUTATION_SYNCH EIGEN_REAL EIGEN_COMPLEX SCALING_MANUAL SCALING_AUTO SCALING_ALWAYS SCALERS_RAW SCALERS_LOG VECTOR_SSE VECTOR_NONE THREADING_CPP THREADING_NONE PROCESSOR_CPU FRAMEWORK_CPU

1 : NVIDIA GeForce RTX 4090
    Global memory (MB): 24210
    Clock speed (Ghz): 2.52
    Number of cores: 16384
    Flags: PRECISION_SINGLE PRECISION_DOUBLE COMPUTATION_SYNCH COMPUTATION_ASYNCH EIGEN_REAL EIGEN_COMPLEX SCALING_MANUAL SCALING_AUTO SCALING_ALWAYS SCALERS_RAW SCALERS_LOG VECTOR_NONE THREADING_CPP THREADING_NONE PROCESSOR_GPU FRAMEWORK_CUDA PARALLELOPS_STREAMS PARALLELOPS_GRID

2 : NVIDIA GeForce RTX 4090
    Global memory (MB): 24210
    Clock speed (Ghz): 2.52
    Number of cores: 16384
    Flags: PRECISION_SINGLE PRECISION_DOUBLE COMPUTATION_SYNCH COMPUTATION_ASYNCH EIGEN_REAL EIGEN_COMPLEX SCALING_MANUAL SCALING_AUTO SCALING_ALWAYS SCALERS_RAW SCALERS_LOG VECTOR_NONE THREADING_CPP THREADING_NONE PROCESSOR_GPU FRAMEWORK_CUDA PARALLELOPS_STREAMS PARALLELOPS_GRID

如果你也遇到報錯CUDA error: "Out of memory" (2) from file </beagle-lib/libhmsbeagle/GPU/GPUInterfaceCUDA.cpp>, line 596. 說明你的xml的數(shù)據(jù)太大了，顯卡資源不足。

會輸出.log文件，

使用Tracer評估上面的MCMC的結(jié)果

使用Tracer github評估結(jié)果是否收斂。
這個是在windows端的軟件。

image.png

File > import 導(dǎo)入上面輸出的log文件即可
我此處的示例因為跑的mcmc次數(shù)太低，所以ESS的值很小，并不是可以使用的結(jié)果。ESS >200, 才表示采樣充分。需要增大chain length

合并后驗樹（TreeAnnotator）

打開 TreeAnnotator設(shè)置：
Burnin（丟棄前 10% 的樹）。
Posterior Probability Limit（如 0.5，僅保留高支持率節(jié)點）。
選擇 Maximum Clade Credibility Tree（最大分支可信樹）。
輸出 .tree 文件（可用于可視化）。
或者是linux命令行操作
treeannotator -burnin 30 -limit 0.5 test-test.trees test.trees

繪制得到的樹文件可視化，可以使用Figtree或beast2自帶的DensiTree或Y叔的ggtree包。

調(diào)參優(yōu)化ESS值

LogCombiner 是一個用于合并多個獨立運行的MCMC鏈（如重復(fù)運行或不同種子下的結(jié)果）或調(diào)整日志文件采樣頻率的工具，可有效提高ESS值（Effective Sample Size）。
如果你的ESS值始終不理想，可以使用Beast2自帶的一個程序LogCombiner來優(yōu)化log文件和tree文件。
注意：前提條件是你生成所有的log的xml參數(shù)除了隨機數(shù)不同之外，其他的都必須一樣。否則無法合并。

合并多個log文件

logcombiner -log run1.log -log run2.log -log run3.log -o combined.log -b 10 -resample 5000

參數(shù)說明：

-b 丟棄的burn-in比例默認是10，即10%。
-log 指定輸入log文件
-resample 重新采樣頻率（步數(shù)）合并后如果文件太大了，可以設(shè)置稍大的采樣數(shù)，比如從1000調(diào)整到5000
-o 輸出log文件名

合并tree文件

logcombiner -log run1.trees -log run2.trees -log run3.trees  -o combined.trees -b 10 -resample 5000

合并后再使用Tracer查看log文件，看ESS值是否大于200。

beast的插件安裝

beast支持win,linux系統(tǒng)。
BEAUti是其中一個程序，如果linux沒有X11，沒有的話沒法用窗口模式操作。這一步可以在windows上進行操作。

beast的windows和linux的插件是相同的。
例如：安裝使用用于SNP數(shù)據(jù)分析SNAPP插件
windows端點擊左側(cè)的File > manage package 就可以看到下面的界面，點擊下面的?就可以出現(xiàn)你的默認的包的安裝位置了。

image.png

上面顯示我的包的位置是C:\Users\XXX\BEAST\2.7

windows安裝SNAPP包

windows你可以直接選擇對應(yīng)的包，然后點擊下面的Install/Update按鈕即可安裝。如果安裝失敗，可能由于網(wǎng)絡(luò)問題，則可以直接去下載對應(yīng)的zip壓縮包，解壓縮后放入上面的包的路徑中。
linux的位置是~/.beast/2.7/
如果你有多個版本的beast,比如2.7.6和2.7.7，實際使用時會優(yōu)先使用最高版本的beast.

linux安裝SNAPP包

packagemanager -add SNAPP

這里的packagemanager是beast2自帶的腳本，如果上面還是報錯，會提示無法訪問網(wǎng)絡(luò)信息。
進入下面的路徑https://raw.githubusercontent.com/CompEvol/CBAN/master/packages2.7.xml ，去里面找到對應(yīng)的beast2版本的SNAPP的下載路徑。

mkdir ~/.beast/2.7/SNAPP
cd ~/.beast/2.7/SNAPP
wget -c https://github.com/BEAST2-Dev/SNAPP/releases/download/v1.6.0/SNAPP.addon.v1.6.0.zip
unzip SNAPP.addon.v1.6.0.zip
rm -rf SNAPP.addon.v1.6.0.zip

windows同樣把上面的zip解壓到文件夾SNAPP,然后把這個文件夾放到C:\Users\XXX\BEAST\2.7路徑下。

下面以使用命令行的SNAPP分析SNP數(shù)據(jù)為例

##刪除原有的beauti.properties
rm -rf ~/.beast/2.7/beauti.properties
##再次運行packagemanager查看安裝情況
packagemanager -list

#從github下載ruby的腳本
wget https://raw.githubusercontent.com/mmatschiner/snapp_prep/master/snapp_prep.rb
#生成xml文件
ruby snapp_prep.rb -a SNAPPER -v NC_031969.f5.sub5.vcf -t individuals.txt -c constraints.txt -m 1000 -l 100000

snapp_prep.rb的命令參數(shù)

-a 指定分析時使用的是SNAPP或SNAPPER,默認是SNAPP
-v vcf文件
-t txt文本文件，第1列是物種信息，第2列是樣本信息，tab分割符號
-c 化石證據(jù)年齡信息文件
-l MCMC的代數(shù)，默認是500000
-m 最多使用多少SNP.無默認值
-x 輸出的xml文件的前綴
-o 輸出的tree和log文件的前綴

individuals.txt的文件格式如下（所有的樣本名都需要在這個文件里，tab分割）

species individual
物種名稱1   樣本名1
物種名稱2   樣本名3
物種名稱1   樣本名2

constraints.txt （指定物種的分化時間，中間的19數(shù)字代表19MYA）

  lognormal(0,19,0.16)  crown   物種名稱1,物種名稱2,物種3

Tracer的界面詳細解析

image.png

使用Tracer導(dǎo)入log文件之后的界面如右側(cè)。

點擊屏幕右側(cè)上面的Trace
可以看到入上面的分布圖，橫坐標(biāo)是訓(xùn)練的chain length,可以看出隨著訓(xùn)練的次數(shù)的增加，左側(cè)對應(yīng)參數(shù)的變化。左側(cè)有對應(yīng)參數(shù)的ESS值，ESS小于100是紅色，100<ESS<200 是黃色，ESS >200是正常的黑色。黑色的ESS值才表示模型可信。如果ESS較小，算出的95% HPD的區(qū)間會非常大。
2 .點擊屏幕右側(cè)上面的Estimates

image.png

此時能看到分布的柱狀圖和上面的表格文件
能夠看到95% HPD interval和ESS,均值和中間值等。
STDEV標(biāo)準(zhǔn)差，
number of samples 樣本數(shù)，這個需要達到至少10000才可信。
如果你有多次運行的log文件，也可以點擊左側(cè)的+導(dǎo)入進來

image.png

現(xiàn)在我們就可以看到combine后的分布，可以看到我的這個是咧，ESS值才51，右側(cè)的Trace分布圖，后面還出現(xiàn)了斷層，說明還需要加大訓(xùn)練次數(shù)，即增加chain length。
一般要求的ESS>200,才被接受，如果某些雜志要求嚴格，可能要ESS>400.

增大ESS的方法

增加采樣頻率()
增加chain的長度（chain length）
多個log文件合并分析
如果你反復(fù)運行同一個模型，比如超過10億次，還不收斂（ESS<200）,此時建議你換一個模型試試。不同的模型的log文件無法合并。

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av

物種分化時間估計工具beast2

物種分化時間估計工具beast2

用于beast并行化運行的beagle的安裝

準(zhǔn)備xml文件(在windows上使用BEAUti)

fa_to_xml.py腳本需要3個參數(shù)，

使用beast2進行分析

使用Tracer評估上面的MCMC的結(jié)果

合并后驗樹（TreeAnnotator）

繪制得到的樹文件可視化，可以使用Figtree或beast2自帶的DensiTree或Y叔的ggtree包。

調(diào)參優(yōu)化ESS值

合并多個log文件

合并tree文件

beast的插件安裝

windows安裝SNAPP包

linux安裝SNAPP包

下面以使用命令行的SNAPP分析SNP數(shù)據(jù)為例

individuals.txt的文件格式如下（所有的樣本名都需要在這個文件里，tab分割）

constraints.txt （指定物種的分化時間，中間的19數(shù)字代表19MYA）

Tracer的界面詳細解析

增大ESS的方法

相關(guān)閱讀更多精彩內(nèi)容

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

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九 欧美,1769亚洲,黄色成人av

物種分化時間估計工具beast2

用于beast并行化運行的beagle的安裝

準(zhǔn)備xml文件(在windows上使用BEAUti)

fa_to_xml.py腳本需要3個參數(shù)，

使用beast2進行分析

使用Tracer評估上面的MCMC的結(jié)果

合并后驗樹（TreeAnnotator）

繪制得到的樹文件可視化，可以使用Figtree或beast2自帶的DensiTree或Y叔的ggtree包。

調(diào)參優(yōu)化ESS值

合并多個log文件

合并tree文件

beast的插件安裝

windows安裝SNAPP包

linux安裝SNAPP包

下面以使用命令行的SNAPP分析SNP數(shù)據(jù)為例

individuals.txt的文件格式如下（所有的樣本名都需要在這個文件里，tab分割）

constraints.txt （指定物種的分化時間，中間的19數(shù)字代表19MYA）

Tracer的界面詳細解析

增大ESS的方法

相關(guān)閱讀更多精彩內(nèi)容

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

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av

繪制得到的樹文件可視化，可以使用Figtree或beast2自帶的DensiTree或Y叔的ggtree包。

individuals.txt的文件格式如下（所有的樣本名都需要在這個文件里，tab分割）

constraints.txt （指定物種的分化時間，中間的19數(shù)字代表19MYA）