Using-divMigrate

Email: kkeenan02 -at- qub.ac.uk

Introduction

Since diveRsity v1.9.5, a new version of the function divMigrate is available. This new function implements C++ code making it up to 20X faster than the old version. The new version also allows users to test the significance of directional differentiation using a 95% confidence interval rejection method. An integrated bootstrapping method is used to estimate the confidence region of relative migration.

Downloading diveRsity

To get the latest version of diveRsity from github, paste the following into the R console:

if ("devtools" %in% rownames(installed.packages()) == FALSE) {
    install.packages("devtools", repo = "http://cran.rstudio.com", dep = TRUE)
}
library("devtools")
install_github("diveRsity", "kkeenan02")

This code will download and install the package devtools (if required), which is then used to install diveRsity from github.

Example

The following example demonstrates how divMigrate might be used to visually explore the direction of differentiation within a given system. The data used in this example have been simulated using fastsimcoal2. The input parameter file is as follows:

//Number of population samples (demes)
2 samples to simulate :
//Population effective sizes (number of genes)
1000
1000
//Samples sizes
100
100
//Growth rates  : negative growth implies population expansion
0
0
//Number of migration matrices : 0 implies no migration between demes
1
//migration matrix
0.000 0.005
0.000 0.000
//historical event: time, source, sink, migrants, new size, new growth rate, migr. matrix 
0  historical event
//Number of independent loci [chromosome] 
1 0
//Per chromosome: Number of linkage blocks
1
//per Block: data type, num loci, rec. rate and mut rate + optional parameters
MICROSAT  10  0.0000  0.0005 0 0

The important information to notice from this file is the migration matrix. This simulates migration from pop1 to pop2 at a rate of 0.005 genes per generation. Since fastsimcoal2 is a coalescent simulator, this is equivalent of migration from pop2 to pop1 forward in time.

fastsimcoal2 generates an arlequin genotype file, which can be converted to a genepop file using the arp2gen function in diveRsity, as follows:

library(diveRsity)
arp2gen("2popSTRmigr.arp")

This will generate a genepop file, named 2popSTRmigr.gen in the current working directory.

Following conversion, we can now read the genepop file and run divMigrate.

run1 <- divMigrate("2popSTRmigr.gen")

Caution! The method used in this function is still under development. 

This command will simply calculate the relative migration according to Sundvqist et al, and plot the resulting network. We can see that the method successfully resolves the migration relationship simulated. The raw number used to draw the network are stored in the results object run1.

The new version of divMigrate also provides a filter mechanism for noisy migration networks. For example, if we were to run divMigrate on Test_data (distributed with diveRsity), the resulting network is complex.

tst <- divMigrate(Test_data)

## Caution! The method used in this function is still under development.

As we can see in the network, many of the migration connections are defined by very weak relative migration. By setting an appropriate filter_threshold we can tidy the network up a bit.

tst2 <- divMigrate(Test_data, filter_threshold = 0.5)

## Caution! The method used in this function is still under development.

Now we see that only relative migration above 0.5 is represented in the network.

The new divMigrate function also provided a method to test for the significance of directional migration between all pairs of population samples. Going back to the original data. This would be done as follows:

system.time({
run2 <- divMigrate("2popSTRmigr.gen", nbs = 1000)
})

## Caution! The method used in this function is still under development.

##    user  system elapsed 
##    3.51    0.00    3.57

Here we see that both the original network is plotted, as well as a network containing only edges for significant migration (i.e. non-overlapping relative migration 95% CI).

We can also see that 1000 bootstrap iterations for these data only takes 3.6 sec. For larger data, an option to use parallel processing is also provided. This can be activated by specifying para = TRUE in the function command.

The colour of network edges can also be controlled using the plot_col argument:

run3 <- divMigrate("2popSTRmigr.gen", nbs = 1000, plot_col = "red")

## Caution! The method used in this function is still under development.

Additional information about the divMigrate function can be found by typing the following into the R console:

?divMigrate