@ Override public double value ( double alpha ) { double logLikelihood = 0.0 ; for ( ReadCounts rc : readCounts . values ( ) ) { double refProb = rc . getRefFreq ( ) ; double pPosition = 0.0 ; for ( Genotype trueGenotype1 : Genotype . values ( ) ) { double pGenotype1 = pGenotype ( trueGenotype1 , refProb ) ; for ( Genotype trueGenotype2 : Genotype . values ( ) ) { double pGenotype2 = pGenotype ( trueGenotype2 , refProb ) ; double pObsGivenGenotype = 1.0 ; for ( ReadQualityCount rqc : rc . getReadQualityCounts ( ) ) { Base base = rqc . getBase ( ) ; double pErr = phredToProb ( rqc . getQuality ( ) ) ; double pObs = ( ( 1.0 - alpha ) * probObsGivenTruth ( base , trueGenotype1 , ReadStatus . CORRECT ) + ( alpha ) * probObsGivenTruth ( base , trueGenotype2 , ReadStatus . CORRECT ) ) * ( 1.0 - pErr ) + ( ( 1.0 - alpha ) * probObsGivenTruth ( base , trueGenotype1 , ReadStatus . ERROR ) + ( alpha ) * probObsGivenTruth ( base , trueGenotype2 , ReadStatus . ERROR ) ) * pErr ; pObsGivenGenotype *= Math . pow ( pObs , rqc . getCount ( ) ) ; } pPosition += pObsGivenGenotype * pGenotype1 * pGenotype2 ; } } logLikelihood += Math . log ( pPosition ) ; } return logLikelihood ; }
Compute the likelihood of a contaminant fraction alpha .