参悟 Estimated Effective Migration Surface(EEMS)

入门

在比较早之前就看过一篇简书

[EEMS推测生境内生物的基因流 - 简书 (jianshu.com)](https://www.jianshu.com/p/9ef469ed7c98)

介绍过EEMS这款软件,点了收藏,直接搁置。最近才慢慢参悟该软件的使用,踩过很多坑,所以直接介绍一下EEMS。

乍练

首先上官网:https://github.com/dipetkov/eems
使用手册在Documentation中,有更详细的说明。

通晓

EMMS的安装比较麻烦,每一个模块都需要单独编译,这里只介绍我用到的三个模块。

#EMMS 主体可以直接git clone 下载 略
#首先是软件主命令 runeems_snps
#需要两款其他软件作为支持 
#Eigen进行线性代数计算
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进行随机数生成和生境几何形状
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转换
*.coord 采样的坐标,经纬度就可
*.outer 栖息地坐标,同经纬度
现在我们一个个准备输入文件:

.diff

此处坑较多,务必注意
#首先要确定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

记录个体采样地理位置的文件:

'经度' \t'纬度'
建议经度在前,当然纬度在前也可以,不过后面需要注意这个的前后顺序

.outer

记录栖息地位置的文件
说人话就是
把所有样本包括进去的一个地理范围,类似下图



采样地点都在范围内。
可以用此网站准备改文件。

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

到此为止,最主要的三个文件也就准备好了,还有一类可选的文件。详细的可以阅读手册。

大成

输入文件准备好了就可以运行软件了。是一个类似配置文件的形式。

#params-run.ini
datapath = ./data/sample  # 输入文件前缀的路径
mcmcpath = ./out   #输出文件路径
#gridpath =   #刚提到的可选的输入文件前缀路径
nIndiv =48  #个体数
nSites =5343964  #位点数
nDemes =700 #grid 网格数700我觉得有点太慢了。。。。
diploid =true  #是否为二倍体
numMCMCIter =2000000
numBurnIter =1000000
numThinIter =9999
因为用到的是SNP数据,三个输入文件内并不能体现位点数,因此我用的是真实的位点数,可不可以瞎编我也不知道。
Demes 我理解的是最后网格数的多少,越多越慢,示例为300好像,我采用700,具体见下图。
剩下的三个建议就默认吧
#运行
runeems_snps --params params-run.ini

圆满

输出文件都有了,最重要的一点就是可视化。

###这里手册里介绍的很详细,不多加赘述
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 指经度在前 纬度在后 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
           )

最后结果大致这样:


大道

判断结果好不好的重要标准是: 模拟的模型是否收敛
可视化后有一个图能用来判断。
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图的样子。
为了达到收敛,有两种方法

1. 多跑几次,选最好的
2. 用上一次的结果作为下一次的输入,循环
我们就草率的选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

跑到能让自己舒服的结果就好。
感谢 软件作者 及 DumplingLucky 还有KZR

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
禁止转载,如需转载请通过简信或评论联系作者。
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 204,921评论 6 478
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 87,635评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 151,393评论 0 338
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,836评论 1 277
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,833评论 5 368
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,685评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,043评论 3 399
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,694评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 42,671评论 1 300
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,670评论 2 321
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,779评论 1 332
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,424评论 4 321
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,027评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,984评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,214评论 1 260
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,108评论 2 351
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,517评论 2 343

推荐阅读更多精彩内容