SFS的基本理解嘱蛋、介紹的一篇文章：https://theg-cat.com/tag/joint-sfs/
如果不做算法吃媒，基本的原理這里應(yīng)該也差不多了，先跑個(gè)試試
fastsimcoal2.6 SFS群體歷史動(dòng)態(tài)模擬教程：
https://speciationgenomics.github.io/fastsimcoal2/
評(píng)估SFS教程：
https://speciationgenomics.github.io/easysfs/
Stairway plot Github網(wǎng)站：
https://github.com/xiaoming-liu/stairway-plot-v2
1雀费、評(píng)估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 ..
后面按照教程來(lái)就行。參照官網(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
看一下有沒(méi)有很多frequency很低的值，需不需要down sample。
我這里沒(méi)有气筋，所以直接對(duì)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)教程有光稼。
不管哪種方法或南，最后希望得到的是這樣的：
這是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建模模擬種群動(dòng)態(tài)歷史
使用Fastsimcoal2.6蝴罪，利用SFS模擬種群動(dòng)態(tài)歷史需要三個(gè)輸入文件：
- .tpl 模版文件（參數(shù)文件）：用來(lái)指定歷史事件董济、migration、樣本大小要门、突變率重組率等等虏肾。注意，這里的突變率重組率都為【per generation】
- .est 參數(shù)評(píng)估文件（prior文件）：在這里對(duì)參數(shù)施加prior欢搜，指定參數(shù)之間關(guān)系等等封豪。
- .obs SFS觀測(cè)文件：就是上一步生成的。
這三個(gè)文件的前綴必須相同炒瘟，在這里為"6.GL1.folded.HN"吹埠。如果為unfolded SFS，SFS觀測(cè)文件后綴應(yīng)該是DAF之類的而不是MAF疮装，具體得看說(shuō)明書(shū)缘琅。
一般來(lái)說(shuō)，有兩種構(gòu)建模型的方法：
1）可以構(gòu)建很多模型（可能3～10個(gè)）來(lái)對(duì)種群歷史建模斩个，然后比較每個(gè)模型的lnL（也就是LRT檢驗(yàn)）胯杭，看那個(gè)模型最合適驯杜。比如第一個(gè)模型是有效群體大小先上升后下降在上升受啥，第二個(gè)是先下降在上升等等。對(duì)于單個(gè)群體（不考慮migration）的簡(jiǎn)單模型，一般也就是設(shè)置幾個(gè)bottleneck時(shí)間點(diǎn)滚局，擴(kuò)張時(shí)間點(diǎn)以及對(duì)應(yīng)的種群大小參數(shù)居暖。
2）只構(gòu)建一個(gè)模型，限制“時(shí)間點(diǎn)”prior的范圍藤肢，但是放松Ne prior的范圍太闺，并保持一致，讓算法自己探索種群在某個(gè)時(shí)間點(diǎn)是擴(kuò)張還是收縮嘁圈。在RULES中只對(duì)時(shí)間做限制省骂，而不對(duì)Ne做限制。這種就需要增加運(yùn)行的次數(shù)最住，讓模型達(dá)到收斂钞澳，比如設(shè)置-L 100。例如：
.est文件：
.tpl文件：
設(shè)置好配置文件以后涨缚，開(kāi)跑：
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 &
# 對(duì)于folded SFS轧粟，記得設(shè)置 -m
生成文件里有l(wèi)ikelihood的評(píng)估，拿過(guò)來(lái)比較模型就行了脓魏。對(duì)最優(yōu)模型可以再跑bootstrap兰吟，說(shuō)明書(shū)里有。
我個(gè)人覺(jué)得對(duì)于單個(gè)種群而言茂翔，用SFS不一定是好的選擇混蔼。PSMC從原理上更加靠譜一些。而fastsimcoal模型的優(yōu)勢(shì)在于其可以考慮多個(gè)種群的migration珊燎、融合與分裂拄丰、群體大小變化等等。個(gè)人見(jiàn)解俐末，歡迎討論料按。
1月21日更新
fastsimcoal模型比對(duì)結(jié)果還可以，但是并不是很直觀卓箫，我試了下Stairway plot v2.1载矿，感覺(jué)挺靠譜，也很簡(jiǎn)單烹卒。
網(wǎng)站：
https://github.com/xiaoming-liu/stairway-plot-v2
很湊巧闷盔，這篇文章，stairway plot2是去年（2020）發(fā)布的
git clone或者通過(guò)他給的鏈接下下來(lái)旅急，然后準(zhǔn)備一個(gè)配置文件逢勾，就可以開(kāi)跑了
SFS就用我們之前得到的。
參數(shù)配置README文件里面都有寫(xiě)藐吮，SFS也直接貼里面就行溺拱，比f(wàn)astsimcoal方便逃贝。
具體算法孰優(yōu)孰劣，得了解算法迫摔，再看看人家比對(duì)的文章沐扳。好像有文章說(shuō)單個(gè)參數(shù)的簡(jiǎn)單模型或者多個(gè)種群用fastsimcoal靠譜，單個(gè)種群可能stairway plot好一些句占。個(gè)人理解沪摄。
tips：如果報(bào)錯(cuò)了很可能是命名問(wèn)題導(dǎo)致沒(méi)有找到j(luò)ava class的路徑，可以修改回原命名試試纱烘。
