Question: Defining populations in PopGenome after running readMS
1
gravatar for ncrouc2
2.2 years ago by
ncrouc220
ncrouc220 wrote:

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

ms popgenome R • 898 views
ADD COMMENTlink modified 2.2 years ago • written 2.2 years ago by ncrouc220

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 REPLYlink written 2.2 years ago by ddiez1.7k

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 REPLYlink written 2.2 years ago by ncrouc220

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 REPLYlink written 2.2 years ago by ddiez1.7k

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

ADD REPLYlink written 2.2 years ago by ncrouc220
1
gravatar for ncrouc2
2.2 years ago by
ncrouc220
ncrouc220 wrote:

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 COMMENTlink written 2.2 years ago by ncrouc220

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 REPLYlink written 2.2 years ago by ddiez1.7k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1013 users visited in the last hour