Hi Adam,
Please can you explain how the window-pi calculation works, and what it means intuitively? I dare say there is a simple formula but I can't find it anywhere. I thought it was a simple average of the pi values for all variable sites within that window, but clearly it isn't, and I can't see how it's incorporating invariant sites.
Sorry for the confusion!
Cheers,
Hilary
Not really a bug-report... :-)
Pi is also known as nucleotide diversity, and is the estimate of the average number of differences between a pair of chromosomes. In a window, there will be lots of sites where the chromosomes match, and hence you need to account for those sites in the calculation. If useful, you can inspect the source code for the calculation. It's in the function 'output_windowed_nucleotide_diversity' in vcf_file_output.cpp .
OK, thanks. I've had a look at the code now and understand what it's doing, but I still have a problem. I was trying to calculate pi for each chromosome by setting the window size equal to the chromosome length. I was getting sensible numbers except for the largest chromosome, for which vcftools told me that pi=7.28949 . This is clearly wrong, and when I calculate it in R using the pi values from the sites.pi file and the N_CHR values from the corresponding frq file, I get something sensible. I think C++ is having trouble for this chromosome because the numbers are too big and the relevant variables are unsigned ints. Specifically, at line 2938:
N_pairs = bins[CHROM][s][N_variant_site_pairs] + (N_monomorphic_sites * N_comparisons);
For my largest chromosome,
total.n.variant.site.pairs = 26527590
N_comparisons = 72
N_monomorphic_sites = 118949600
Thus, N_pairs ~8.6 billion, but it was declared as unsigned int at line 2922 . Is this the problem? If so, is it possible to declare the variables differently so one can calculate pi for a large window (e.g. a whole chromosome)?
Many thanks,
Hilary
From: Adam Auton auton1@users.sf.netamp#97;amp#117;amp#116;amp#111;amp#110;amp#49;amp#64;amp#117;amp#115;amp#101;amp#114;amp#115;amp#46;amp#115;amp#102;amp#46;amp#110;amp#101;amp#116;
Reply-To: "[vcftools:bugs]" 42@bugs.vcftools.p.re.sf.netamp#52;amp#50;amp#64;amp#98;amp#117;amp#103;amp#115;amp#46;amp#118;amp#99;amp#102;amp#116;amp#111;amp#111;amp#108;amp#115;amp#46;amp#112;amp#46;amp#114;amp#101;amp#46;amp#115;amp#102;amp#46;amp#110;amp#101;amp#116;
Date: Mon, 25 Feb 2013 20:17:52 +0000
To: "[vcftools:bugs]" 42@bugs.vcftools.p.re.sf.netamp#52;amp#50;amp#64;amp#98;amp#117;amp#103;amp#115;amp#46;amp#118;amp#99;amp#102;amp#116;amp#111;amp#111;amp#108;amp#115;amp#46;amp#112;amp#46;amp#114;amp#101;amp#46;amp#115;amp#102;amp#46;amp#110;amp#101;amp#116;
Subject: [vcftools:bugs] #42 window-pi calculation
Not really a bug-report... :-)
Pi is also known as nucleotide diversity, and is the estimate of the average number of differences between a pair of chromosomes. In a window, there will be lots of sites where the chromosomes match, and hence you need to account for those sites in the calculation. If useful, you can inspect the source code for the calculation. It's in the function 'output_windowed_nucleotide_diversity' in vcf_file_output.cpp .
[bugs:#42]http://sourceforge.net/p/vcftools/bugs/42/ window-pi calculation
Status: closed
Created: Mon Feb 25, 2013 06:14 PM UTC by Hilary
Last Updated: Mon Feb 25, 2013 06:14 PM UTC
Owner: nobody
Hi Adam,
Please can you explain how the window-pi calculation works, and what it means intuitively? I dare say there is a simple formula but I can't find it anywhere. I thought it was a simple average of the pi values for all variable sites within that window, but clearly it isn't, and I can't see how it's incorporating invariant sites.
Sorry for the confusion!
Cheers,
Hilary
Sent from sourceforge.net because you indicated interest in https://sourceforge.net/p/vcftools/bugs/42/
To unsubscribe from further messages, please visit https://sourceforge.net/auth/prefs/
Related
Bugs:
#42I've made a change to the version in the SVN. Please check out and see if this fixes your issue.
Thanks,
Adam
Yes, seems to be OK now – thank you!
Hilary
From: Adam Auton auton1@users.sf.netamp#97;amp#117;amp#116;amp#111;amp#110;amp#49;amp#64;amp#117;amp#115;amp#101;amp#114;amp#115;amp#46;amp#115;amp#102;amp#46;amp#110;amp#101;amp#116;
Reply-To: "[vcftools:bugs]" 42@bugs.vcftools.p.re.sf.netamp#52;amp#50;amp#64;amp#98;amp#117;amp#103;amp#115;amp#46;amp#118;amp#99;amp#102;amp#116;amp#111;amp#111;amp#108;amp#115;amp#46;amp#112;amp#46;amp#114;amp#101;amp#46;amp#115;amp#102;amp#46;amp#110;amp#101;amp#116;
Date: Tue, 26 Feb 2013 15:16:02 +0000
To: "[vcftools:bugs]" 42@bugs.vcftools.p.re.sf.netamp#52;amp#50;amp#64;amp#98;amp#117;amp#103;amp#115;amp#46;amp#118;amp#99;amp#102;amp#116;amp#111;amp#111;amp#108;amp#115;amp#46;amp#112;amp#46;amp#114;amp#101;amp#46;amp#115;amp#102;amp#46;amp#110;amp#101;amp#116;
Subject: [vcftools:bugs] #42 window-pi calculation
I've made a change to the version in the SVN. Please check out and see if this fixes your issue.
Thanks,
Adam
[bugs:#42]http://sourceforge.net/p/vcftools/bugs/42/ window-pi calculation
Status: closed
Created: Mon Feb 25, 2013 06:14 PM UTC by Hilary
Last Updated: Mon Feb 25, 2013 08:17 PM UTC
Owner: nobody
Hi Adam,
Please can you explain how the window-pi calculation works, and what it means intuitively? I dare say there is a simple formula but I can't find it anywhere. I thought it was a simple average of the pi values for all variable sites within that window, but clearly it isn't, and I can't see how it's incorporating invariant sites.
Sorry for the confusion!
Cheers,
Hilary
Sent from sourceforge.net because you indicated interest in https://sourceforge.net/p/vcftools/bugs/42/
To unsubscribe from further messages, please visit https://sourceforge.net/auth/prefs/
Related
Bugs:
#42