入門
在比較早之前就看過一篇簡書
[EEMS推測生境內(nèi)生物的基因流 - 簡書 (jianshu.com)](http://www.reibang.com/p/9ef469ed7c98)
介紹過EEMS這款軟件咽瓷，點了收藏召烂，直接擱置。最近才慢慢參悟該軟件的使用沼琉，踩過很多坑北苟，所以直接介紹一下EEMS。
乍練
首先上官網(wǎng)：https://github.com/dipetkov/eems
使用手冊在Documentation中打瘪，有更詳細(xì)的說明友鼻。
通曉
EMMS的安裝比較麻煩，每一個模塊都需要單獨編譯闺骚，這里只介紹我用到的三個模塊彩扔。
#EMMS 主體可以直接git clone 下載 略
#首先是軟件主命令 runeems_snps
#需要兩款其他軟件作為支持
#Eigen進(jìn)行線性代數(shù)計算
https://eigen.tuxfamily.org/index.php?title=Main_Page
#下載后解壓 安裝過程直接 粘貼下來：
Let's call this directory 'source_dir' (where this INSTALL file is).
Before starting, create another directory which we will call 'build_dir'.
Do:
cd build_dir
cmake source_dir
make install
The "make install" step may require administrator privileges.
You can adjust the installation destination (the "prefix")
by passing the -DCMAKE_INSTALL_PREFIX=myprefix option to cmake, as is
explained in the message that cmake prints at the end.
#Boost進(jìn)行隨機數(shù)生成和生境幾何形狀
https://www.boost.org/
(不記得了 可以google 一下 how to install boost)
#兩款軟件裝好后，找到/eems/runeems_snps/src Makefile
#修改依賴路徑 make 就好 例：
EIGEN_INC = /data/01/user152/software/eigen-3.2.10/eigeneigen/include/eigen3/
BOOST_LIB = /data/01/user152/software/boost/lib
BOOST_INC = /data/01/user152/software/boost/include
#然后是bed2diffs模塊的編譯
#需要先安裝 libplinkio (https://github.com/fadern/libplinkio)
#這個比較簡單僻爽，略
#同樣找到Makefile 修改依賴路徑（在src里） 例：
PLINKIO = /data/01/user152/software/plinkio
EXE = bed2diffs_v1
OBJ = bed2diffs_v1.o data.o
CXXFLAGS = -I/data/01/user152/software/plinkio/include -O3 -Wall -Werror -fopenmp
LDFLAGS = -lplinkio
#最后是plot模塊虫碉，跟著手冊裝就好
#遇到rgos裝不上的情況，去跟管理員py一下胸梆，讓root裝一下就ok了敦捧。
至此，軟件就搞定了碰镜。
小成
runeems_snp 需要三個主要輸入文件：
*.diffs 遺傳距離矩陣绞惦，可以通過VCF轉(zhuǎn)換
*.coord 采樣的坐標(biāo)，經(jīng)緯度就可
*.outer 棲息地坐標(biāo)洋措，同經(jīng)緯度
現(xiàn)在我們一個個準(zhǔn)備輸入文件：
.diff
此處坑較多济蝉，務(wù)必注意
#首先要確定VCF中的個體是不是每個都有采樣信息，沒有的話需要刪掉，外群也需要刪掉王滤。
#可以用vcftools 刪除個體
#我在做的過程中做過很多次贺嫂，只介紹跑通的這次的做法，要問我為什么這么做雁乡，別問第喳，問就是這么做能跑通
1. 需要將VCF處理為下種模樣
1. 第一列需要都改為1
2. POS 順延
3. ID 為snp1 順延
#之后運行PLINK
plink --vcf sample.vcf --allow-extra-chr --maf 0.05 --recode --out sample
plink --noweb --file sample --allow-extra-chr --make-bed --maf 0.05 --allele1234 --out sample
生成的smple.[bed/bim/fam]就是下一步的生成文件。
bed2diffs_v1 --bfile ./sample
會生成第一個輸入文件踱稍， sample.diffs
.coord
記錄個體采樣地理位置的文件：
'經(jīng)度' \t'緯度'
建議經(jīng)度在前曲饱，當(dāng)然緯度在前也可以，不過后面需要注意這個的前后順序
.outer
記錄棲息地位置的文件
說人話就是
把所有樣本包括進(jìn)去的一個地理范圍珠月，類似下圖
采樣地點都在范圍內(nèi)扩淀。
可以用此網(wǎng)站準(zhǔn)備改文件。
http://www.birdtheme.org/useful/v3tool.html
到此為止啤挎，最主要的三個文件也就準(zhǔn)備好了驻谆，還有一類可選的文件。詳細(xì)的可以閱讀手冊庆聘。
大成
輸入文件準(zhǔn)備好了就可以運行軟件了胜臊。是一個類似配置文件的形式。
#params-run.ini
datapath = ./data/sample # 輸入文件前綴的路徑
mcmcpath = ./out #輸出文件路徑
#gridpath = #剛提到的可選的輸入文件前綴路徑
nIndiv =48 #個體數(shù)
nSites =5343964 #位點數(shù)
nDemes =700 #grid 網(wǎng)格數(shù)700我覺得有點太慢了伙判。象对。。宴抚。
diploid =true #是否為二倍體
numMCMCIter =2000000
numBurnIter =1000000
numThinIter =9999
因為用到的是SNP數(shù)據(jù)织盼，三個輸入文件內(nèi)并不能體現(xiàn)位點數(shù)，因此我用的是真實的位點數(shù)酱塔，可不可以瞎編我也不知道。
Demes 我理解的是最后網(wǎng)格數(shù)的多少危虱，越多越慢羊娃，示例為300好像，我采用700埃跷，具體見下圖蕊玷。
剩下的三個建議就默認(rèn)吧
#運行
runeems_snps --params params-run.ini
圓滿
輸出文件都有了，最重要的一點就是可視化弥雹。
###這里手冊里介紹的很詳細(xì)垃帅，不多加贅述
library("rEEMSplots")
eems_results <- file.path("out-1") #輸出文件路徑
name_figures <- file.path("out") #輸出圖片文件名前綴
eems.plots(mcmcpath = eems_results,
plotpath = paste0(name_figures, "-output-PDFs"), #輸出為PDF
longlat = TRUE, # TRUE 指經(jīng)度在前 緯度在后 FALSE反之
out.png = FALSE, #輸出為PDF
#add.grid = TRUE, #可選輸入文件才有
#col.grid = "gray90", #可選輸入文件才有
#lwd.grid = 2, #可選輸入文件才有
#add.outline = TRUE, # 是否顯示邊界 （我覺得太丑了）
#col.outline = "blue",
#lwd.outline = 5,
add.demes = TRUE, #顯示demes
col.demes = "red",
pch.demes = 5,
min.cex.demes = 0.5,
max.cex.demes = 1.5
)
最后結(jié)果大致這樣：
大道
判斷結(jié)果好不好的重要標(biāo)準(zhǔn)是： 模擬的模型是否收斂
可視化后有一個圖能用來判斷。
doc里是這么解釋的：
The eems.plots function in the rEEMSplots package generates a posterior trace plot. (a) This MCMC
chain obviously has not converged. (b) There is no indication that the MCMC chain has not converged. This is
not quite the same as there is evidence that the MCMC chain has converged. (c) To be confident that EEMS has
converged, we can run the MCMC sampler for more iterations. (d) Alternatively, to be confident that EEMS has
converged, we can run the MCMC sampler several times, each time starting from a different randomly initialized
parameter state.
說人話呢就是 最差也得是b圖的樣子剪勿。
為了達(dá)到收斂贸诚，有兩種方法
1. 多跑幾次，選最好的
2. 用上一次的結(jié)果作為下一次的輸入，循環(huán)
我們就草率的選2哈酱固！
#params-pre.ini
datapath = ./data/sample
mcmcpath = ./out
prevpath = ./out-1
#gridpath =
nIndiv =48
nSites =5343964
nDemes =700
diploid =true
numMCMCIter =1000000
numBurnIter =0
numThinIter =9999
跑到能讓自己舒服的結(jié)果就好械念。
感謝 軟件作者 及 DumplingLucky 還有KZR
