search

grunwaldlab / poppr

🌶 An R package for genetic analysis of populations with mixed (clonal/sexual) reproduction
https://grunwaldlab.github.io/poppr
68 stars 26 forks source link

Optional disabling of statistics in poppr table #248

Open Neato-Nick opened 3 years ago

Neato-Nick commented 3 years ago

Please place an "x" in all the boxes that apply

Would it be possible to enable or disable some of the diversity statistics calculated in the poppr table? In a population I'm analyzing, .ia is getting an error because of negative values. It's already been shown that this population is very highly clonal, but I am interested in the other diversity indices aside from index of association. A workaround I see for this is to just not run the index of association, maybe this looks like a parameter that takes a character vector of 'calculations to run' and by default it just runs all of them?

zkamvar commented 3 years ago

In a population I'm analyzing, .ia is getting an error because of negative values.

That's.... odd. This function shouldn't be erroring on negative values of IA. Can you include a copy of this error?

Would it be possible to enable or disable some of the diversity statistics calculated in the poppr table?

I might be able to implement that. The poppr() function is one of those weird legacy functions whose scope changed over time.

FWIW, You can get most of the way to the table with the diversity_stats() function, which gives you clonality statistics from the MLG counts and I think I could write in the rarefaction in there (but I'd need to check). As for stats like Hexp, that one I really should just export since it calculates the corrected version (whereas summary() calculates the uncorrected version).

If all you need are MLG stats, then diversity_table() is your best bet.

Neato-Nick commented 3 years ago

That's.... odd. This function shouldn't be erroring on negative values of IA. Can you include a copy of this error?

I don't have a reprex for the data that caused the error, but this is the error I'm getting. I've gotten negative ia results before and poppr reported them, I'm guessing this is maybe a combination of large data (>=130 individuals w/40k snps) and very high clonality? Negative length vector is a weird error and my google-fu didn't get me very far. I kind of was just thinking this was due to my very highly clonal data. Error in .Ia.Rd(popx, missing) : negative length vectors are not allowed

diversity_stats()

Thanks, this is helpful! Clonality stats are good. I do want Hexp eventually too though. Looks like it doesn't split out populations nicely like the poppr table does though, guess I'd need to seppop and run mlg.table() and diversity_stats() for each one?

zkamvar commented 3 years ago

Error in .Ia.Rd(popx, missing) : negative length vectors are not allowed

This is definitely a bug and I'm not sure where it's coming from. Would you be willing to send me a copy of your data so I could check it out?

Neato-Nick commented 3 years ago

Hi @zkamvar , sorry I ghosted you. I was tracking down other potential issues in my data to make sure it was not a bug upstream and I feel much more confident in my data now. I'm e-mailing you at the address in your bio.

zkamvar commented 3 years ago

I realize I responded via email and not on GitHub.

This is happening because there are too many loci being fed into the index of association function. It needs to construct a vector of n(n-1)/2 (where n is the number of loci) to calculate the denominator and the number of loci presented is just too large. The recommended alternative is to convert your source data to genlight or snpclone and use something like samp.ia() to calculate the index off of a random sample of loci.

Details for nerds

The reason this is happening is because there is an integer overflow in the construction of the covariance matrix (the count variable on line 111):

https://github.com/grunwaldlab/poppr/blob/94c745aa2e6a806a125f48b86f9c07d8b6ec95aa/src/poppr_distance.c#L95-L116

I can avoid this by specifying a long int (or whatever solution I came up with to solve #234), adding a warning about too many loci, and adding a meaningful error when the counter overflows.

Neato-Nick commented 3 years ago

Ah, that totally makes sense. It worked for ~22k loci but crashed at ~56k loci, but I recognize that even ~22k loci is a lot to use. I'll definitely subset my data.

Thanks for digging into the weeds to figure out what's going on. I must be a huge nerd, 234 and #235 were interesting!