Examples of org.broadinstitute.gatk.utils.GenomeLoc

org.broadinstitute.gatk.utils.GenomeLoc
Created by IntelliJ IDEA. User: mdepristo Date: Mar 2, 2009 Time: 8:50:11 AM Genome location representation. It is *** 1 *** based closed. Note that GenomeLocs start and stop values can be any positive or negative number, by design. Bound validation is a feature of the GenomeLocParser, and not a fundamental constraint of the GenomeLoc

        Genotype sam_2_2_1_eval = GenotypeBuilder.create("test2_sample2", Arrays.asList(reference_A,alt_C));


        Genotype sam_2_1_1_truth = GenotypeBuilder.create("test2_sample1", Arrays.asList(reference_A,reference_A));
        Genotype sam_2_2_1_truth = GenotypeBuilder.create("test2_sample2", Arrays.asList(reference_A,reference_A));


        GenomeLoc loc = genomeLocParser.createGenomeLoc("chr1", 3, 3);
        VariantContextBuilder eval_1_builder = new VariantContextBuilder();
        VariantContextBuilder truth_1_builder = new VariantContextBuilder();


        eval_1_builder.alleles(Arrays.asList(reference_A,alt_C));
        truth_1_builder.alleles(Arrays.asList(reference_A,alt_C));
        eval_1_builder.genotypes(Arrays.asList(sam_2_1_1_eval,sam_2_2_1_eval));
        truth_1_builder.genotypes(Arrays.asList(sam_2_1_1_truth,sam_2_2_1_truth));


        eval_1_builder.loc(loc.getContig(),loc.getStart(),loc.getStop());
        truth_1_builder.loc(loc.getContig(),loc.getStart(),loc.getStop());


        Pair<VariantContext,VariantContext> testDataSite1 = new Pair<VariantContext, VariantContext>(eval_1_builder.make(),truth_1_builder.make());


        reference_A = Allele.create(BaseUtils.Base.A.base,true);
        Allele alt_T = Allele.create(BaseUtils.Base.T.base);


        // site 2 -
        //  sample 1: no-call/hom-ref
        //  sample 2: hom-var/hom-var


        Genotype sam_2_1_2_eval = GenotypeBuilder.create("test2_sample1",Arrays.asList(Allele.NO_CALL,Allele.NO_CALL));
        Genotype sam_2_2_2_eval = GenotypeBuilder.create("test2_sample2",Arrays.asList(alt_T,alt_T));
        Genotype sam_2_1_2_truth = GenotypeBuilder.create("test2_sample1",Arrays.asList(reference_A,reference_A));
        Genotype sam_2_2_2_truth = GenotypeBuilder.create("test2_sample2",Arrays.asList(alt_T,alt_T));


        loc = genomeLocParser.createGenomeLoc("chr1", 4, 4);
        eval_1_builder = new VariantContextBuilder();
        truth_1_builder = new VariantContextBuilder();


        eval_1_builder.loc(loc.getContig(),loc.getStart(),loc.getStop());
        truth_1_builder.loc(loc.getContig(),loc.getStart(),loc.getStop());
        eval_1_builder.alleles(Arrays.asList(reference_A,alt_T));
        truth_1_builder.alleles(Arrays.asList(reference_A,alt_T));
        eval_1_builder.genotypes(Arrays.asList(sam_2_1_2_eval,sam_2_2_2_eval));
        truth_1_builder.genotypes(Arrays.asList(sam_2_1_2_truth,sam_2_2_2_truth));


        Pair<VariantContext,VariantContext> testDataSite2 = new Pair<VariantContext, VariantContext>(eval_1_builder.make(),truth_1_builder.make());


        Allele alt_G = Allele.create(BaseUtils.Base.G.base);


        // site 3 -
        //  sample 1: alleles do not match
        //  sample 2: het/het
        Genotype sam_2_1_3_eval = GenotypeBuilder.create("test2_sample1",Arrays.asList(alt_G,alt_T));
        Genotype sam_2_2_3_eval = GenotypeBuilder.create("test2_sample2",Arrays.asList(reference_A,alt_T));
        Genotype sam_2_1_3_truth = GenotypeBuilder.create("test2_sample1",Arrays.asList(alt_T,alt_T));
        Genotype sam_2_2_3_truth = GenotypeBuilder.create("test2_sample2",Arrays.asList(reference_A,alt_T));


        loc = genomeLocParser.createGenomeLoc("chr1",5,5);
        eval_1_builder = new VariantContextBuilder();
        truth_1_builder = new VariantContextBuilder();
        eval_1_builder.loc(loc.getContig(),loc.getStart(),loc.getStop());
        truth_1_builder.loc(loc.getContig(),loc.getStart(),loc.getStop());
        eval_1_builder.alleles(Arrays.asList(reference_A,alt_T,alt_G));
        truth_1_builder.alleles(Arrays.asList(reference_A,alt_T));
        eval_1_builder.genotypes(Arrays.asList(sam_2_1_3_eval,sam_2_2_3_eval));
        truth_1_builder.genotypes(Arrays.asList(sam_2_1_3_truth,sam_2_2_3_truth));


        Pair<VariantContext,VariantContext> testDataSite3 = new Pair<VariantContext, VariantContext>(eval_1_builder.make(),truth_1_builder.make());


        // site 4 -
        //  sample 1: unavailable/het
        //  sample 2: unavailable/ref
        Genotype sam_2_1_4_eval = GenotypeBuilder.create("test2_sample1",new ArrayList<Allele>(0));
        Genotype sam_2_2_4_eval = GenotypeBuilder.create("test2_sample2",new ArrayList<Allele>(0));
        Genotype sam_2_1_4_truth = GenotypeBuilder.create("test2_sample1",Arrays.asList(reference_A,alt_T));
        Genotype sam_2_2_4_truth = GenotypeBuilder.create("test2_sample2",Arrays.asList(reference_A,reference_A));


        loc = genomeLocParser.createGenomeLoc("chr1",6,6);
        eval_1_builder = new VariantContextBuilder();
        truth_1_builder = new VariantContextBuilder();
        eval_1_builder.loc(loc.getContig(),loc.getStart(),loc.getStop());
        truth_1_builder.loc(loc.getContig(),loc.getStart(),loc.getStop());
        eval_1_builder.alleles(Arrays.asList(reference_A,alt_T));
        truth_1_builder.alleles(Arrays.asList(reference_A,alt_T));
        eval_1_builder.genotypes(Arrays.asList(sam_2_1_4_eval,sam_2_2_4_eval));
        truth_1_builder.genotypes(Arrays.asList(sam_2_1_4_truth,sam_2_2_4_truth));


        Pair<VariantContext,VariantContext> testDataSite4 = new Pair<VariantContext, VariantContext>(eval_1_builder.make(),truth_1_builder.make());


        // site 5 -
        //  sample 1: hom-var/no-call
        //  sample 2: het/het
        Genotype sam_2_1_5_eval = GenotypeBuilder.create("test2_sample1",Arrays.asList(alt_C,alt_C));
        Genotype sam_2_2_5_eval = GenotypeBuilder.create("test2_sample2",Arrays.asList(reference_A,alt_C));
        Genotype sam_2_1_5_truth = GenotypeBuilder.create("test2_sample1",Arrays.asList(Allele.NO_CALL,Allele.NO_CALL));
        Genotype sam_2_2_5_truth = GenotypeBuilder.create("test2_sample2",Arrays.asList(reference_A,alt_C));


        loc = genomeLocParser.createGenomeLoc("chr1",7,7);
        eval_1_builder = new VariantContextBuilder();
        truth_1_builder = new VariantContextBuilder();
        eval_1_builder.loc(loc.getContig(),loc.getStart(),loc.getStop());
        truth_1_builder.loc(loc.getContig(),loc.getStart(),loc.getStop());
        eval_1_builder.alleles(Arrays.asList(reference_A,alt_C));
        truth_1_builder.alleles(Arrays.asList(reference_A,alt_C));
        eval_1_builder.genotypes(Arrays.asList(sam_2_1_5_eval,sam_2_2_5_eval));
        truth_1_builder.genotypes(Arrays.asList(sam_2_1_5_truth,sam_2_2_5_truth));

View Full Code Here

    public void testTrimTo(final Map<Haplotype,AssemblyResult> haplotypesAndResultSets, final ActiveRegion original) {
        final AssemblyResultSet subject = new AssemblyResultSet();
        for (final Map.Entry<Haplotype,AssemblyResult> entry : haplotypesAndResultSets.entrySet())
            subject.add(entry.getKey(),entry.getValue());
        subject.setRegionForGenotyping(original);
        final GenomeLoc originalLocation = original.getExtendedLoc();
        final int length = originalLocation.size();
        final GenomeLoc newLocation = originalLocation.setStop(originalLocation.setStart(originalLocation,originalLocation.getStart() + length / 2),originalLocation.getStop() - length / 2);
        final ActiveRegion newRegion = original.trim(newLocation);


        final Map<Haplotype,Haplotype> originalHaplotypesByTrimmed = new HashMap<>(haplotypesAndResultSets.size());
        for (final Haplotype h : haplotypesAndResultSets.keySet())
            originalHaplotypesByTrimmed.put(h.trim(newRegion.getExtendedLoc()), h);


        final AssemblyResultSet trimmed = subject.trimTo(newRegion);


        Assert.assertFalse(subject.wasTrimmed());
        Assert.assertTrue(trimmed.wasTrimmed());


        for (final Haplotype h : trimmed.getHaplotypeList()) {
            Assert.assertEquals(h.getGenomeLocation(),newLocation);
            Assert.assertEquals(h.getBases().length,newLocation.size());
        }
    }

View Full Code Here

        if ( refHaplotype.length() != activeRegion.getExtendedLoc().size() ) throw new IllegalArgumentException("refHaplotype " + refHaplotype.length() + " and activeRegion location size " + activeRegion.getLocation().size() + " are different");
        if ( ploidyModel == null) throw new IllegalArgumentException("the ploidy model cannot be null");
        if ( model == null) throw new IllegalArgumentException("the genotyping model cannot be null");
        final int ploidy = ploidyModel.samplePloidy(0); // the first sample = the only sample in reference-confidence mode.


        final GenomeLoc refSpan = activeRegion.getLocation();
        final List<ReadBackedPileup> refPileups = getPileupsOverReference(refHaplotype, calledHaplotypes, paddedReferenceLoc, activeRegion, refSpan, readLikelihoods);
        final byte[] ref = refHaplotype.getBases();
        final List<VariantContext> results = new ArrayList<>(refSpan.size());
        final String sampleName = readLikelihoods.sampleAt(0);


        final int globalRefOffset = refSpan.getStart() - activeRegion.getExtendedLoc().getStart();
        for ( final ReadBackedPileup pileup : refPileups ) {
            final GenomeLoc curPos = pileup.getLocation();
            final int offset = curPos.getStart() - refSpan.getStart();


            final VariantContext overlappingSite = getOverlappingVariantContext(curPos, variantCalls);
            if ( overlappingSite != null && overlappingSite.getStart() == curPos.getStart() ) {
                    results.add(overlappingSite);
            } else {
                // otherwise emit a reference confidence variant context
                final int refOffset = offset + globalRefOffset;
                final byte refBase = ref[refOffset];
                final RefVsAnyResult homRefCalc = calcGenotypeLikelihoodsOfRefVsAny(sampleName,ploidy,model,pileup, refBase, (byte)6, null);
                homRefCalc.capByHomRefLikelihood();


                final Allele refAllele = Allele.create(refBase, true);
                final List<Allele> refSiteAlleles = Arrays.asList(refAllele, GATKVariantContextUtils.NON_REF_SYMBOLIC_ALLELE);
                final VariantContextBuilder vcb = new VariantContextBuilder("HC", curPos.getContig(), curPos.getStart(), curPos.getStart(), refSiteAlleles);
                final GenotypeBuilder gb = new GenotypeBuilder(sampleName, GATKVariantContextUtils.homozygousAlleleList(refAllele, ploidy));
                gb.AD(homRefCalc.AD_Ref_Any);
                gb.DP(homRefCalc.getDP());


                // genotype likelihood calculation

View Full Code Here

        if ( read.getReferenceIndex() == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX ) {
            cleanAndCallMap(ref, read, metaDataTracker, null);
            return 0;
        }


        GenomeLoc readLoc = getToolkit().getGenomeLocParser().createGenomeLoc(read);
        // hack to get around unmapped reads having screwy locations
        if ( readLoc.getStop() == 0 )
            readLoc = getToolkit().getGenomeLocParser().createGenomeLoc(readLoc.getContig(), readLoc.getStart(), readLoc.getStart());


        if ( readLoc.isBefore(currentInterval) ) {
            if ( !sawReadInCurrentInterval )
                emit(read);
            else
                readsNotToClean.add(read);
        }
        else if ( readLoc.overlapsP(currentInterval) ) {
            sawReadInCurrentInterval = true;


            if ( doNotTryToClean(read) ) {
                readsNotToClean.add(read);
            } else {

View Full Code Here

        // find the alternate allele(s) that we should be using
        final List<Allele> alleles = getFinalAllelesToUse(tracker, ref, allAllelesToUse, GLs);
        if (alleles == null || alleles.isEmpty() || (alleles.size() == 1 && alleles.get(0).isReference()))
            return null;
        // start making the VariantContext
        final GenomeLoc loc = ref.getLocus();
        final int endLoc = getEndLocation(tracker, ref, alleles);


        final VariantContextBuilder builder = new VariantContextBuilder("UG_call", loc.getContig(), loc.getStart(), endLoc, alleles);
        builder.alleles(alleles);


        final HashMap<String, Object> attributes = new HashMap<String, Object>();


        if (UAC.referenceSampleName != null && perLaneErrorModels != null)

View Full Code Here

        // TODO -- it would be nice if we could use indels from 454/Ion reads as alternate consenses
    }


    private void cleanAndCallMap(ReferenceContext ref, GATKSAMRecord read, RefMetaDataTracker metaDataTracker, GenomeLoc readLoc) {
        if ( readsToClean.size() > 0 ) {
            GenomeLoc earliestPossibleMove = getToolkit().getGenomeLocParser().createGenomeLoc(readsToClean.getReads().get(0));
            if ( manager.canMoveReads(earliestPossibleMove) )
                clean(readsToClean);
        }
        knownIndelsToTry.clear();
        indelRodsSeen.clear();

View Full Code Here

        return sum + value;
    }


    public void onTraversalDone(Integer result) {
        if ( readsToClean.size() > 0 ) {
            GenomeLoc earliestPossibleMove = getToolkit().getGenomeLocParser().createGenomeLoc(readsToClean.getReads().get(0));
            if ( manager.canMoveReads(earliestPossibleMove) )
                clean(readsToClean);
            emitReadLists();
        } else if ( readsNotToClean.size() > 0 ) {
            emitReadLists();

View Full Code Here

        final Metrics noData = new Metrics(0.5, 0.0, 0.0, 0, 20);
        final Metrics unknown = new Metrics(0.5, 30.0, 120000.0, 2500, 20);


        final Metrics[] array = {unmappable, highGC, lowGC, unsequenceable, noData, unknown};


        final GenomeLoc testInterval = new UnvalidatingGenomeLoc("chr1", 0, 10000, 20000);
        final GenomeLoc smallInterval = new UnvalidatingGenomeLoc("chr1", 0, 1, 4);




        Assert.assertNotEquals(tool.checkMappability(unmappable), "");
        Assert.assertNotEquals(tool.checkGCContent(highGC), "");
        Assert.assertNotEquals(tool.checkGCContent(lowGC), "");

View Full Code Here

    // TODO -- WARNING DOESN'T WORK WITH REDUCED READS
    //
    private Map<String, Integer> countConsensusAlleles(ReferenceContext ref,
                                                       Map<String, AlignmentContext> contexts,
                                                       AlignmentContextUtils.ReadOrientation contextType) {
        final GenomeLoc loc = ref.getLocus();
        HashMap<String, Integer> consensusIndelStrings = new HashMap<String, Integer>();


        int insCount = 0, delCount = 0;
        // quick check of total number of indels in pileup
        for ( Map.Entry<String, AlignmentContext> sample : contexts.entrySet() ) {
            final AlignmentContext context = AlignmentContextUtils.stratify(sample.getValue(), contextType);


            final ReadBackedPileup indelPileup = context.getBasePileup();
            insCount += indelPileup.getNumberOfInsertionsAfterThisElement();
            delCount += indelPileup.getNumberOfDeletionsAfterThisElement();
        }


        if ( insCount < minIndelCountForGenotyping && delCount < minIndelCountForGenotyping )
            return Collections.emptyMap();


        for (Map.Entry<String, AlignmentContext> sample : contexts.entrySet()) {
            // todo -- warning, can be duplicating expensive partition here
            AlignmentContext context = AlignmentContextUtils.stratify(sample.getValue(), contextType);


            final ReadBackedPileup indelPileup = context.getBasePileup();


            final int nIndelReads = indelPileup.getNumberOfInsertionsAfterThisElement() + indelPileup.getNumberOfDeletionsAfterThisElement();
            final int nReadsOverall = indelPileup.getNumberOfElements();


            if ( nIndelReads == 0 || (nIndelReads / (1.0 * nReadsOverall)) < minFractionInOneSample ) {
                continue;
            }


            for (PileupElement p : indelPileup) {
                final GATKSAMRecord read = ReadClipper.hardClipAdaptorSequence(p.getRead());
                if (read == null)
                    continue;
                if (ReadUtils.is454Read(read)) {
                    continue;
                }


                if ( p.isBeforeInsertion() ) {
                    final String insertionBases = p.getBasesOfImmediatelyFollowingInsertion();
                    // edge case: ignore a deletion immediately preceding an insertion as p.getBasesOfImmediatelyFollowingInsertion() returns null [EB]
                    if ( insertionBases == null )
                        continue;


                    boolean foundKey = false;
                    // copy of hashmap into temp arrayList
                    ArrayList<Pair<String,Integer>> cList = new ArrayList<Pair<String,Integer>>();
                    for (Map.Entry<String, Integer> s : consensusIndelStrings.entrySet()) {
                        cList.add(new Pair<String, Integer>(s.getKey(), s.getValue()));
                    }


                    if (read.getAlignmentEnd() == loc.getStart()) {
                        // first corner condition: a read has an insertion at the end, and we're right at the insertion.
                        // In this case, the read could have any of the inserted bases and we need to build a consensus


                        for (int k=0; k < cList.size(); k++) {
                            String s = cList.get(k).getFirst();
                            int cnt = cList.get(k).getSecond();
                            // case 1: current insertion is prefix of indel in hash map
                            if (s.startsWith(insertionBases)) {
                                cList.set(k,new Pair<String, Integer>(s,cnt+1));
                                foundKey = true;
                            }
                            else if (insertionBases.startsWith(s)) {
                                // case 2: indel stored in hash table is prefix of current insertion
                                // In this case, new bases are new key.
                                foundKey = true;
                                cList.set(k,new Pair<String, Integer>(insertionBases,cnt+1));
                            }
                        }
                        if (!foundKey)
                            // none of the above: event bases not supported by previous table, so add new key
                            cList.add(new Pair<String, Integer>(insertionBases,1));


                    }
                    else if (read.getAlignmentStart() == loc.getStart()+1) {
                        // opposite corner condition: read will start at current locus with an insertion
                        for (int k=0; k < cList.size(); k++) {
                            String s = cList.get(k).getFirst();
                            int cnt = cList.get(k).getSecond();
                            if (s.endsWith(insertionBases)) {

View Full Code Here

        final int stepSize = 200;
        final int nReadsToUse = 5;


        for ( int startI = start; startI < end; startI += stepSize) {
            final int endI = startI + windowSize;
            final GenomeLoc refLoc = genomeLocParser.createGenomeLoc(contig, startI, endI);
            tests.add(new Object[]{new ReadThreadingAssembler(), refLoc, nReadsToUse});
        }


        return tests.toArray(new Object[][]{});
    }

View Full Code Here

0 1 2 3 4 5 6 7 8 9

TOP

Related Classes of org.broadinstitute.gatk.utils.GenomeLoc

org.broad.igv.plugin.mongovariant.VariantReviewSource

org.broadinstitute.gatk.engine.contexts.AlignmentContextUtils

org.broadinstitute.gatk.engine.datasources.providers.AllLocusViewUnitTest

org.broadinstitute.gatk.engine.datasources.providers.CoveredLocusViewUnitTest

org.broadinstitute.gatk.engine.datasources.providers.IntervalOverlappingRODsFromStream

org.broadinstitute.gatk.engine.datasources.providers.IntervalReferenceOrderedViewUnitTest

org.broadinstitute.gatk.engine.datasources.providers.IntervalReferenceOrderedViewUnitTest$TribbleIteratorFromCollection

org.broadinstitute.gatk.engine.datasources.providers.LocusReferenceView

org.broadinstitute.gatk.engine.datasources.providers.LocusReferenceViewUnitTest

org.broadinstitute.gatk.engine.datasources.providers.LocusViewTemplate

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.