1. 群體演化歷史分析的案例：
群體演化包括有效群體大小(Ne)變化、基因流，遷徙，分化等，會(huì)對(duì)等位基因頻率產(chǎn)生顯著影響，塑造了現(xiàn)有遺傳多樣性的模式和水平。群體演化、遺傳漂變和自然選擇共同決定了基因組遺傳多樣性的命運(yùn)。
Example: 蘇格蘭罕見(jiàn)的食肉褐鱒（rare piscivorous brown trout: ferox）與普通褐鱒（normal brown trout）在進(jìn)化過(guò)程中是否有基因交流？
2. 如何通過(guò)測(cè)序數(shù)據(jù)進(jìn)行群體演化歷史的推斷？
（1）群體的系譜進(jìn)化樹(shù)的形狀包含了群體演化歷史的模式
（2）在不同的群體演化歷史的模式下，SFS有不同的分布
什么是SFS(site frequency spectrum)：
請(qǐng)參考位點(diǎn)頻譜(SFS)計(jì)算
單個(gè)群體
兩個(gè)群體
顏色越偏紅，表示數(shù)量越多，越偏藍(lán)表示數(shù)量越少。如果兩個(gè)群體完全分開(kāi)，那它們derived allele頻率相同的交集就越少，表現(xiàn)在2D SFS上就是密度偏向各自的坐標(biāo)軸，如果群體交流混合，它們derived allele頻率相同的交集就越多，表現(xiàn)在2D SFS上就是密度偏向x=y的這條對(duì)角線。后面兩個(gè)模型對(duì)SFS的影響很像，都是使兩群體的SFS趨同，可能結(jié)合群體分化時(shí)間，核苷酸突變速率等推斷具體是哪種模型。
（3）通過(guò)最大似然估計(jì)得到最佳參數(shù)
個(gè)人理解是，通過(guò)推斷參數(shù)，給出模型，使expected SFS更合理，也就是更加類(lèi)似于兩個(gè)分離的群體。模型參數(shù)包括：群推分歧時(shí)間（也可以自己提供），基因流，以及有效群體大小。
（4）計(jì)算軟件：fastsimcoal2
fastsimcoal2是由伯爾尼大學(xué)的Laurent Excoffier小組2016年開(kāi)發(fā)的一種非常靈活的人口統(tǒng)計(jì)(Demography)建模軟件。 它通過(guò)執(zhí)行合并模擬，使用位點(diǎn)頻譜（SFS），推斷最適合所觀察數(shù)據(jù)的模型參數(shù)。
軟件下載
官網(wǎng)：http://cmpg.unibe.ch/software/fastsimcoal2/
使用fastsimcoal2和模擬數(shù)據(jù)在簡(jiǎn)單模型下推斷參數(shù)
輸入文件
3. fastsimcoal2實(shí)戰(zhàn)分析
使用fastsimcoal2和軟件自帶的測(cè)試數(shù)據(jù)在簡(jiǎn)單模型下推斷參數(shù)
（1）準(zhǔn)備輸入文件：
輸入文件1：SFS文件
文件2PopDiv20Mb_jointDAFpop1_0.obs是兩個(gè)群體的觀測(cè)2D-SFS.
# set the working directory
setwd("/Software/fsc26_mac64/example files/2PopDiv20Mb")
# read the observed SFS
# 2D from two populations with different effective sizes that diverged some time ago
pop2d <- read.table("2PopDiv20Mb_jointDAFpop1_0.obs", skip=1, stringsAsFactors = F, header=T, row.names = 1)
head(pop2d)
d0_0 d0_1 d0_2 d0_3 d0_4 d0_5
d1_0 19985747 8350 1628 360 62 8
d1_1 966 0 0 0 0 0
d1_2 479 0 0 0 0 0
d1_3 328 0 0 0 0 0
d1_4 249 0 0 0 0 0
d1_5 1760 13 18 13 19 0
輸入文件2：定義模型文件
每個(gè)模型都在TPL文件中定義，我們要推斷的參數(shù)都有名稱(chēng)標(biāo)簽（例如NPOP，TDIV）
//Parameters for the coalescence simulation program : fsimcoal2.exe
2 samples to simulate :
//Population effective sizes (number of genes)
NPOP1
NPOP2
//Samples sizes and samples age
5
5
//Growth rates : negative growth implies population expansion
0
0
//Number of migration matrices : 0 implies no migration between demes
0
//historical event: time, source, sink, migrants, new deme size, new growth rate, migration matrix index
1 historical event
TDIV 1 0 1 RESIZE0 0 0
//Number of independent loci [chromosome]
1 0
//Per chromosome: Number of contiguous linkage Block: a block is a set of contiguous loci
1
//per Block:data type, number of loci, per generation recombination and mutation rates and optional parameters
FREQ 1 0 2.5e-8 OUTEXP
EST文件中定義了每個(gè)參數(shù)的搜索范圍。 對(duì)于每個(gè)參數(shù)，我們使用相應(yīng)的參數(shù)名稱(chēng)標(biāo)簽指定搜索范圍。
// Search ranges and rules file
// ****************************
[PARAMETERS]
//#isInt? #name #dist.#min #max
//all Ns are in number of haploid individuals
1 NPOP1 logunif 10 1e7 output
1 NPOP2 logunif 10 1e7 output
1 NANC logunif 10 1e7 output
1 TDIV unif 10 1e5 output
[RULES]
[COMPLEX PARAMETERS]
0 RESIZE0 = NANC/NPOP1 hide
注意：TPL和EST文件需要具有相同的文件名，但具有不同的擴(kuò)展名：[filename] .est和[filename] .tpl。
（2）運(yùn)行fastsimcoal2：
#參數(shù)估計(jì)（推測(cè)群體演化歷史）
fsc26 -t PopDiv_diff.tpl -e PopDiv_diff.est -d -0 -n 100000 -L 40 -s 0 -M -q -c 80
#還可以產(chǎn)生模擬群體數(shù)據(jù)（另一個(gè)功能）
#使用輸入?yún)?shù)文件中定義的參數(shù)值來(lái)模擬進(jìn)化模型下的數(shù)據(jù)
fsc26 -i test.par -n 100
#使用從先驗(yàn)隨機(jī)抽取的參數(shù)值來(lái)模擬進(jìn)化模型下的數(shù)據(jù)
fsc26 -t test.tpl -n 10 -e test.est -E 100
#使用在外部文件中定義的參數(shù)值來(lái)模擬進(jìn)化模型下的數(shù)據(jù)
fsc26 -t test.tpl -n 100 -f test.def
參數(shù)說(shuō)明：
-n: 模擬次數(shù)，該值應(yīng)大于100,000。建議使用200,000到1,000,000。
-L: 迭代次數(shù)，該值至少應(yīng)為20，建議使用50和100之間。
-M: 使用似然優(yōu)化推斷參數(shù)。
-d: 對(duì)于derived SFS使用-d，對(duì)于MAF SFS使用-m。
-0: 說(shuō)明觀察到的SFS中沒(méi)有單態(tài)位。
-q: 快速模式，不打印所有信息。
-C: 為具有至少1個(gè)SNP的所有輸入計(jì)算似然。如果指定-Cx，則所有少于x個(gè)SNP的條目將匯集在一起。當(dāng)觀察到SFS中很多位點(diǎn)SNP很少時(shí)，此選項(xiàng)用以避免過(guò)度擬合。
-c: 指定多線程的選項(xiàng)。 -c1 -B1用于單核，-c4 -B4用于4核。
（3）軟件輸出：
對(duì)我們最重要的三個(gè)文件：
* .bestlhoods：具有最大似然參數(shù)估計(jì)值和相應(yīng)似然性的文件。 這就是我們想要的-參數(shù)估計(jì)！
* _DAFpop0.txt：具有通過(guò)優(yōu)化過(guò)程中使似然性最大化的參數(shù)獲得的預(yù)期SFS的文件，對(duì)于檢查SFS是否合適是必需的。
* .simparam：文件中包含運(yùn)行模擬的設(shè)置示例，用于錯(cuò)誤時(shí)檢查。
查看結(jié)果
#讀取最大似然估計(jì)參數(shù)文件
maxlhoodEstimates <- read.table(paste(settings$pathTo_InputFile, "/",
settings$TPL_EST_file_tag, "/",
settings$TPL_EST_file_tag, ".bestlhoods",
sep=""), header=T)
#查看估計(jì)參數(shù)
maxlhoodEstimates
其中，MaxObsLhood指如果期望值與觀察到的SFS完全匹配，即期望的SFS是相對(duì)觀察到的SFS，則值越大。MaxEstLhood是根據(jù)模型參數(shù)估計(jì)的最大似然，擬合度越高，MaxObsLhood和MaxEstLhood之間的差異就越小。
#獲取觀察到的和預(yù)期的SFS
# Fit of the model expected SFS to the observed SFS
# Tag for the end of the observed SFS file
obsfileend <- "_jointDAFpop1_0"
# Read the observed SFS - SNP counts
obssfs <- read.table(paste(settings$pathTo_InputFile,
settings$TPL_EST_file_tag, obsfileend, ".obs",
sep=""), skip=1, stringsAsFactors = F, header=T)
# Read the expected SFS - PROBABILITIES
expsfs <- read.table(paste(settings$pathTo_InputFile,
settings$TPL_EST_file_tag, "/",
settings$TPL_EST_file_tag, obsfileend, ".txt",
sep=""), header=T)
# Plot the fit of the 2D SFS, including of the marginal 1D SFS
# the function plot2dSFS is defined in the utilfunctions.r
# you need to give as input the following arguments
# obsSFS : matrix with observed SFS (counts)
# expSFS : matrix with expected SFS (probabilities)
# xtag : string with the label of x-axis
# ytag : string with the label of y-axis
# minentry : number with the minimum entry in the SFS (all entries with less than this are pooled together)
plot2dSFS(obsSFS=obssfs, expSFS=expsfs, xtag="Pop2", ytag="Pop1", minentry=1)
#繪制模型
#maxL.par file 是一個(gè)最大似然估計(jì)參數(shù)文件
path_to_maxL_file <- paste(settings$pathTo_InputFile, settings$TPL_EST_file_tag, "/",settings$TPL_EST_file_tag, "_maxL", sep="")
parFileInterpreter(args=path_to_maxL_file, pop.names=c("Pop1","Pop2"), gentime=1, printPDF=FALSE)
參考：
