SFS的基本理解、介紹的一篇文章：https://theg-cat.com/tag/joint-sfs/
如果不做算法，基本的原理這里應(yīng)該也差不多了，先跑個(gè)試試

fastsimcoal2.6 SFS群體歷史動態(tài)模擬教程：
https://speciationgenomics.github.io/fastsimcoal2/
評估SFS教程：
https://speciationgenomics.github.io/easysfs/
Stairway plot Github網(wǎng)站：
https://github.com/xiaoming-liu/stairway-plot-v2

1、評估SFS (estimate Site Allele Frequency)

使用ANGSD進(jìn)行SFS Estimation。這個(gè)軟件考慮了missing data和低depth位點(diǎn)。以bam文件為input。
https://github.com/ANGSD/angsd

安裝

cd ~/software
wget http://popgen.dk/software/download/angsd/angsd0.930.tar.gz
tar xf angsd0.930.tar.gz
cd htslib;make;cd ..
cd angsd
make HTSSRC=../htslib
cd ..

后面按照教程來就行。參照官網(wǎng)example：http://www.popgen.dk/angsd/index.php/SFS_Estimation

如果missing data不多，可以直接用easySFS。
https://github.com/isaacovercast/easySFS

python3 /home/software/easySFS/easySFS.py -i vcf \
-p ./pops_file.txt -a -f --preview

看一下有沒有很多frequency很低的值，需不需要down sample。

我這里沒有，所以直接對VCF中第一、第二個(gè)種群分別取20、16個(gè)sample（其實(shí)是10、8個(gè)個(gè)體，二倍體所以sample數(shù)量乘二）：

python3 /home/software/easySFS/easySFS.py -i vcf \
-p ./pops_file.txt -a -f --proj 20,16

這個(gè) ./pops_file.txt 是告訴軟件VCF里哪幾個(gè)個(gè)體是一個(gè)種群的文件。格式官網(wǎng)教程有。

不管哪種方法，最后希望得到的是這樣的：

.obs文件格式

這是SFS observation文件，后綴為.obs。一共應(yīng)該有2N+1列，N為樣本個(gè)體數(shù)。后一半都是0是因?yàn)檫@是folded SFS。因?yàn)槲也恢廊后w的祖先狀態(tài)序列。如何設(shè)置folded或者unfolded，幾個(gè)軟件的教程里都有。

2、Fastsimcoal2.6 ， SFS建模模擬種群動態(tài)歷史

Fastsimcoal2.6輸入文件

使用Fastsimcoal2.6，利用SFS模擬種群動態(tài)歷史需要三個(gè)輸入文件：

• .tpl 模版文件（參數(shù)文件）：用來指定歷史事件、migration、樣本大小、突變率重組率等等。注意，這里的突變率重組率都為【per generation】
• .est 參數(shù)評估文件（prior文件）：在這里對參數(shù)施加prior，指定參數(shù)之間關(guān)系等等。
• .obs SFS觀測文件：就是上一步生成的。

這三個(gè)文件的前綴必須相同，在這里為"6.GL1.folded.HN"。如果為unfolded SFS，SFS觀測文件后綴應(yīng)該是DAF之類的而不是MAF，具體得看說明書。

一般來說，有兩種構(gòu)建模型的方法：
1）可以構(gòu)建很多模型（可能3～10個(gè)）來對種群歷史建模，然后比較每個(gè)模型的lnL（也就是LRT檢驗(yàn)），看那個(gè)模型最合適。比如第一個(gè)模型是有效群體大小先上升后下降在上升，第二個(gè)是先下降在上升等等。對于單個(gè)群體（不考慮migration）的簡單模型，一般也就是設(shè)置幾個(gè)bottleneck時(shí)間點(diǎn)，擴(kuò)張時(shí)間點(diǎn)以及對應(yīng)的種群大小參數(shù)。
2）只構(gòu)建一個(gè)模型，限制“時(shí)間點(diǎn)”prior的范圍，但是放松Ne prior的范圍，并保持一致，讓算法自己探索種群在某個(gè)時(shí)間點(diǎn)是擴(kuò)張還是收縮。在RULES中只對時(shí)間做限制，而不對Ne做限制。這種就需要增加運(yùn)行的次數(shù)，讓模型達(dá)到收斂，比如設(shè)置-L 100。例如：
.est文件：

.tpl文件：

image.png

設(shè)置好配置文件以后，開跑：

fsc26 -t ${i}.GL1.folded.HN.tpl -e ${i}.GL1.folded.HN.est \
-m -0 -C 10 -n 100000 -L 40 -s 0 -M -q -c 5 &

# 對于folded SFS，記得設(shè)置 -m

生成文件里有l(wèi)ikelihood的評估，拿過來比較模型就行了。對最優(yōu)模型可以再跑bootstrap，說明書里有。

我個(gè)人覺得對于單個(gè)種群而言，用SFS不一定是好的選擇。PSMC從原理上更加靠譜一些。而fastsimcoal模型的優(yōu)勢在于其可以考慮多個(gè)種群的migration、融合與分裂、群體大小變化等等。個(gè)人見解，歡迎討論。

1月21日更新

fastsimcoal模型比對結(jié)果還可以，但是并不是很直觀，我試了下Stairway plot v2.1，感覺挺靠譜，也很簡單。
網(wǎng)站：
https://github.com/xiaoming-liu/stairway-plot-v2

很湊巧，這篇文章，stairway plot2是去年（2020）發(fā)布的

git clone或者通過他給的鏈接下下來，然后準(zhǔn)備一個(gè)配置文件，就可以開跑了

SFS就用我們之前得到的。

配置文件

參數(shù)配置README文件里面都有寫，SFS也直接貼里面就行，比fastsimcoal方便。

具體算法孰優(yōu)孰劣，得了解算法，再看看人家比對的文章。好像有文章說單個(gè)參數(shù)的簡單模型或者多個(gè)種群用fastsimcoal靠譜，單個(gè)種群可能stairway plot好一些。個(gè)人理解。

tips：如果報(bào)錯(cuò)了很可能是命名問題導(dǎo)致沒有找到j(luò)ava class的路徑，可以修改回原命名試試。

image.png