jgx65 / hierfstat

the hierfstat package
24 stars 14 forks source link

allelic.richness rounded numbers #58

Closed Gillous closed 2 years ago

Gillous commented 2 years ago

I am contacting you since I am encountering an issue with the Package : Hierfstat ; and the Fonction : allelic.richness My guess is that my input file contains too many entries, as if cut down to ~200 entries, the function is working great again.

We seek to calculate with "allelic.richness" an index of allelic richness diversity by rarefaction, for 4 populations, called DAPC 1 to 4. We obtain :

          1  2 3         4

VNTR02 6.284283 9 5 8.216615 VNTR03 1.000000 7 6 2.976744 VNTR04 4.506138 16 2 15.008352 VNTR06 1.777778 4 1 2.486105 VNTR07 1.777778 11 9 5.642302 VNTR08 3.506631 8 4 6.959234 VNTR09 1.000000 3 1 1.488372 VNTR10 8.717586 18 4 17.260108 VNTR12 1.000000 6 4 4.605293 VNTR13 1.000000 2 2 2.000000 VNTR17 1.951076 6 3 5.483883 VNTR18 2.777285 6 3 3.976742 VNTR19 1.951076 3 2 1.932133

However, what troubles us is that population 2 only shows whole numbers; which is normal for population 3 which has the fewest individuals. Thus, it is believed that the results of this population 2 are not good.

Do you have time to look at this issue? Should I post it on github or somewhere else?

Input file is in enclosure and script is here :

library(hierfstat) vn13 = read.csv("Hierfstat_inputfile.csv",sep=";",h=T, quote="") allelic.richness(vn13, diploid=FALSE)

Best regards,

Gilles Hierfstat_inputfile.csv

jgx65 commented 2 years ago

Hi @Gillous , there is indeed an issue, an overflow in the calculation that is not caught. I will fix it, but meanwhile, if you enter

allelic.richness(vn13,diploid=FALSE,min.n=250) You'll get something more meaningfull:

> allelic.richness(dat,diploid=FALSE,min.n=250)
$min.all
[1] 250

$Ar
              1        2        3         4
VNTR02 5.868205 1.357910 4.956064  8.041457
VNTR03 1.000000 1.178955 5.678717  2.830565
VNTR04 4.204454 2.279455 2.000000 14.719167
VNTR06 1.661376 1.325990 1.000000  2.407629
VNTR07 1.661376 1.000000 8.999545  5.316090
VNTR08 3.208679 1.886684 4.000000  6.791307
VNTR09 1.000000 1.000000 1.000000  1.415282
VNTR10 8.503459 6.807752 3.999555 17.029658
VNTR12 1.000000 2.352424 3.978036  4.458817
VNTR13 1.000000 1.325990 1.850340  2.000000
VNTR17 1.885928 2.253210 2.978036  5.402067
VNTR18 2.657151 1.749308 2.999992  3.830534
VNTR19 1.885928 1.325990 1.996841  1.883937

Cheers

jgx65 commented 2 years ago

just uploaded a new version of hierfstat with the issue fixed (commit #97987d3), let me know if this works and if so, please close the issue.

Gillous commented 2 years ago

Wonderful! Everything seems in order!

Great work.

Gilles