This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Popgenome nvalid sites

Hi all,

I would like to use the r program popgenome to calculate some neutrality stats. However, everytime that I read the data I have 0 or NaN valid sites.

For example, if I read all chromosome 1 like this:

GENOME.class_all <- readVCF("File.vcf.gz", numcols=100000, tid="Chr01", frompos=1, topos=159207788, include.unknown=T)

I have this:

get.sum.data(GENOME.class_all)

                           n.sites     n.biallelic.sites n.gaps n.unknowns n.valid.sites
1 - 159207788            15920778                 8542        0          0             0
                       n.polyallelic.sites trans.transv.ratio
1 - 159207788                           0                 NaN

Playing with the numcols, frompos, and topos parameters sometimes I have n.valid.sites 0 or NaN but apparently popgenome is not recognizing any valid site.

My vcf file was imputed with beagle5 and for Chr01 I have 91955 sites for 25 individuals, where the last SNP has the position 159207788. I am wondering if I am doing something wrong.

vcf popgenome
• 1,717 views • link
written by shinken123 shinken123

1 answer

The reason why some slots are empty it is because the readVCF function only fills the next slots:

Slot Description

  1. n.sites total number of sites
  2. n.biallelic.sites number of biallelic sites
  3. region.data some detailed information about the data read
  4. region.names names of regions

This info is here just in case someone else has doubts.

The function readData it include n.valid.sites

Slot Description

  1. n.sites total number of sites
  2. n.biallelic.sites number of biallelic sites
  3. n.gaps number of sites with gaps
  4. n.unknowns number of sites with unknown nucleotides
  5. n.valid.sites number of valid sites
  6. n.polyallelic.sites number of sites with >2 nucleotides
  7. trans.transv.ratio transition/transversion ratio of biallelic sites
  8. region.names names of regions
  9. region.data some detailed information about the data read

So with this funtion I could obtain that value.

• 0 views • link
written by shinken123 shinken123

Log in to answer this question.

More posts like this