Defining populations in PopGenome after running readMS

Question

[ncrouc2]

Hi,

I am running simulations in ms and then (wanting to) analyze the output using PopGenome in R, however I am having trouble defining the different populations when I load the data.

For example, say I run the following command line: ./ms 6 1 -t 6 -I 3 2 2 2 -ej 1 3 2 -ej 2 2 1 > test_output

I can load the output into R:

library(PopGenome)

x <- readMS(test_output)

I then run into trouble trying to define the fact that there are 3 populations. For example I tried (following the tutorial) : pop.1 <- c("seq1", "seq2")

pop.2 <- c("seq3", "seq4")

pop.3 <- c("seq5", "seq6")

x <- set.populations(x, list(pop.1, pop.2, pop.3))

I think that this is probably due to how the sequences are named once they are loaded into PopGenome, but I haven't been able to figure it out. Any pointers would be greatly appreciated

Thanks

[ddiez]

Are you getting some error when running set.populations()? Why do you say you are having trouble? Not familiar with PopGenome package but, what does x look like?

[ncrouc2]

Running the code as above doesn't throw errors, but when I try to perform additional steps, e.g. neutrality.stats(x), then I get Error in populations[[xx]] : subscript out of bounds. x here is class GENOME and has 110 slots corresponding to different statistics about the data.

[ddiez]

My only guess is that you might be right about the ames of the sequences. How can you know the names? What information get.details() gives you? What I am thinking is some method that will tell you the correct names. Also, can you run a minimal example (i.e. not with your own data; maybe there are some in the documentation) without errors?

[ncrouc2]

So far I have just followed the example in the documentation and tinkered without success

Answer 1

[ncrouc2]

I have found the solution. The get.individuals function provides a list of names of the sequences. In this instance, it was "1", "2", "3", "4", "5", "6". To get the example above to work is, therefore:

pop.1 <- c("1", "2")

pop.2 <- c("3", "4")

pop.3 <- c("5", "6")

x <- set.populations(x, list(pop.1, pop.2, pop.3))

[ddiez]

Great you found a solution! One comment: here and in your question, use the 101010 button to format chunks of code so that they look better. Take a look to other questions/answers to see the difference.