Comparison of Bayesian and maximum-likelihood inference of population genetic parameters

Migrate-n 方法原始文獻(xiàn)，2006 年發(fā)表在 Bioinformatics 上。作者為 Peter Beerli。

使用 coalescence theory，結(jié)合了最大似然法（ML）和貝葉斯推斷，推斷群體大小、migration rates。

Parameters are scaled effective population sizes Θi?(4× effective population size × mutation rate per site per generation), and scaled immigration rates?Mji?(immigration rate divided by mutation rate).?

θ = 4Neμ，M = m/μ。m 是每一代遷徙的個(gè)體的比例，μ 是每個(gè)位點(diǎn)的突變率，Ne 是有效種群大小。

記一下我總結(jié)的 migrate-n 的用法：

migrate-n 是為桑格測(cè)序數(shù)據(jù)設(shè)計(jì)的一款進(jìn)行種群大小、遷移率推斷的軟件，可以首先設(shè)定幾個(gè)可能的模型（例如兩個(gè)種群，可以設(shè)置 full migration model 即兩個(gè)種群間有雙向基因流，也可以設(shè)置單向基因流的模型），運(yùn)行 migrate-n 得到每個(gè)模型的 marginal likelihood。

對(duì)于二代測(cè)序數(shù)據(jù)，migrate-n 有將 vcf 轉(zhuǎn)換為 migration 格式的腳本，見 tutorial。但是這個(gè)轉(zhuǎn)換格式的腳本不支持過多的染色體數(shù)目，例如我自己的數(shù)據(jù)是 scaffold level 的參考基因組，就無法完整地進(jìn)行轉(zhuǎn)換。而如果是染色體水平的參考基因組，雖然這一步可以成功轉(zhuǎn)換，但到了運(yùn)行的步驟也很可能失敗；例如我選擇了最大的 20 條 scaffolds，無法成功運(yùn)行，日志文件里并不明說理由，程序中止，很可能是因?yàn)閿?shù)據(jù)量太大。在我自己的例子里，減小數(shù)據(jù)量可以解決問題?？梢噪S機(jī)地選取一些不太大的 scaffolds。

設(shè)置模型：

我自己的二群體模型，簡(jiǎn)單的三種情況。

軟件的 tutorial 里詳細(xì)地介紹了模型的設(shè)置方法，可以非常復(fù)雜。除了migration model 以外還有 divergence model。我這里是 migration model。

運(yùn)行 migrate-n，在 menu 里設(shè)置好輸入文件、parameter 等，運(yùn)行。

我遇到了一個(gè)很奇怪的 bug，有幾個(gè)個(gè)體的序列就是讀不進(jìn)去。后來改了個(gè)體名稱就好了。

獲得的 outfile?？梢詫⒉煌Ｐ偷?outfile 分別命名為 model1.out 等。作者在 tutorial 里寫到用以下命令獲取最佳模型：

grep "All" model*.out | sort -n -k 4,4 | python bf.py

注意這個(gè) python 腳本不在軟件包里，需要用 tutorial 里的鏈接下載。此外在我自己的情況里，上述腳本不能直接用于計(jì)算，因?yàn)榘?All 的行不僅有模型總結(jié)參數(shù)，還有 θ 和 M 的估計(jì)，所以我自己又在 sort 之后加上了 head -n N，這個(gè) N 是比較的模型的數(shù)目。獲得的結(jié)果如下表：

軟件 tutorial 里給的示例

LBF 是貝葉斯因子的 log 值。等于 0 時(shí)意味著貝葉斯因子為 1，在其他地方讀到說這意味著完全否定原假設(shè) 。這里我想需要進(jìn)一步了解 migrate-n 的內(nèi)容才能讀明白了?？傊@里給出了最佳模型?？梢阅米罴涯Ｐ偷?outfile，選取出估計(jì)得到的 θ 和 M，進(jìn)一步得到有效種群大小和遷移率。

可以讀一讀其他用到這一軟件的文章。

Genetic diversity and structure of the Chinese lake gudgeon (?Sarcocheilichthys sinensis)

一種魚的線粒體 DNA 分析。用了許多模擬的方法，包括 ABC、S-DIVA、migrate-n、TMRCA，研究種群的地理歷史。證明了這個(gè)物種在長江、岷江間的擴(kuò)散歷史。

寫到用 migrate-n 估計(jì)了 mutation-scaled population size θ = 2Neμ，以及 mutation-scaled immigration rate M = m/μ。（對(duì)于線粒體的數(shù)據(jù)）

mitochondrial DNA?θ?=?2Neμ?(where Ne is the effective population size of females, and μ is the mutation rate).