Question: PopGenome package for R, problem calculating neutrality stats
gravatar for aga.lipinska
3.6 years ago by
aga.lipinska0 wrote:

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)



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


snp R • 1.7k views
ADD COMMENTlink written 3.6 years ago by aga.lipinska0
Please log in to add an answer.


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