Forum:PopGenome nucleotide diversity

Question

[yp]

Hi Group,

I would like to calculate theta watterson using genome-wide SNPs. I am using PopGenome software. I find this software is very useful for my study, but very difficult to work with.

I have 10291 SNPs distributed on 17 chromosomes. How do I upload?

I could upload my data in the VCF format, but unable to calculate Theta Watterson. I would like to use sliding window approach to perform nucleotide diversity.

Could you please suggest me suitable code? In addition, could you please send me a link to download Arabidopsis sub-directory and SNP input format?

Here are the codes I used:

GENOME.class<-readVCF("Bc_test.vcf.gz",tid="1",from=1,to=10291,approx=FALSE,gffpath=FALSE,numcols=10291)  # for tid=" 1:17", I want to use SNP from all the 17 chromosomes

get.neutrality((GENOME.class,theta=TRUE,stats=TRUE))
get.neutrality(GENOME.class)[[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
1 - 10291 0.7972621                   1        NA 0.6281656 0.3860124    NaN      NaN    NaN        NaN

#Sliding widow
GENOME.class.slide<-sliding.window.transform(GENOME.class,width=10000,jump=10000)
|            :            |            :            | 100 %
|

Error in SLIDE.POS[[zz]] <- window + snp.begin : 
  attempt to select less than one element

Answer 1

[Kolea Zimmerman]

To get Watterson's theta, you first need to calculate neutrality stats using the neutrality.stats() function, then use the get.neutrality() function to return the stats. See the documentation here: http://www.inside-r.org/node/132345

I think you need to provide the readVCF function with a GFF file if you are going to calculate sliding window stats based on reference genome position. As it is now, there is no information provided to the function that indicates the relative positions within the genome of the SNPs in the VCF file, so it is calculating window stats based on your SNPs (almost all of your SNPs would be in 1 window). To calculate genome based sliding windows, you would indicate this when you call the function by specifying: type = 2 (documentation: http://www.inside-r.org/node/132345).

I recommend this tutorial so you can familiarize yourself with how the package works: https://cran.r-project.org/web/packages/PopGenome/vignettes/Whole_genome_analyses_using_VCF_files.pdf

[yp]

Thanks for you help. How can I get theta Watterson per site?

[Kolea Zimmerman]

I'm not sure I understand your question. If you are asking how to calculate a different Watterson's Estimator for each polymorphic site, that is not really reccomended because watterson's estimator is an estimate of the mutation rate for an average site in a genome/gene, the estimate does not differ site by site. The calculation requires multiple polymorphic sites.

If you are asking how to just return the Watterson's theta values for your genes/sliding windows, then just use the get.neutrality() function to return the stats.

Hope that helps!