Defining populations in PopGenome after running readMS
1
1
Entering edit mode
7.1 years ago
ncrouc2 ▴ 20

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

R PopGenome ms • 3.2k views
ADD COMMENT
0
Entering edit mode

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?

ADD REPLY
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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?

ADD REPLY
0
Entering edit mode

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

ADD REPLY
1
Entering edit mode
7.1 years ago
ncrouc2 ▴ 20

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))

ADD COMMENT
0
Entering edit mode

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.

ADD REPLY

Login before adding your answer.

Traffic: 1696 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6