Package org.broadinstitute.gatk.tools.walkers.indels

Examples of org.broadinstitute.gatk.tools.walkers.indels.PairHMMIndelErrorModel


        this.maxQualityScore = UAC.maxQualityScore;
        this.minQualityScore = UAC.minQualityScore;
        this.phredScaledPrior = UAC.phredScaledPrior;
        log10minPower = Math.log10(UAC.minPower);

        PairHMMIndelErrorModel pairModel = null;
        LinkedHashMap<Allele, Haplotype> haplotypeMap = null;
        double[][] perReadLikelihoods = null;

        double[] model = new double[maxQualityScore+1];
        Arrays.fill(model,Double.NEGATIVE_INFINITY);

        boolean hasCalledAlleles = false;

        final PerReadAlleleLikelihoodMap perReadAlleleLikelihoodMap = new PerReadAlleleLikelihoodMap();
        if (refSampleVC != null) {

            for (Allele allele : refSampleVC.getAlleles()) {
                if (allele.isCalled()) {
                    hasCalledAlleles = true;
                    break;
                }
            }
            haplotypeMap = new LinkedHashMap<Allele, Haplotype>();
            if (refSampleVC.isIndel()) {
                pairModel = new PairHMMIndelErrorModel(UAC.INDEL_GAP_OPEN_PENALTY, UAC.INDEL_GAP_CONTINUATION_PENALTY,
                        UAC.OUTPUT_DEBUG_INDEL_INFO, UAC.pairHMM);
                IndelGenotypeLikelihoodsCalculationModel.getHaplotypeMapFromAlleles(refSampleVC.getAlleles(), refContext, refContext.getLocus(), haplotypeMap); // will update haplotypeMap adding elements
            }
        }

        double p = QualityUtils.qualToErrorProbLog10((byte)(maxQualityScore-minQualityScore));
        if (refSamplePileup == null || refSampleVC == null  || !hasCalledAlleles) {
            for (byte q=minQualityScore; q<=maxQualityScore; q++) {
                // maximum uncertainty if there's no ref data at site
                model[q] = p;
            }
            this.refDepth = 0;
        }
        else {
            hasData = true;
            int matches = 0;
            int coverage = 0;

            Allele refAllele = refSampleVC.getReference();

            if ( refSampleVC.isIndel()) {
                //perReadLikelihoods = new double[readCounts.length][refSampleVC.getAlleles().size()];
                final int eventLength = IndelGenotypeLikelihoodsCalculationModel.getEventLength(refSampleVC.getAlleles());
                if (!haplotypeMap.isEmpty())
                    perReadLikelihoods = pairModel.computeGeneralReadHaplotypeLikelihoods(refSamplePileup,haplotypeMap,refContext, eventLength, perReadAlleleLikelihoodMap);
            }
            int idx = 0;
            for (PileupElement refPileupElement : refSamplePileup) {
                if (DEBUG)
                    System.out.println(refPileupElement.toString());
View Full Code Here


    protected GeneralPloidyIndelGenotypeLikelihoodsCalculationModel(final UnifiedArgumentCollection UAC, final Logger logger) {
        super(UAC, logger);


        pairModel = new PairHMMIndelErrorModel(UAC.INDEL_GAP_OPEN_PENALTY, UAC.INDEL_GAP_CONTINUATION_PENALTY,
                UAC.OUTPUT_DEBUG_INDEL_INFO, UAC.pairHMM);
        haplotypeMap = new LinkedHashMap<Allele, Haplotype>();
    }
View Full Code Here


    protected IndelGenotypeLikelihoodsCalculationModel(final UnifiedArgumentCollection UAC,
                                                       final Logger logger) {
        super(UAC, logger);
        pairModel = new PairHMMIndelErrorModel(UAC.INDEL_GAP_OPEN_PENALTY, UAC.INDEL_GAP_CONTINUATION_PENALTY,
                UAC.OUTPUT_DEBUG_INDEL_INFO, UAC.pairHMM);
        DEBUG = UAC.OUTPUT_DEBUG_INDEL_INFO;
        haplotypeMap = new LinkedHashMap<Allele, Haplotype>();
        ignoreSNPAllelesWhenGenotypingIndels = UAC.IGNORE_SNP_ALLELES;
    }
View Full Code Here

TOP

Related Classes of org.broadinstitute.gatk.tools.walkers.indels.PairHMMIndelErrorModel

Copyright © 2018 www.massapicom. All rights reserved.
All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.