dipetkov/eems

Estimating Effective Migration Surfaces

EEMS

This repository contains an implementation of the EEMS method for analyzing and visualizing spatial population structure from geo-referenced genetic samples. EEMS uses the concept of effective migration to model the relationship between genetics and geography, and it outputs an estimated effective migration surface (hence, EEMS) - a visual representation of population structure that can highlight potential regions of higher-than-average and lower-than-average historic gene flow.

Please consider reading the paper EEMS-article.pdf (or the published version here) and the documentation EEMS-doc.pdf first.

• The directory bed2diffs contains a small program that reads genotypes in Plink binary format and computes the matrix of average pairwise differences, to pass to EEMS.
• The directory pipeline contains a pipeline with the goal to combine all scripts.
• The directory plotting contains an R package, rEEMSplots, to visualize EEMS results.
• The directory runeems_sats contains a C++ implementation of EEMS to use with microsatellite data.
• The directory runeems_snps contains a C++ implementation of EEMS to use with SNP data.

C++ implementation

The C++ implementation uses the Eigen template library for linear algebra computations and the Boost libraries for random number generation and the habitat geometry. EEMS has been tested with Eigen 3.2.2 3.3.9 and Boost 1.57 1.76.0 and might not be compatible with newer versions of Boost/Eigen. After downloading Eigen (which does not need installation) and installing Boost, update the variables EIGEN_INC, BOOST_LIB, BOOST_INC in the Makefile. The dynamic Boost libraries are linked to slightly differently on a Mac and a Linux machine, so run make darwin on a Mac or make linux on a Linux machine. This C++ implementation has not been tested on Windows.

EEMS for SNPs and microsatellites

There are two versions of EEMS: runeems_snps for SNP data and runeems_sats for microsatellite data. The data input format and the EEMS model are somewhat different for SNPs and microsatellites, hence the two versions. The source code can be found in runeems_snps/src and runeems_sats/src, respectively. The directories runeems_snps/data and runeems_sats/data contain data, simulated with ms, to illustrate the input file format and how EEMS is run.

For SNP data, the input files are:

• datapath.diffs: the matrix of average pairwise genetic differences (which can be computed with the program bed2diffs)

• datapath.coord: the sample coordinates (two coordinates per sample, one sample per line)

• datapath.outer: the habitat coordinates (as a sequence of vertices that outline a closed polygon)

Here datapath is the full path + the file name (but without the extension).

And for microsatellite data, the input files are:

• datapath.sites: genotype data (one sample per line, a negative number indicates the allele is missing)

• datapath.coord: as described above for SNP data.

• datapath.outer: as described above for SNP data.

EEMS also requires a configuration file where various program options can be specified. For example, consider runeems_snps/src/params-simno1.ini which contains the following information:

datapath = ./data/barrier-schemeZ-nIndiv300-nSites3000
mcmcpath = ./data/barrier-schemeZ-nIndiv300-nSites3000-EEMS-nDemes200-simno1
nIndiv = 300
nSites = 3000
nDemes = 200
diploid = false
numMCMCIter = 2000000
numBurnIter = 1000000
numThinIter = 9999

This file specifies the following required input arguments: the path to the input data (datapath), the path to the output data (mcmcpath), the number of samples (nIndiv), the number of markers (nSites), the density of the population grid (nDemes), is the species diploid or haploid (diploid), the number of MCMC and burn-in iterations (numMCMCIter, numBurnIter), and the thinning interval (numThinIter).

./runeems_snps --params params-simno1.ini --seed 123

(Specifying the random seed is optional.)

Finally, the EEMS results can be visualized with the function eems.plots defined in the R package rEEMSplot. The package is not on CRAN, so install it from source instead. (The code is in the directory plotting.)

## Part 1: Install rEEMSplots
## Check that the current directory contains the rEEMSplots source directory
if (file.exists("./rEEMSplots")) {
  install.packages("rEEMSplots", repos = NULL, type = "source")
} else {
  stop("Move to the directory that contains the rEEMSplots source to install the package.")
}


## Possibly change the working directory with setwd()


## Part 2: Generate graphics
library(rEEMSplots)

mcmcpath = "./data/barrier-schemeX-nIndiv300-nSites3000-EEMS-nDemes200-simno1"
plotpath = "./plot/barrier-schemeX-nIndiv300-nSites3000-EEMS-nDemes200-simno1-rEEMSplots"

eems.plots(mcmcpath, plotpath, longlat = TRUE)

The function eems.plots generates several figures automatically (to encourage looking at all the figures). There are examples in EEMS-doc.pdf, with captions that explain each figure. eems.plots also saves several objects to an RData file, which can be read back from the file with load.

load(paste0(plotpath, "-rdist.RData"))

ls()
#> [1] "B.component" "G.component" "W.component" "xym.values"  "xyq.values"

library("ggplot2")
library("dplyr")

## Reproduce plotpath-rdist01.png,
## which plots observed vs fitted dissimilarities between demes
ggplot(B.component %>% filter(size > 1),
       aes(fitted, obsrvd)) +
  geom_point()

The data frame B.component contains information about the dissimilarities between observed pairs of demes (alpha, beta). "Observed" means that each deme has at least one sampled individual assigned to it. Furthermore, each deme has two coordinates (x and y, longitude and latitude) and these are labeled alpha.x, alpha.y for deme alpha and beta.x, beta.y for deme beta.

## Reproduce plotpath-rdist02.png,
## which plots observed vs fitted dissimilarities within demes
ggplot(W.component %>% filter(size > 1),
       aes(fitted, obsrvd)) +
  geom_point()

## Reproduce plotpath-rdist03.png,
## which plots observed dissimilarities against great circle distances between demes
ggplot(W.component %>% filter(size > 1),
       aes(fitted, obsrvd)) +
  geom_point()

## Empty matrices unless you have specified additional coordinates
## at which to estimate the migration and diversity rates
## with the `xy.coords` argument to `eems.plots`
xym.values
xyq.values

Shiny app

It is straightforward to make the scatter plots which plot fitted vs observed dissimilarities between demes. But it is even more useful to link the dissimilarities to the corresponding locations in the habitat. The between-demes dissimilarities saved in plotpath-rdist.RData can be visualized with a little interactive Shiny app available here. It can help to identify outliers in the pairwise scatter plots, if any, with specific geographic locations on the map.

The MATLAB/Octave implementation

This version is less efficient and it is provided here only for completeness, as it is the original implementation.

eems's Issues

Combining results from different runs

I have a question about combining the results from different runs. The manual and the manuscript both recommend running the program with different numbers of demes, and then averaging the results from those different runs. However, I can't seem to figure out how to average different runs. How would I go about doing this?

Thanks!

Segfault error in rEEMSplots using example data

Hi, I have successfully installed, but have not been able to run the R package rEEMSplots. Whether I use data files generated myself by running runeems_snps, or if I use the example files that come with the rEEMSplots package, I always encounter a segmentation fault.

I do not understand what the error message means or how to fix it. Searching online said this type of error might indicate a problem with the code. I am using R version 3.4.0 on Mac OS X 10.11.6 El Capitan.

Here is an example of the error output to the R console using the example scripts and example files provided with the rEEMSplots documentation. Any help or suggestions would be greatly appreciated. Thank you very much! - Jordan

The error message:

library(rEEMSplots)
eems.results.to.plot = paste(path.package("rEEMSplots"),"/extdata/EEMS-example",sep="")
name.figures.to.save = "EEMS-example-rEEMSplots"

eems.plots(mcmcpath = eems.results.to.plot, plotpath = paste(name.figures.to.save,"-default",sep=""), longlat = TRUE)

Input projection:
Output projection:

Using the default DarkOrange to Blue color scheme, with 'white' as the midpoint color.
It combines two color schemes from the 'dichromat' package, which itself is based on
a collection of color schemes for scientific data graphics:
Light A and Bartlein PJ (2004). The End of the Rainbow? Color Schemes for Improved Data
Graphics. EOS Transactions of the American Geophysical Union, 85(40), 385.
See also http://geog.uoregon.edu/datagraphics/color_scales.htm

Using 'euclidean' distance to assign interpolation points to Voronoi tiles.

Processing the following EEMS output directories :
/Library/Frameworks/R.framework/Versions/3.4/Resources/library/rEEMSplots/extdata/EEMS-example
Plotting effective migration surface (posterior mean of m rates)
/Library/Frameworks/R.framework/Versions/3.4/Resources/library/rEEMSplots/extdata/EEMS-example

*** caught segfault ***
address 0x18, cause 'memory not mapped'

Traceback:
1: .Call("rEEMSplots_tiles2contours_standardize", PACKAGE = "rEEMSplots", tiles, rates, seeds, marks, distm)
2: tiles2contours_standardize(tiles, rates, cbind(xseed, yseed), marks, dimns$dist.metric)
3: transform.rates(dimns, tiles, rates, xseed, yseed, zero_mean)
4: average.eems.contours(mcmcpath, dimns, longlat, plot.params, is.mrates = TRUE, plot.xy = m.plot.xy)
5: eems.plots(mcmcpath = eems.results.to.plot, plotpath = paste(name.figures.to.save, "-default", sep = ""), longlat = TRUE)

Possible actions:
1: abort (with core dump, if enabled)
2: normal R exit
3: exit R without saving workspace
4: exit R saving workspace
Selection:

The diffs matrix is not full rand.

Dear author,

I met a problem that shows the diffs matrix is not full rank. Because my snps matrix only included the mtDNA bi-allelic SNPs, and they were just converted into the 0|0 or 1|1 format, I'm not sure whether it was related to this. if not, how to solve the non-full rank diff matrix issue? Could you please help? Thank you.

diffs not a full-rank matrix

Dr. Petkova, I am trying to run eems_snps and encountered an error:

[Diffs::initialize]
Error: The dissimilarity matrix is not a full-rank distance matrix

I am wondering if you could help me resolve the issue?

I created the.diffs file using bed2diffs, and see that the file looks good (nrows = ncols, diagonals = 0).
I estimated the rank using rankMatrix from the 'Matrix' R package, and found it is indeed different (slightly) from 0 (rank = 2.065015e-14). I checked the rank of a matrix I had successfully used in the past, but it was also different from 0 (rank = 4.174439e-14).
The only difference between the matrices is that the matrix that worked was computed with ~4x more SNPs.
Is there a simple way to modify the matrix so that eems_snps will accept it, without modifying eems results?

bed2diff error: [Data::getsize] Error opening plink files

Dear Dr Petkova,

Using bed2diffs on my Plink files (.bed, .fam and .bim, all at the same directory), I got the following message:

./bed2diffs_v1 --bfile palma2 --nthreads 4
Compute the average genetic differences according to:
Dij = (1/|Mij|) sum_{m in Mij} (z_{im} - z_{jm})^2
where Mij is the set of SNPs where both i and j are called

[Data::getsize] Error opening plink files palma2.[bed/bim/fam]

Please see the input files attached. I am wondering what might be wrong with the dataset, which is extensive (it's a product from GBS data) but seems alright. Could you please help me with that?
Thanks
test.zip

Issue while plotting my data in R

Hi,

When trying to plot my graph I get this R error: "Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function ‘coordinates’ for signature ‘"numeric". What could lead to this error and how to fix this?

additional info that doesn't seem right when running the eems analysis in the terminal:
There are 1 observed demes (out of 172 demes)
The population grid has 172 demes and 444 edges
There are 792 samples assigned to 1 observed demes

Many thanks

Error in make_eems_plots

I'm able to run the program and visualize the results fine for the most part, however when I try to use m.plot.xy or q.plot.xy to add graphical elements to the plots like so

make_eems_plots(mcmcpath = eems_results,
                longlat = TRUE,
                m.plot = { points(coord, col = "purple") },
                q.plot = { points(coord, col = "purple") }) 

I get:

Error in make_eems_plots(mcmcpath = eems_results, longlat = TRUE, m.plot.xy = { : 
  unused arguments (m.plot.xy = {
    points(coord, col = "purple")
}, q.plot.xy = {
    points(coord, col = "purple")
})

Any suggestions?

Thanks,
Tom

EEMS on ANGSD's ibsMat

Hi Desislava - I am trying to make EEMS work on the IBS-distance matrix from ANGSD ( http://www.popgen.dk/angsd/index.php/PCA_MDS ), and I am having a funny problem: with increasing number of individuals (>300) the determinant of the IBS matrix becomes essentially zero and runeems throws an error. Is there a way around it, do you think?
cheers
Misha

Complex habitats

Hi -

I'm running EEMS using a species has a long, narrow distribution. We get different results when we use a habitat that approximately follows the coastline (complex) compared to one where we smooth the habitat perimeter (simple). It appears that one inlet in the complex habitat is causing different inferences of effective migration rate. This holds true across a number of different grid sizes. I'm wondering why this is and if you could provide any further insight.

Here's the simple habitat:

Here's the complex habitat:

Improving the plotting script

I played around a bit with the idea of overlaying mean and variance in one plot. I ended up with ggplot. I currently have a working example, but it does not look nearly as nice as your figures:

See https://github.com/BenjaminPeter/eems/tree/master/runeems_snps/plot and the tmp*pdf plots for what I ended up with.

Also, just to re-iterate map projections from my email: They can be used with the mapproj package. But everything has to be passed through a function to transform coordinates.

require(mapproj)
require(maps)

m <-function(x,y){
   return(mapproject(x,y,
                      orientation=c(90,150,0),
                      projection="mercator"))
}
limits <- m(c(-20,-40), c(-50,70))
plot(NA, xlim=limits$x, ylim=limits$y)
map(orientation=c(90,150,0),projection="mercator", wrap=T, add=T)

pt <- m(-87.6847, 41.8369)
points(pt$x, pt$y,col="red", pch=".", cex=10)

Problem compiling runeems_sats

Hi

On a Red Hat 6 server with boost 1.59 and eigen3.2.5.2 installed, I'm trying to compile runeems_sats, checkout a328992.
I get the following error:

$ make linux
g++ -I/usr/local/boost_1.59/include -I/usr/include/eigen3 -O3 -DNDEBUG -c runeems_sats.cpp -o runeems_sats.o
g++ -I/usr/local/boost_1.59/include -I/usr/include/eigen3 -O3 -DNDEBUG -c eems.cpp -o eems.o
eems.cpp: In member function ‘void EEMS::initialize_diffs()’:
eems.cpp:64: error: ‘to_string’ was not declared in this scope
make: *** [eems.o] Error 1

Compiling runeems_snps worked without problem using the same compiler, paths to boost, eigen3 etc.

$ g++ --version
g++ (GCC) 4.4.7 20120313 (Red Hat 4.4.7-16)

Any idea what the problem might be? I can provide more info if required.

Thank you and best regards
Chris

Question on karyotype

Good day,
I'm writing you because I've got an issue with processing a dataset with eems. I have to analyse a dataset of livestock with a number of chromosomes larger than the human karyotype. This creates a problem with libplinkio, which cannot calculate the difference matrix from my bed file. Is there a way to work around the problem? For example, converting the chromosome number to 0 will somehow affect the results?
Thank you in advance for your help
Andrea

Errors in compiling runeems_snps

Hi Dr. Petkova, I am trying to implement EEMS using Eigen 3.2.2 and Boost 1.57.

After installing Eigen and Boost, I face the error when trying to link them to the EEMS program (scroll to bottom for full error message):

ld: symbol(s) not found for architecture x86_64
clang: error: linker command failed with exit code 1 (use -v to see invocation)
make: *** [darwin] Error 1

I have seen similar questions and have tried various solutions: changing order of libraries/programs during compilation; change from libc++ to C++ Standard Library libstdc++ to build boost; adding -lboost_serialization, but none have worked. I appreciate your help and advice.

Installation steps:

1. Download EEMS
2. Download Eigen 3.2.2 & copy to /usr/local/include
3. Download & install Boost 1_57

cd /Users/fuzzywuzzy/Downloads
tar --bzip2 -xf boost_1_57_0.tar.bz2 #extract boost files
cd /Users/fuzzywuzzy/Downloads/boost_1_57_0 #go into the new folder
./bootstrap.sh #create the Boost build utilities
./b2 install --prefix=/usr/local #build Boost libraries

4. Update the variables EIGEN_INC, BOOST_INC, BOOST_LIB in the EEMS Makefile
5. Link dependencies to program

cd /Users/fuzzywuzzy/Downloads/eems-master/runeems_snps/src/ #enter source code folder
make darwin    
c++ -I/usr/local -I/usr/local/eigen -O3 -DNDEBUG -o runeems_snps runeems_snps.o eems.o util.o mcmc.o draw.o habitat.o graph.o -L/usr/local/bin/lib -lboost_system -lboost_program_options -lboost_filesystem

This is where the error occurs - error message below:

Undefined symbols for architecture x86_64:
  "boost::program_options::to_internal(std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > const&)", referenced from:
      std::__1::vector<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, std::__1::allocator<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > > > boost::program_options::to_internal<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > >(std::__1::vector<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, std::__1::allocator<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > > > const&) in runeems_snps.o
  "boost::program_options::validation_error::get_template(boost::program_options::validation_error::kind_t)", referenced from:
      std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > const& boost::program_options::validators::get_single_string<char>(std::__1::vector<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, std::__1::allocator<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > > > const&, bool) in runeems_snps.o
  "boost::program_options::basic_parsed_options<char> boost::program_options::parse_config_file<char>(std::__1::basic_istream<char, std::__1::char_traits<char> >&, boost::program_options::options_description const&, bool)", referenced from:
      _main in runeems_snps.o
  "boost::program_options::options_description::options_description(std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > const&, unsigned int, unsigned int)", referenced from:
      _main in runeems_snps.o
  "boost::program_options::invalid_option_value::invalid_option_value(std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > const&)", referenced from:
      void boost::program_options::validate<long, char>(boost::any&, std::__1::vector<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, std::__1::allocator<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > > > const&, long*, long) in runeems_snps.o
      void boost::program_options::validate<int, char>(boost::any&, std::__1::vector<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, std::__1::allocator<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > > > const&, int*, long) in runeems_snps.o
      void boost::program_options::validate<double, char>(boost::any&, std::__1::vector<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, std::__1::allocator<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > > > const&, double*, long) in runeems_snps.o
  "boost::program_options::error_with_option_name::error_with_option_name(std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > const&, std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > const&, std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > const&, int)", referenced from:
      std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > const& boost::program_options::validators::get_single_string<char>(std::__1::vector<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, std::__1::allocator<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > > > const&, bool) in runeems_snps.o
  "boost::program_options::arg", referenced from:
      boost::program_options::typed_value<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, char>::name() const in runeems_snps.o
      boost::program_options::typed_value<long, char>::name() const in runeems_snps.o
      boost::program_options::typed_value<int, char>::name() const in runeems_snps.o
      boost::program_options::typed_value<bool, char>::name() const in runeems_snps.o
      boost::program_options::typed_value<double, char>::name() const in runeems_snps.o
  "boost::program_options::detail::cmdline::set_additional_parser(boost::function1<std::__1::pair<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > >, std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > const&>)", referenced from:
      boost::program_options::basic_command_line_parser<char>::extra_parser(boost::function1<std::__1::pair<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > >, std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > const&>) in runeems_snps.o
  "boost::program_options::detail::cmdline::cmdline(std::__1::vector<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, std::__1::allocator<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > > > const&)", referenced from:
      boost::program_options::basic_command_line_parser<char>::basic_command_line_parser(int, char const* const*) in runeems_snps.o
  "boost::program_options::validate(boost::any&, std::__1::vector<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, std::__1::allocator<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > > > const&, std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >*, int)", referenced from:
      boost::program_options::typed_value<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, char>::xparse(boost::any&, std::__1::vector<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, std::__1::allocator<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > > > const&) const in runeems_snps.o
  "boost::program_options::validate(boost::any&, std::__1::vector<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, std::__1::allocator<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > > > const&, bool*, int)", referenced from:
      boost::program_options::typed_value<bool, char>::xparse(boost::any&, std::__1::vector<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, std::__1::allocator<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > > > const&) const in runeems_snps.o
  "boost::program_options::operator<<(std::__1::basic_ostream<char, std::__1::char_traits<char> >&, boost::program_options::options_description const&)", referenced from:
      _main in runeems_snps.o
  "boost::program_options::error_with_option_name::substitute_placeholders(std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > const&) const", referenced from:
      vtable for boost::exception_detail::clone_impl<boost::exception_detail::error_info_injector<boost::program_options::validation_error> > in runeems_snps.o
      vtable for boost::exception_detail::error_info_injector<boost::program_options::validation_error> in runeems_snps.o
      vtable for boost::program_options::validation_error in runeems_snps.o
      vtable for boost::exception_detail::clone_impl<boost::exception_detail::error_info_injector<boost::program_options::invalid_option_value> > in runeems_snps.o
      vtable for boost::exception_detail::error_info_injector<boost::program_options::invalid_option_value> in runeems_snps.o
      vtable for boost::program_options::invalid_option_value in runeems_snps.o
  "boost::program_options::value_semantic_codecvt_helper<char>::parse(boost::any&, std::__1::vector<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, std::__1::allocator<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > > > const&, bool) const", referenced from:
      vtable for boost::program_options::typed_value<std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, char> in runeems_snps.o
      vtable for boost::program_options::typed_value<long, char> in runeems_snps.o
      vtable for boost::program_options::typed_value<int, char> in runeems_snps.o
      vtable for boost::program_options::typed_value<bool, char> in runeems_snps.o
      vtable for boost::program_options::typed_value<double, char> in runeems_snps.o
ld: symbol(s) not found for architecture x86_64
clang: error: linker command failed with exit code 1 (use -v to see invocation)
make: *** [darwin] Error 1

Thank you!

Error in rEEMSplots

Hello,

I am trying to work through the example to get a handle on using the program.

I had issues installing the plotting package following the instructions in the manual. But was able to successfully install the package using the following
library(devtools)
install_github("dipetkov/eems/plotting/rEEMSplots")

When I tried making plots with one of the provided example data sets I get the following....
Error in bitmap(paste0(plotpath, "%02d.png"), type = "png16m", res = plot.params$res, :
GhostScript was not found

I have since installed GhostScript on my computer but was wondering if you could provide advice on how to get R to recognize it as I still get the error when trying to make pngs. I haven't been able to find a good solution through a google search. I am however able to create pdfs only.

Thank you

Error: Check that the genotype matrix is a ..

Hi, trying to use runeems_sats and am getting "Error: Check that the genotype matrix is a nIndiv x 128 matrix. I'm confused by this. Shouldn't the matrix be nIndiv x nLoci (=n x 2L for L diploid loci)? Pasted below is the output:

Tom-Devitts-MBP:src tom$ ./runeems_sats --params params-chain1.ini
Using Boost 1_61 (EEMS was tested with version 1_57)
Eigen 3.2.9 (EEMS was tested with version 3.2.2)

Will save EEMS output to ../data/gopher_eems_200demes_chain1/
[Habitat::initialize]
Loaded habitat points from ../data/gopher_eems.outer
Input habitat:
POLYGON((30.9264 -87.7875,30.7943 -87.7436,30.6999 -87.6777,30.5108 -87.6777,30.3592 -87.6777,30.2454 -87.5019,30.2264 -87.3261,30.2644 -87.1283,30.3023 -86.9306,30.3023 -86.7768,30.3023 -86.5351,30.3023 -86.2494,30.1695 -86.0077,29.9603 -85.832,29.7125 -85.6342,29.5979 -85.3925,29.5979 -85.0849,29.617 -84.997,29.6552 -84.8652,29.7125 -84.5356,29.7888 -84.4477,29.8269 -84.2939,29.8841 -84.206,29.9222 -84.0302,29.8841 -83.9423,29.846 -83.8984,29.8079 -83.8324,29.7316 -83.7445,29.5979 -83.6127,29.5023 -83.5468,29.4067 -83.4369,29.2918 -83.3271,29.2343 -83.2172,29.1767 -83.1293,29.0616 -83.0195,29.0039 -82.9096,28.8886 -82.8437,28.7731 -82.8437,28.6767 -82.8876,28.5417 -82.8876,28.4258 -82.9316,28.3291 -82.9535,28.2711 -82.9535,28.1356 -82.9535,27.9998 -82.9755,27.9222 -82.9755,27.8639 -82.9755,27.6695 -82.9755,27.5722 -82.8656,27.4748 -82.8656,27.4163 -82.8437,27.221 -82.7558,27.0842 -82.602,26.8492 -82.4482,26.7315 -82.3603,26.5744 -82.2944,26.4564 -82.2284,26.4367 -82.1845,26.3974 -82.0966,26.358 -82.0966,26.2595 -82.0087,26.2201 -81.9867,26.0623 -81.8988,25.9636 -81.8769,25.8845 -81.7231,25.7659 -81.6352,25.6471 -81.5033,25.4091 -81.3495,25.2105 -81.3495,25.0514 -81.3495,24.9717 -81.2177,24.9717 -81.0859,24.892 -80.7782,24.9717 -80.4706,25.0912 -80.1191,25.4091 -80.0751,25.6867 -79.9872,26.0426 -79.9433,26.4367 -79.8993,26.7511 -79.8993,27.2992 -79.8993,27.5332 -79.8993,27.9222 -80.163,28.1162 -80.2509,28.2711 -80.3827,28.5031 -80.2948,28.696 -80.2948,29.2726 -80.7343,29.5406 -80.998,29.7697 -80.998,30.1505 -81.0419,30.3403 -81.0419,30.3687 -81.1078,30.5013 -81.1957,30.6148 -81.2946,30.7188 -81.3166,30.8038 -81.4814,30.8226 -81.7011,30.8792 -81.8,30.8792 -81.9318,30.8792 -82.0307,30.8226 -82.0966,30.7188 -82.1076,30.6054 -82.1076,30.6621 -82.3053,30.6621 -82.4591,30.6905 -82.6899,30.6905 -82.8986,30.6999 -83.2611,30.7188 -83.6347,30.7566 -84.0082,30.7754 -84.4147,30.8038 -84.7443,30.9264 -84.8432,31.0206 -84.8652,31.0582 -84.953,31.077 -85.1508,31.0958 -85.5463,31.0958 -85.8759,31.1147 -86.2934,31.1711 -86.7548,31.1335 -87.3261,31.1147 -87.5898,31.077 -87.6996,31.0206 -87.7216,30.9829 -87.7216,30.9264 -87.7436,30.9264 -87.7875))
[Habitat::initialize] Done.

[Graph::initialize]
Generate population grid and sample assignment
Loaded sample coordinates from ../data/gopher_eems.coord
There are 31 observed demes (out of 184 demes)
The population grid has 184 demes and 470 edges
There are 1191 samples assigned to 31 observed demes
[Graph::initialize] Done.

[Diffs::initialize]
Error: Check that the genotype matrix is a nIndiv-by-128 matrix.

Error regarding dissimilarity calculation on rEEMSplots.R

Dear Dr Petkova

I am running runeems_snps using the following parameters:
datapath = /home//sonia/Documents/pp
mcmcpath = /home/sonia/Documents/pp/test
nIndiv = 140
nSites = 1020
nDemes = 50
diploid = false
numMCMCIter = 10000000
numBurnIter = 8000000
numThinIter = 100000

The run finishes without warnings, but when I use the R script I get the following msg:

Plotting average dissimilarities within and between demes
/home/sonia/Documents/pp/
Error in if (ncol(x) != 2) { : argument is of length zero

When I look at my files rdist*, all of them seem incomplete (please see attached), and after doing some tests I could not figure out what might going on with my matrix. Could you please help me with that?
Thanks

test.zip

parameterization of negative-bi loglikelihood

Here in the code is the log-likelihood of the NegBi up to a constant (where prob, size) are treated as constants

*
 Log probability mass function of the negative binomial distribution, with
 parameters size (number of failures) and prob (probability of success)
 Also -- log pmf of the zero-truncated negative binomial, with the same parameters,
 as long as k>0, since truncating only changes the constant of proportionality
 Uses gamma(n) = (n-1)!
 */
double dnegbinln(const int k, const int size, const double prob) {
    double pln = -Inf;
    if ( (k>=0) && (size>0) && (prob>0) && (prob<1) ) {
        pln = lgamma(size+k) - lgamma(k+1) + k*log(prob);
    }
    return (pln);
}

However, shouldn't it be

pln = lgamma(size+k) - lgamma(k+1) + k*log(1-prob);

So I think the eems implementation has p as the probability of failure (not the probability of success).

Unknown loading error

Hi,

I used this package a number of months ago and everything worked perfectly. Now, however, I am receiving this error, and have no idea as to why the package won't load as it has in the past.

Any advise would be greatly appreciated.

** testing if installed package can be loaded from temporary location
Error: package or namespace load failed for ‘rEEMSplots’ in dyn.load(file, DLLpath = DLLpath, ...):
 unable to load shared object '/Users/joshuahallas/Library/R/4.0/library/00LOCK-rEEMSplots/00new/rEEMSplots/libs/rEEMSplots.so':
  dlopen(/Users/joshuahallas/Library/R/4.0/library/00LOCK-rEEMSplots/00new/rEEMSplots/libs/rEEMSplots.so, 6): Library not loaded: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libR.dylib
  Referenced from: /Users/joshuahallas/Library/R/4.0/library/00LOCK-rEEMSplots/00new/rEEMSplots/libs/rEEMSplots.so
  Reason: image not found
Error: loading failed
Execution halted
ERROR: loading failed
* removing ‘/Users/joshuahallas/Library/R/4.0/library/rEEMSplots’
Warning in install.packages :
  installation of package ‘rEEMSplots’ had non-zero exit status

error make linux command

Dear Dr. Petkova,

I got the following error message after updating the variables EIGEN_INC, BOOST_LIB, BOOST_INC in the Makefile.

/bin/ld: cannot find -lboost_system
/bin/ld: cannot find -lboost_program_options
/bin/ld: cannot find -lboost_filesystem
collect2: error: ld returned 1 exit status
make: *** [linux] Error 1

this is how it looks like in the Makefile

EIGEN_INC = /global/scratch/damasco/eems/eigen/
BOOST_LIB = /global/scratch/damasco/eems/boost_1_57_0/
BOOST_INC = /global/scratch/damasco/eems/boost_1_57_0/

I am not sure if I am doing it correctly. I saw on the boost program website that these libraries need to be built separately, but I have no clue how. Will it prevent to run eems properly?

thank you for your help
Gabriel

Problem compiling runeems_snps

Hi Dr. Petkova, I am running into a problem setting up runeems_snps due to its requirement of boost 1.57 or higher. I have an older ubuntu version that only allows up to boost 1.54, and setting up a newer ubuntu version would hardly be feasible for us.

This is the error report I get when trying to compile the program:

habitat.cpp: In member function ‘void Habitat::generate_outer(const string&)’:
habitat.cpp:27:8: error: ‘is_valid’ is not a member of ‘boost::geometry’
if (!boost::geometry::is_valid(domain)) {
^
make: *** [habitat.o] Error 1

I'm wondering if you have a work-around for this - e.g. a way to compile the runeems_snps program with boost 1.54?
Thanks for your help,
-Mike

Segmentation fault 11 Error

Hi Dr. Petkova,

When trying to run the sample files from the master package, I get Segmentation fault 11 error. I start with ./runeems_sats --params params-chain1.ini --seed 123, but the error shows up pretty much immediately. Here is the output:

calvisitor-10-105-166-192:src sasha$ ./runeems_sats --params params-chain1.ini --seed 123
Using Boost 1_57 (EEMS was tested with version 1_57)
Eigen 3.3.7 (EEMS was tested with version 3.2.2)

Segmentation fault: 11

Do you know why this might be?

Thanks a lot!

question about updates of hyperparameters

eems/runeems_snps/src/eems.cpp

Line 602 in 7c3e9c0

nowmrateS2 = draw.rinvgam(params.mrateShape_2 + 0.5 * nowmtiles, params.mrateScale_2 + 0.5 * SSm);

is this line and the line above it correct? The likelihoods here are truncated normal, not normal, which seems to mess up the usual conjugacy?

Grid input format

It might be useful at some point to define custom grid points outside the eems executable.

We were talking about grid input formats, I suppose the easiest would be to use the same format as is currently output in demes.txt and edges.txt, i.e.

demes.txt:

x and y coordiante of deme location, implicit indexing (n_demes rows, 2 cols)

       0        0
0.741935        0
 2.22581        0
 2.96774        0
 6.67742        0
 8.90323        0
 9.64516        0
  11.129        0
0.370968   0.7625
 1.85484   0.7625

edges.txt:

line i contains indeces of demes which are connected to i, (n_demes rows) max_n_deme cols (max_n_deme=6?); 0 implies no connection

 0   0   9   2   0   0
 1   9  80  72   0   0
72  10  11   4   0   0
 3  11  81  73   0   0
76  85  86  77   0   0
78  13  14   7   0   0
 6  14  87  79   0   0
79  15  16   0   0   0
 0  17  18  80   2   1
80  88  19  11   3  72

I guess the easiest way to implement this would be to reimplement

void Graph : : make _triangular_grid(...)

as the other two functions you described should be exactly the same as if the grid was constructed within eems.

Unable to plot dissimalirty plots

Hi, I have successfully created the main EEMS plot, however I an unable to plot the dissimilarity plots using -rdist in the shiny app.

When I load the -rdist.RData into R studio it shows that all the *.components are NULL.
How will I access the goof=dness of fit?
I combined 3 runs each from 200-1000 in step sizes of 100 demes each time. I have attached some files including the rdist and mrates01 from my run_eems_plot

test.zip.zip

full-rank

Plot EEMS with unique pch and col for each Population

Pipeline requires more libraries

Hello,
I am using EEMS in a server where I don't have sudo priviledges and it turns out EEMS needs more libraries than you say in the readme.
When I tried running it it also couldn't find fiona, collections, mpl_toolkit and descartes along with pandas and shapely.
Could you add those to the list of needed packages?

RAD distance matrix: choosing number of demes

Dear Desislava -

I am working on a analytical procedure for 2bRAD for non-model species where I try to avoid calling SNPs altogether and calculate between-individual genetic distances from the fraction of shared unique RAD tags. This seems to work quite well for PCA and hierarchical clustering, and I just made my first poke at EEMS (runeems_snps) and got some results! However I struggle to understand them. I only see color in my plots if I use low number of demes (25), if I go to 50 or 100 I just get zero across my whole seascape. This is of course not implausible - maybe there are really no barriers and no variation in genetic diversity in the Baltic sea for this species - but my 25-deme results seem to be reproducible from run to run.
I also see only three points in my rdist plots, this does not look quite normal...

cheers

Misha Matz, UT Austin
(attaching my polygon and populations - so far analyzed only the named ones)

error starting EEMS

Dear author,

I was trying to use runeems_sats and I have run into some issues. I am working with DArTseq data, using the silicoDArT which in essence is a dataset of microsatellites. In this specific case, I have 57326 loci for 158 samples. I exported all the files I need it to work with EEMS, and after loading them I get the following error (see below).

Thanks in advance,

[Habitat::initialize] Done.

[Graph::initialize]
Generate population grid and sample assignment
Loaded sample coordinates from ../data/qfly-nIndiv158-nSites57326.coord
There are 20 observed demes (out of 170 demes)
The population grid has 170 demes and 446 edges
There are 158 samples assigned to 20 observed demes
[Graph::initialize] Done.

[Diffs::initialize]
Read genotype data from ../data/qfly-nIndiv158-nSites57326.diffs
[Diffs::initialize] Done.

Initialize EEMS random state

[EEMS::initialize_state]
EEMS starts with 30 qtiles and 25 mtiles
[EEMS::initialize_state] Done.

Input parameters:
datapath = ../data/qfly-nIndiv158-nSites57326
mcmcpath = ../data/qfly-nIndiv158-nSites57326-EEMS-nDemes200-chain1
prevpath =
gridpath =
distance = euclidean
diploid = 0
nIndiv = 158
nSites = 57326
nDemes = 200
seed = 1518919729
numMCMCIter = 2000000
numBurnIter = 1000000
numThinIter = 9999
negBiSize = 10
negBiProb = 0.670000
qVoronoiPr = 0.250000
mrateShape = 0.000500
qrateShape = 0.002000
sigmaShape = 0.001000
qrateScale = 0.500000
mrateScale = 2.000000
sigmaScale = 1.000000
mSeedsProposalS2 = 0.010000
qSeedsProposalS2 = 0.100000
mEffctProposalS2 = 0.100000
qEffctProposalS2 = 0.001000
mrateMuProposalS2 = 0.010000

Initial log prior: -32923.86
Initial log llike: inf

Error while plotting map in eems.plots

Hello,

I am trying to add a map to my eems plot.

This is the code I ran:
library(rworldxtra)
library(rworldmap)
library(rgdal)
cmcpath='/eems/runeems_snps/run1'
plotpath='/eems/runeems_snps/run1'
coord__long_lat <- read.table(paste0('inputeems', ".coord"))
projection_none<-"+proj=longlat +datum=WGS84"
projection_mercator<-"+proj=merc +datum=WGS84"`

eems.plots(mcmcpath='/Users/nitinr/eems/runeems_snps/run1/', plotpath= '/Users/nitinr/eems/runeems_snps/run1/eems', longlat=TRUE, m.plot.xy = { points(coord__long_lat, col = "purple", pch=18) }, q.plot.xy = { points(coord__long_lat, col = "purple", pch=18) },out.png=FALSE),
projection.in = projection_none, projection.out = projection_mercator,
add.map =TRUE, col.map ="black", lwd.map =5)

This is the error I get:
Error in h(simpleError(msg, call)) : error in evaluating the argument 'x' in selecting a method for function 'nlayers': cannot do this transformation In addition: Warning message: In rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, : 701 projected point(s) not finite

How do I get around it? Thank you!

Unable to add points in rEEMSplots

Hello,

I'm using rEEMSplots to plot results, and I'd like to add the original points to the migration surface graph. I use the following syntax and do not get any errors, but the points are not appearing on my output graph.

plot_p <- rEEMSplots::eems.plots(mcmcpath = eems_results_P, plotpath = names_figures_P, add.grid = TRUE, projection.in = "+proj=longlat +datum=WGS84", projection.out = "+proj=longlat +datum=WGS84", add.outline = TRUE, longlat = FALSE, m.plot.xy = { points(coord_P, col = "purple", pch=18) }, q.plot.xy = { points(coord_P, col = "purple", pch=18) })

In addition, is it possible to overlay a streams/waterbodies shapefile on top of the migration/diversity plots?

Thank you!

Error starting EEMS

Hi Dr. Petkova,

I'm trying to run eems_sats. I've got it working with one set of data, so it seems to be running fine, but not with the other set. All my input files and input parameters read in properly but then I see:

Initial log prior: 97.3
Initial log llike: inf

[RunEEMS] Error starting EEMS.

I realize this might be difficult to disentangle without my files, but I've compared the inputs with those that have worked and do not see any differences. The genotype files look very similar, the size of the .sites matrix and .coord file matches the parameters, and the .outer is a simple rectangle. Is there perhaps something with the data itself that is causing this problem? Is 'initial log llike: inf' a problem? In runs that have worked for me this number has been a large negative number.

Any help would be greatly appreciated.

Many Thanks,
Emma

Error in grid.Call polygon edge not found

When I include the 'add_outline' argument in reemsplots2:

make_eems_plots(mcmcpath = "barrier-schemeZ-nIndiv300-nSites3000-EEMS-nDemes200-chain1", add_grid = TRUE, add_demes = TRUE, add_outline = TRUE)

I get the following error:

Error in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, :
polygon edge not found

Any idea what's causing this?

Thanks,
Tom

Error (rgeos):

Hello,

I am receiving the following error:

"The habitat geometry is not a valid ring (a ring is both simple and closed)
Let the habitat be the rectangle defined by (xmin, ymin) and (xmax, yxmax)
where xmin, xmax = range(longitude) and ymin, ymax = range(latitude).
Plotting posterior probability trace"

I have checked my outer.txt and the first and last lat long are the same. (attached)
outer.txt

Here is my script:

## Check that the current directory contains the rEEMSplots source directory (from GitHub)
if (file.exists("./rEEMSplots/")) {
  install.packages("rEEMSplots", repos=NULL, type="source")
} else {
  stop("Move to the directory that contains the rEEMSplots source to install the package.")
}


extdata_path <- ("/Users/joshuahallas/Documents/unr/projects/thamnophis_population/th_elegans/ELTE_subset/eems/rEEMSplots/results")

eems_results_elegans <- file.path(extdata_path, c("EL_25apr_Indiv97-nSites10121-EEMS-nDemes600-chain1", 
                                                  "EL_25apr_Indiv97-nSites10121-EEMS-nDemes600-chain2",
                                                  "EL_25apr_Indiv97-nSites10121-EEMS-nDemes600-chain3"))

eems_results_terrestris <- file.path(extdata_path, c("TE_29apr_Indiv90-nSites10121-EEMS-nDemes600-chain1", 
                                                  "TE_29apr_Indiv90-nSites10121-EEMS-nDemes600-chain2",
                                                  "TE_29apr_Indiv90-nSites10121-EEMS-nDemes600-chain3"))

name_figures <- file.path(path.expand("~"), "/Users/joshuahallas/Documents/unr/projects/thamnophis_population/th_elegans/ELTE_subset/eems/rEEMSplots/results")

library("rEEMSplots")
library("rworldmap")
library("rworldxtra")
library("rgdal")
library("Rcpp")
library("RcppEigen")
library("raster")
library("rgeos")
library("sp")

projection_none <- "+proj=longlat +datum=WGS84"
projection_mercator <- "+proj=merc +datum=WGS84"

eems.plots(mcmcpath = eems_results_elegans,
           plotpath = paste0(name_figures, "-elegans_geographic-map"),
           longlat = F,
           plot.width = 11, 
           plot.height = 11,
           projection.in = projection_none,
           projection.out = projection_mercator,
           add.map = TRUE,
           col.map = "black",
           lwd.map = 2,
           add.grid = FALSE,
           col.grid = "gray90",
           lwd.grid = 2,
           add.outline = TRUE,
           col.outline = "black",
           lwd.outline = 5,
           add.demes = TRUE,
           col.demes = "black",
           pch.demes = 5,
           min.cex.demes = 0.5,
           max.cex.demes = 1.5)

I have rerun older analyses and I received no error. I have also been able to run the current data using the function make_eems_plots from reemsplots2 and received no error. I am unable to figure out the issue. Any suggestions? Thank you for your time.

Error in *.coord

I am a new student in bioinformatics, so When I tried eems, some problems occurred that I could not solve well. I input the geographic coordinates of the samples into the *.coord file, but it reported an error: Check that *.coord is a list of locations, two coordinates per row. Maybe I just started so I didn’t find the problem in format description.

Some of my input coordinates are as follows:

15.90530 7.75420
15.90530 7.75420
15.90530 7.75420
15.90530 7.75420
15.90530 7.75420
8.56460 3.45940
8.56460 3.45940
8.56460 3.45940
8.56460 3.45940
8.56460 3.45940
8.56460 3.45940
8.56460 3.45940
4.85190 6.01280
4.85190 6.01280
4.85190 6.01280
4.85190 6.01280
4.85190 6.01280

please help me
thanks a lot

R_squared

Dear Dr Petkova:
How can I calculate the R2 statistics? I'm using reemsplots2 to plot my results.
Thank you!
Marina

shiny app error

Hello,

I am trying to use the shiny app. I uploaded the .RData file , however I got this error: An error has occurred. Check your logs or contact the app author for clarification.
I have no extra information about what is possibly wrong, I don't know if I am doing something wrong or if the app is not available at the moment? Any idea? Thank you!

Error: Check that the dissimilarity matrix is a square nIndiv-by-nIndiv matrix.

I am receiving this error message in the nDemes log file when trying to run EEMS, in the [Diffs::initialize] section.

"Error: Check that the dissimilarity matrix is a square nIndiv-by-nIndiv matrix."

I checked that the dissimilarity matrix (named datapath.diffs) is in fact a square, and checked that the end-of-line characters are in unix. The original structure file does have missing data (0s), but I do not think that should be a problem. There are no issues with the .outer and .coord files, and the program successfully initialized and completed the habitat and graph runs.

Thanks.

Posterior trace - nan in output

Hi there!

While I was using eems and plotting functions, i managed to get migration maps but no diagnostic figures are being generated. I got the following error:

Error in plot.window(...) : need finite 'ylim' values
In addition: Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf

I figured posteriors might be in mcmcpilogl.txt file and mine has only 'nan' in second column. What could cause such behaviour? Any help would be appreciated. In general, I couldn't find the annotation for output files in the manual. It would be super handy!

Thanks!
Rafal

Check that nDemes > 1

Hi Dr. Petkova,

I'm attempting to use eems_sats for a 21 loci microsatellite set with 357 individuals. I installed and compiled on a Mac (10.13.6) using Homebrew and Boost. I have prepped the input files as described in the manual and have tried specifying several different values for nDemes (20,50,100,200,500,1000) but no matter what number I specify I seem to get the error: "Check that nDemes >1."

Below are the contents of my params-chain1.ini modified file. My files are named datapath.coord, datapath.outer, and datapath.sites. I'm happy to email these files to you if they are helpful.

datapath = Documents/BlackBearMicroPopStruct/eemsinput/datapath
mcmcpath = Documents/BlackBearMicroPopStruct/eemsoutput/datapathEEMSchain1
nIndiv = 357
nSites = 21
nDemes = 100
diploid = true
numMCMCIter = 2000000
numBurnIter = 1000000
numThinIter = 9999

Any help that you can provide would be much appreciated!

Best,

Lauren

Question about the nSites parameter

Dear Desislava,
I would like to ask what the "nSites" parameter will be used for?

Based on the above formula in the instruction Manual(Formula 1 in the "Observed vs fitted dissimilarities" section), I assume that the genetic dissimilarity matrix was already normalised by the number of Sites.
In this case, why you would require an input of nSites in the *.ini files?
My genetic dissimilarity matrix was computed using 3 million common SNPs since it is whole genome low-pass sequencing data. I am a bit afraid that such a large nSites will result in erroneous modelling in the default setting.
Could you please advise?

Error: The population grid is not connected.

Hello,
I'm trying to use runeems_snp for some data and I get the above error saying the population grid is not connected.
I did a search online after looking at the source code and I can't figure out what could cause that error or how I can fix it.

My polygon coordinates start and end at the same point on the map(the locations for the polygon vertices were found by looking up the appropriate coordinates on google maps and are tab-delimited geographical coordinates. The ones for the samples were in the data. )

here's the stderror and eemsrun.txt file:

Using Boost 1_53 (EEMS was tested with version 1_57)
Eigen 3.2.2 (EEMS was tested with version 3.2.2)

Will save EEMS output to ./0-0.output/
[Habitat::initialize]
Loaded habitat points from ./0-0/0-0.outer
[Habitat::initialize] Done.

[Graph::initialize]
Generate population grid and sample assignment
Loaded sample coordinates from ./0-0/0-0.coord
Error: The population grid is not connected.

parameters:

Input parameter values:
datapath = ./0-0/0-0
mcmcpath = ./0-0.output
prevpath =
gridpath =
distance = euclidean
diploid = 0
nIndiv = 2582
nSites = 1233553
nDemes = 200
seed = 1519055359
numMCMCIter = 2000000
numBurnIter = 1000000
numThinIter = 9999
negBiSize = 10
negBiProb = 0.670000
qVoronoiPr = 0.250000
mrateShape = 0.000500
qrateShape = 0.002000
sigmaShape = 0.001000
qrateScale = 0.500000
mrateScale = 2.000000
sigmaScale = 1.000000
mSeedsProposalS2 = 0.010000
qSeedsProposalS2 = 0.100000
mEffctProposalS2 = 0.100000
qEffctProposalS2 = 0.001000
mrateMuProposalS2 = 0.010000

Question about adding geographic map

Hi Desislava,
thank you very much for your detailed introduction of EEMS. I have a question about adding geographic map. I can add borders to countries. Can boundaries be added within countries?For example, add the provinces of China.
Looking forward to your reply.
Best,
Hanna

Issue with installation

Hi Dr. Petkova,

I'm having trouble with installing bed2diffs properly. I've downloaded the EEMS files and the plinkio files from Github (I put both master directories in one folder labeled "EEMS"). Then I edited the PLINKIO variable in the makefile in "src-wout-openmp" (I'm working on a Mac) to where I've downloaded the "libplinkio-master" directory. However, I'm getting this error when I run make darwin in terminal within the src-wout-openmp directory:

Nehas-MacBook-Pro:src-wout-openmp nehasavant$ make darwin
c++ -I/Users/nehasavant/Dropbox/Columbia_New_York/Salamanders/Downstream_Analyses/EEMS/libplinkio-master//include -O3 -Wall -Werror -c bed2diffs_v1.cpp -o bed2diffs_v1.o
In file included from bed2diffs_v1.cpp:3:
./data.hpp:7:10: fatal error: 'plinkio/plinkio.h' file not found
#include <plinkio/plinkio.h>
         ^~~~~~~~~~~~~~~~~~~
1 error generated.
make: *** [bed2diffs_v1.o] Error 1

I'm not quite sure what I'm doing wrong. Any thoughts or advice would be much appreciated!

Thank you,
Neha

Problem running bed2diffs: Error opening plink files

Hi Dr. Petkova, I get an error message when trying to run bed2diffs on my data. Here is the command and error message:

crossley@denali:~/EEMS_analysis/eems-master/bed2diffs/src$ ./bed2diffs_v1 --bfile b33
Set the number of OpenMP threads to 1

Compute the average genetic differences according to:
Dij = (1/|Mij|) sum_{m in Mij} (z_{im} - z_{jm})^2
where Mij is the set of SNPs where both i and j are called

[Data::getsize] Error opening plink files b33.[bed/bim/fam]

I am running EEMS from an Ubuntu terminal, which has bed2diffs compiled/configured. I generated my BED file from PED/MAP files using PLINK. I understand you may need more information from me before you can diagnose the problem, but wondered if it is already apparent to you what the issues are.

Thank you very much for any insight you can give,
-Mike

Error in bitmap

Dear Dr. Petkova,

I am trying to use eems.plots but I keep getting this error message :
Error in bitmap(paste0(plotpath, "%02d.png"), type = "png16m", res = plot.params$res, :
GhostScript was not found
I have tried running the script with the example data using the same input commands from the EEMS doc. Any help on the matter is greatly appreciated.

Best,
Zach

Error when attempting to plot runeems_snps results

Hi Dr. Petkova, I came across an error trying to plot my EEMS results in R that I'm not sure how to resolve. Here is some of my last output:
Processing the following EEMS output directories :
./run1
Plotting effective migration surface (posterior mean of m rates)
./run1
Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false, :
TopologyException: Input geom 1 is invalid: Self-intersection at or near point -90.404795806411528 47.0755877112326 at -90.404795806411528 47.0755877112326

I looked into this a bit online, and some have suggested the issue arises when dealing with unprojected coordinates. But I'm not sure how to specify a projection in the data I'm plotting. I tried your code for specifying a projection (page 15 in the EEMS manual), but got the same error.
Do you have any insight on how to resolve this error?

-Mike

Modifying maps from EEMS

Is it possible to modify the maps in any way-- and how I would go about doing that? For example, I would like to produce a map with the EEMS surface and my sample points on it. I know how to produce such a map in R with a plain map, but the EEMS program just outputs a PDF of the maps, and therefore they aren't modifiable as they are. Would I need to go in and edit the R code for the program to be able to do this? Or is there an easier way?

Thanks!

Population grid problem - ring is not simple or complete - change of coordinate system and horizontal line due to split polygons solution

Hello,

I previously had a problem where I was trying to use a polygon that incorporated samples from the Chukchi peninsula where their coordinates are negative with others from Siberia having positive coordinates. (For example 179.0 and then -179.0 ) which created problems and I ended up not using that area, but now I must use it.

One way I have resolved this in the past is to add 360 to those negative coordinates(for the environment and the individual coordinates) and have 181 for the second one for example. But in EEMS this gives an error saying the “ring is not simple or complete”.

Do you have a way of including such individuals into EEMS from across Siberia? How can I include those individuals when the coordinates make a quick change from 179 to -179 degrees? Is something like that even possible in EEMS?

Dingbats on pdf output of rEEMS plot

Just a small issue: When saving files as pdf, R by default renders circles as characters of the Dingbats font to reduce size. This generates problems when editing pdfs with some vector image programs (e.g. Inkscape and Illustrator), since the dots don't always show up. An easy fix is adding "useDingbats=F" in the pdf command. In the package code this fixed it:

save.graphics <- function(plotpath, plot.params) {
if (plot.params$out.png) {
bitmap(paste0(plotpath, '%02d.png'), type = 'png16m', res = plot.params$res,
height = plot.params$height, width = plot.params$width, units = 'in')
} else {
pdf(paste0(plotpath, '%02d.pdf'),
height = plot.params$height, width = plot.params$width, onefile = FALSE, useDingbats=F)
}
}

Perhaps future releases could have an argument to chose whether using Dingbats in the eems.plots function?

Thanks!

Roberto