Minor bug fix in empiricalMu prior calculation in VQSR.

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@5412 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
rpoplin 2011-03-10 00:42:38 +00:00
parent 0b45de14ed
commit ff7edc4493
3 changed files with 11 additions and 6 deletions

View File

@ -155,9 +155,9 @@ public class RecalDataManager {
} }
/** /**
* Loop over all the collapsed tables and turn the recalDatums found there into an empricial quality score * Loop over all the collapsed tables and turn the recalDatums found there into an empirical quality score
* that will be used in the sequential calculation in TableRecalibrationWalker * that will be used in the sequential calculation in TableRecalibrationWalker
* @param smoothing The smoothing paramter that goes into empirical quality score calculation * @param smoothing The smoothing parameter that goes into empirical quality score calculation
* @param maxQual At which value to cap the quality scores * @param maxQual At which value to cap the quality scores
*/ */
public final void generateEmpiricalQualities( final int smoothing, final int maxQual ) { public final void generateEmpiricalQualities( final int smoothing, final int maxQual ) {

View File

@ -231,16 +231,21 @@ public final class VariantGaussianMixtureModel extends VariantOptimizationModel
} }
} }
double sumWeight = 0.0;
for(int iii = 0; iii < numVariants; iii++) { for(int iii = 0; iii < numVariants; iii++) {
sumWeight += data[iii].weight;
for( int jjj = 0; jjj < numAnnotations; jjj++ ) { for( int jjj = 0; jjj < numAnnotations; jjj++ ) {
empiricalMu[jjj] += data[iii].annotations[jjj] / ((double) numVariants); empiricalMu[jjj] += data[iii].annotations[jjj] * data[iii].weight;
} }
} }
for( int jjj = 0; jjj < numAnnotations; jjj++ ) {
empiricalMu[jjj] /= sumWeight;
}
for(int iii = 0; iii < numVariants; iii++) { for(int iii = 0; iii < numVariants; iii++) {
for( int jjj = 0; jjj < numAnnotations; jjj++ ) { for( int jjj = 0; jjj < numAnnotations; jjj++ ) {
for( int ppp = 0; ppp < numAnnotations; ppp++ ) { for( int ppp = 0; ppp < numAnnotations; ppp++ ) {
if( jjj == ppp ) { if( jjj == ppp ) {
sigmaVals[jjj][ppp] = 1.0; //0.01 * numAnnotations; sigmaVals[jjj][ppp] = 1.0; //0.01 * numAnnotations;
} else { } else {
sigmaVals[jjj][ppp] = 0.0; sigmaVals[jjj][ppp] = 0.0;

View File

@ -29,13 +29,13 @@ public class VariantRecalibrationWalkersIntegrationTest extends WalkerTest {
} }
VRTest yriTrio = new VRTest("yri.trio.gatk_glftrio.intersection.annotated.filtered.chr1.vcf", VRTest yriTrio = new VRTest("yri.trio.gatk_glftrio.intersection.annotated.filtered.chr1.vcf",
"e14079c3a02c112665a6c194fd4f5d5c", // cluster file "f65c27ee40053adc72dd0bfbb628e4d7", // cluster file
"dce581b880ffb6ea39cbada1ecc95915", // tranches "dce581b880ffb6ea39cbada1ecc95915", // tranches
"c3e8a2f43656eab7d847dbf850f844a6", // recalVCF "c3e8a2f43656eab7d847dbf850f844a6", // recalVCF
"50f752a72643db9ad0aa94b3fc4e23d6"); // cut VCF "50f752a72643db9ad0aa94b3fc4e23d6"); // cut VCF
VRTest lowPass = new VRTest("lowpass.N3.chr1.raw.vcf", VRTest lowPass = new VRTest("lowpass.N3.chr1.raw.vcf",
"b0c0f8c8d9fe3d1ed2bde0b1eb82a22d", // cluster file "bda8f17cfc19d23e7e51f99e547f4b3d", // cluster file
"66edae83c50f4e8601fef7fafba774af", // tranches "66edae83c50f4e8601fef7fafba774af", // tranches
"0123537e373657386068a534c0f5c91b", // recalVCF "0123537e373657386068a534c0f5c91b", // recalVCF
"2172368e8585841e5ad96c95d0827c4b"); // cut VCF "2172368e8585841e5ad96c95d0827c4b"); // cut VCF