f766ff086ef77a4730d568e9343a990c0033cf4b tdreszer Mon May 20 12:58:37 2013 -0700 Fixed accounting bug where chrX reference haplotypes could have >100 percent frequency. diff --git src/lib/vcfBits.c src/lib/vcfBits.c index d0c5da6..c1dc612 100644 --- src/lib/vcfBits.c +++ src/lib/vcfBits.c @@ -450,44 +450,49 @@ if (bitReadOne(vBits->bits,vBitsSlot(vBits,genoIx,haploIx,1))) { bitSetOne( bits,hBitsSlot(hBits,variantIx, 1)); hBits->bitsOn++; } } } // gonna keep it? so flesh it out (ignoreRef means > 0 && < all) if (!ignoreReference || hBits->bitsOn) { hBits->genomeIx = genoIx; hBits->haploidIx = haploIx; struct vcfRecord *record = vBitsList->record; // any will do struct vcfGenotype *gt = &(record->genotypes[genoIx]); + + if (hBits->bitsOn // if including reference, then chrX could + || haploIx == 0 || !gt->isHaploid) // have unused diploid positions! + { if (gt->isHaploid || vBitsList->haplotypeSlots == 1) { // chrX will have haplotypeSlots==2 but be haploid for this subject. // Meanwhile if vBits were for homozygous only, haplotypeSlots==1 //assert(haploIx == 0); hBits->ids = lmCloneString(lm,gt->id); } else { int sz = strlen(gt->id) + 3; hBits->ids = lmAlloc(lm,sz); safef(hBits->ids,sz,"%s-%c",gt->id,'a' + haploIx); } slAddHead(&hBitsList,hBits); } + } //else // haploBitsFree(&hBits); // lm so just abandon } } slReverse(&hBitsList); return hBitsList; } /* int vcfHaploBitsOnCmp(const void *va, const void *vb) // Compare to sort haploBits based upon how many bits are on { const struct haploBits *a = *((struct haploBits **)va); const struct haploBits *b = *((struct haploBits **)vb);