入門

在比較早之前就看過一篇簡書

[EEMS推測生境內(nèi)生物的基因流 - 簡書 (jianshu.com)](http://www.itdecent.cn/p/9ef469ed7c98)

介紹過EEMS這款軟件，點(diǎn)了收藏，直接擱置。最近才慢慢參悟該軟件的使用，踩過很多坑，所以直接介紹一下EEMS。

乍練

首先上官網(wǎng)：https://github.com/dipetkov/eems
使用手冊在Documentation中，有更詳細(xì)的說明。

通曉

EMMS的安裝比較麻煩，每一個模塊都需要單獨(dú)編譯，這里只介紹我用到的三個模塊。

#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)行隨機(jī)數(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 順延

#之后運(yùn)行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)去的一個地理范圍，類似下圖

采樣地點(diǎn)都在范圍內(nèi)。
可以用此網(wǎng)站準(zhǔn)備改文件。

http://www.birdtheme.org/useful/v3tool.html

到此為止，最主要的三個文件也就準(zhǔn)備好了，還有一類可選的文件。詳細(xì)的可以閱讀手冊。

大成

輸入文件準(zhǔn)備好了就可以運(yùn)行軟件了。是一個類似配置文件的形式。

#params-run.ini
datapath = ./data/sample  # 輸入文件前綴的路徑
mcmcpath = ./out   #輸出文件路徑
#gridpath =   #剛提到的可選的輸入文件前綴路徑
nIndiv =48  #個體數(shù)
nSites =5343964  #位點(diǎn)數(shù)
nDemes =700 #grid 網(wǎng)格數(shù)700我覺得有點(diǎn)太慢了。。。。
diploid =true  #是否為二倍體
numMCMCIter =2000000
numBurnIter =1000000
numThinIter =9999

因?yàn)橛玫降氖荢NP數(shù)據(jù)，三個輸入文件內(nèi)并不能體現(xiàn)位點(diǎn)數(shù)，因此我用的是真實(shí)的位點(diǎn)數(shù)，可不可以瞎編我也不知道。
Demes 我理解的是最后網(wǎng)格數(shù)的多少，越多越慢，示例為300好像，我采用700，具體見下圖。
剩下的三個建議就默認(rèn)吧

#運(yùn)行
runeems_snps --params params-run.ini

圓滿

輸出文件都有了，最重要的一點(diǎn)就是可視化。

###這里手冊里介紹的很詳細(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