PopGenome package for R, problem calculating neutrality stats

Question

[aga.lipinska]

Dear all,

I have been using PopGenome with VCF files to get some neutrality statistics for RAD-Seq data. I used mpileup and bcftools to get the variant calls and have individual VCF files per chromosome with genotype data for 3 populations and 1 outgroup. Everything seems to work, but I am getting these strange values for Tajima's D above 5 and more. These values only appear when I use the option set.populations:

get.neutrality(data)[[1]]     #(population 1)
     Tajima.D n.segregating.sites Rozas.R_2  Fu.Li.F  Fu.Li.D Fu.F_S Fay.Wu.H    Zeng.E Strobeck.S
test 9.555944                 183        NA 5.185575 1.276137     NA 4.054317 -6.170585         NA

If I treat my samples as one populations I get:

get.neutrality(data)[[1]]     #(all populations)
     Tajima.D n.segregating.sites Rozas.R_2  Fu.Li.F  Fu.Li.D Fu.F_S Fay.Wu.H Zeng.E Strobeck.S
test 1.353448                1018        NA 1.400123 0.999412     NA      NaN    NaN         NA

Below are the commands I use:

data <- readData(path="./vcf", format="VCF", include.unknown=TRUE, gffpath="./gff")
data <- set.populations(data, list(c("ind1", "ind2"..), c("ind3", "ind4"..), c("ind5", "ind6"..), c("ind7", "ind8"..)),  diploid=TRUE)
data <- set.outgroup(data, c("ind1", "ind2"..), diploid=TRUE)
data <- neutrality.stats(data)
get.neutrality(data)[[1]]

What could be the problem? If anyone has any suggestions, I would be very grateful!