search

ksamuk / pixy

Software for painlessly estimating average nucleotide diversity within and between populations
https://pixy.readthedocs.io/
MIT License
117 stars 14 forks source link

pi output gives NA values when there is viable data #45

Closed anitato closed 3 years ago

anitato commented 3 years ago

I'm not 100% sure this is a bug, but I noticed that I'm getting NA values for certain sites that contain data. Not sure why this is happening?

For example, a NA occurs at site 1539: main jcf7190000000318 1530 1530 0.0 1 0 10 15 main jcf7190000000318 1531 1531 0.0 1 0 10 15 main jcf7190000000318 1532 1532 0.0 1 0 10 15 main jcf7190000000318 1533 1533 0.0 1 0 10 15 main jcf7190000000318 1534 1534 0.0 1 0 10 15 main jcf7190000000318 1535 1535 0.0 1 0 10 15 main jcf7190000000318 1536 1536 0.0 1 0 10 15 main jcf7190000000318 1537 1537 0.0 1 0 10 15 main jcf7190000000318 1538 1538 0.0 1 0 10 15 main jcf7190000000318 1539 1539 NA 0 NA NA NA main jcf7190000000318 1540 1540 0.6 1 6 10 15 main jcf7190000000318 1541 1541 0.0 1 0 10 15 main jcf7190000000318 1542 1542 0.6 1 6 10 15 main jcf7190000000318 1543 1543 0.0 1 0 10 15 main jcf7190000000318 1544 1544 0.0 1 0 45 10 main jcf7190000000318 1545 1545 0.0 1 0 45 10 main jcf7190000000318 1546 1546 0.0 1 0 45 10

My command: pixy --stats pi --vcf /path/allsites.g.vcf.gz --populations /path/pixy_population --window_size 1 --output_folder /path/pi_div --output_prefix pixy_output_test

VCF subset: #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT DDP10PL:illumina DDP1PL:illumina DDP2PL:illumina DDP3PL:illumina DDP4PL:illumina DDP5PL:illumina DDP6PL:illumina DDP7PL:illumina DDP8PL:illumina DDP9PL:illumina jcf7190000000318 1530 . A . . . AN=5;DP=97 GT:AD:DP:RGQ 0:7:7:99 .:14:14:0 .:10:10:0 0:12:12:99 .:12:12:0 .:9:9:0 0:6:6:99 0:9:9:99 .:6:6:0 0:12:12:99 jcf7190000000318 1531 . A . . . AN=5;DP=97 GT:AD:DP:RGQ 0:7:7:99 .:14:14:0 .:10:10:0 0:12:12:99 .:12:12:0 .:9:9:0 0:6:6:99 0:9:9:99 .:6:6:0 0:12:12:99 jcf7190000000318 1532 . A . . . AN=5;DP=99 GT:AD:DP:RGQ 0:7:7:99 .:14:14:0 .:10:10:0 0:14:14:99 .:12:12:0 .:9:9:0 0:6:6:99 0:9:9:99 .:6:6:0 0:12:12:99 jcf7190000000318 1533 . A . . . AN=5;DP=99 GT:AD:DP:RGQ 0:7:7:99 .:14:14:0 .:10:10:0 0:14:14:99 .:11:12:0 .:9:9:0 0:6:6:99 0:9:9:99 .:6:6:0 0:12:12:99 jcf7190000000318 1534 . T . . . AN=5;DP=101 GT:AD:DP:RGQ 0:7:7:99 .:13:14:0 .:10:10:0 0:14:14:99 .:11:12:0 .:9:9:0 0:6:6:99 0:10:11:99 .:6:6:0 0:11:12:99 jcf7190000000318 1535 . T . . . AN=5;DP=101 GT:AD:DP:RGQ 0:7:7:99 .:14:14:0 .:10:10:0 0:14:14:99 .:12:12:0 .:9:9:0 0:6:6:99 0:10:11:99 .:6:6:0 0:12:12:99 jcf7190000000318 1536 . G . . . AN=5;DP=101 GT:AD:DP:RGQ 0:7:7:99 .:14:14:0 .:10:10:0 0:14:14:99 .:12:12:0 .:9:9:0 0:6:6:99 0:11:11:99 .:6:6:0 0:12:12:99 jcf7190000000318 1537 . A . . . AN=5;DP=100 GT:AD:DP:RGQ 0:7:7:99 .:14:14:0 .:9:9:0 0:14:14:99 .:12:12:0 .:9:9:0 0:6:6:99 0:11:11:99 .:6:6:0 0:12:12:99 jcf7190000000318 1538 . A . . . AN=5;DP=100 GT:AD:DP:RGQ 0:7:7:99 .:14:14:0 .:9:9:0 0:12:14:99 .:12:12:0 .:9:9:0 0:6:6:99 0:11:11:99 .:6:6:0 0:12:12:99 jcf7190000000318 1539 . GAAT G 1528.68 . AC=4;AF=0.444;AN=9;DP=92;FS=0;MLEAC=4;MLEAF=0.444;MQ=34.84;QD=32.77;SOR=0.869 GT:AD:DP:GQ:PL 0:5,2:7:99:0,125 1:0,11:11:99:495,0 1:0,8:8:99:360,0 0:14,0:14:99:0,573 1:0,9:9:99:405,0 1:0,7:7:99:315,0 0:6,0:6:99:0,270 0:10,1:11:99:0,432 .:0,0:.:.:. 0:13,0:13:99:0,553 jcf7190000000318 1540 . A G 1043.8 . AC=2;AF=0.4;AN=5;DP=101;FS=0;MLEAC=2;MLEAF=0.4;MQ=28.95;QD=30.77;SOR=0.859 GT:AD:DP:GQ:PL 0:5,2:7:99:0,154 .:0,0:.:.:. .:0,0:.:.:. 0:14,0:14:99:0,578 .:0,0:.:.:. .:0,0:.:.:. 0:6,0:6:99:0,266 1:0,10:10:99:450,0 .:0,0:.:.:. 1:0,14:14:99:619,0 jcf7190000000318 1541 . A . . . AN=5;DP=102 GT:AD:DP:RGQ 0:5:7:99 .:0:14:0 .:0:9:0 0:14:14:99 .:0:11:0 .:0:9:0 0:6:6:99 0:11:11:99 .:0:6:0 0:15:15:99 jcf7190000000318 1542 . T G 1043.8 . AC=2;AF=0.4;AN=5;DP=100;FS=0;MLEAC=2;MLEAF=0.4;MQ=29.09;QD=28.07;SOR=0.859 GT:AD:DP:GQ:PL 0:5,2:7:99:0,159 .:0,0:.:.:. .:0,0:.:.:. 0:14,0:14:99:0,547 .:0,0:.:.:. .:0,0:.:.:. 0:6,0:6:99:0,234 1:0,10:10:99:450,0 .:0,0:.:.:. 1:0,14:14:99:619,0 jcf7190000000318 1543 . G . . . AN=5;DP=102 GT:AD:DP:RGQ 0:5:7:99 .:0:14:0 .:0:9:0 0:14:14:99 .:0:11:0 .:0:9:0 0:6:6:99 0:11:11:99 .:0:6:0 0:15:15:99 jcf7190000000318 1544 . A . . . AN=10;DP=102 GT:AD:DP:RGQ 0:7:7:99 0:14:14:99 0:9:9:99 0:14:14:99 0:11:11:99 0:9:9:99 0:6:6:99 0:11:11:99 0:5:6:99 0:15:15:99 jcf7190000000318 1545 . C . . . AN=10;DP=100 GT:AD:DP:RGQ 0:7:7:99 0:14:14:99 0:9:9:99 0:14:14:99 0:11:11:99 0:9:9:99 0:6:6:99 0:11:11:99 0:6:6:99 0:13:13:99 jcf7190000000318 1546 . T . . . AN=10;DP=102 GT:AD:DP:RGQ 0:7:7:99 0:14:14:99 0:9:9:99 0:14:14:99 0:11:11:99 0:9:9:99 0:6:6:99 0:10:11:99 0:6:6:99 0:15:15:99

Populations file: DDP10PL:illumina main DDP1PL:illumina main DDP2PL:illumina main DDP3PL:illumina main DDP4PL:illumina main DDP5PL:illumina main DDP6PL:illumina main DDP7PL:illumina main DDP8PL:illumina main DDP9PL:illumina main

Using: Linux

ksamuk commented 3 years ago

Hi there, I can certainly take a look, but can you post a link to a formatted version of those files?

anitato commented 3 years ago

Yes, sorry, the formatting disappeared once I submitted! Let me make an edit to the original post.

ksamuk commented 3 years ago

I think I see the issue, Site 1539 is an INDEL (REF = GAAT, ALT= G), which are filtered out prior to running pixy (pixy only looks at biallelic SNPs). So there would be no data at that site post-filtration. I've been looking at some options to allow for non-SNP variation, but there aren't many (any?) published algorithms that account for differences in sequence length with structural variants, apologies!

anitato commented 3 years ago

Ah, I don't know why I missed that! I must be blind. Thank you so much! Yes, indels should definitely be ignored in that case, and the larger problem of accounting for indel variation is non-trivial for sure.

ksamuk commented 3 years ago

You're welcome!