Search in sources :

Example 26 with VariantContextBuilder

use of htsjdk.variant.variantcontext.VariantContextBuilder in project gatk by broadinstitute.

the class CoverageUnitTest method makeVC.

private VariantContext makeVC() {
    final GenotypesContext testGC = GenotypesContext.create(2);
    final Allele refAllele = Allele.create("A", true);
    final Allele altAllele = Allele.create("T");
    return (new VariantContextBuilder()).alleles(Arrays.asList(refAllele, altAllele)).chr("1").start(15L).stop(15L).genotypes(testGC).make();
}
Also used : Allele(htsjdk.variant.variantcontext.Allele) GenotypesContext(htsjdk.variant.variantcontext.GenotypesContext) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder)

Example 27 with VariantContextBuilder

use of htsjdk.variant.variantcontext.VariantContextBuilder in project gatk by broadinstitute.

the class EventMap method makeBlock.

/**
     * Create a block substitution out of two variant contexts that start at the same position
     *
     * vc1 can be SNP, and vc2 can then be either a insertion or deletion.
     * If vc1 is an indel, then vc2 must be the opposite type (vc1 deletion => vc2 must be an insertion)
     *
     * @param vc1 the first variant context we want to merge
     * @param vc2 the second
     * @return a block substitution that represents the composite substitution implied by vc1 and vc2
     */
protected VariantContext makeBlock(final VariantContext vc1, final VariantContext vc2) {
    Utils.validateArg(vc1.getStart() == vc2.getStart(), () -> "vc1 and 2 must have the same start but got " + vc1 + " and " + vc2);
    Utils.validateArg(vc1.isBiallelic(), "vc1 must be biallelic");
    if (!vc1.isSNP()) {
        Utils.validateArg((vc1.isSimpleDeletion() && vc2.isSimpleInsertion()) || (vc1.isSimpleInsertion() && vc2.isSimpleDeletion()), () -> "Can only merge single insertion with deletion (or vice versa) but got " + vc1 + " merging with " + vc2);
    } else {
        Utils.validateArg(!vc2.isSNP(), () -> "vc1 is " + vc1 + " but vc2 is a SNP, which implies there's been some terrible bug in the cigar " + vc2);
    }
    final Allele ref, alt;
    final VariantContextBuilder b = new VariantContextBuilder(vc1);
    if (vc1.isSNP()) {
        // we have to repair the first base, so SNP case is special cased
        if (vc1.getReference().equals(vc2.getReference())) {
            // we've got an insertion, so we just update the alt to have the prev alt
            ref = vc1.getReference();
            alt = Allele.create(vc1.getAlternateAllele(0).getDisplayString() + vc2.getAlternateAllele(0).getDisplayString().substring(1), false);
        } else {
            // we're dealing with a deletion, so we patch the ref
            ref = vc2.getReference();
            alt = vc1.getAlternateAllele(0);
            b.stop(vc2.getEnd());
        }
    } else {
        final VariantContext insertion = vc1.isSimpleInsertion() ? vc1 : vc2;
        final VariantContext deletion = vc1.isSimpleInsertion() ? vc2 : vc1;
        ref = deletion.getReference();
        alt = insertion.getAlternateAllele(0);
        b.stop(deletion.getEnd());
    }
    return b.alleles(Arrays.asList(ref, alt)).make();
}
Also used : Allele(htsjdk.variant.variantcontext.Allele) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) VariantContext(htsjdk.variant.variantcontext.VariantContext)

Example 28 with VariantContextBuilder

use of htsjdk.variant.variantcontext.VariantContextBuilder in project gatk by broadinstitute.

the class EventMap method processCigarForInitialEvents.

protected void processCigarForInitialEvents() {
    final Cigar cigar = haplotype.getCigar();
    final byte[] alignment = haplotype.getBases();
    int refPos = haplotype.getAlignmentStartHapwrtRef();
    if (refPos < 0) {
        return;
    }
    // Protection against SW failures
    final List<VariantContext> proposedEvents = new ArrayList<>();
    int alignmentPos = 0;
    for (int cigarIndex = 0; cigarIndex < cigar.numCigarElements(); cigarIndex++) {
        final CigarElement ce = cigar.getCigarElement(cigarIndex);
        final int elementLength = ce.getLength();
        switch(ce.getOperator()) {
            case I:
                {
                    if (refPos > 0) {
                        // protect against trying to create insertions/deletions at the beginning of a contig
                        final List<Allele> insertionAlleles = new ArrayList<>();
                        final int insertionStart = refLoc.getStart() + refPos - 1;
                        final byte refByte = ref[refPos - 1];
                        if (BaseUtils.isRegularBase(refByte)) {
                            insertionAlleles.add(Allele.create(refByte, true));
                        }
                        if (cigarIndex == 0 || cigarIndex == cigar.numCigarElements() - 1) {
                        // if the insertion isn't completely resolved in the haplotype, skip it
                        // note this used to emit SYMBOLIC_UNASSEMBLED_EVENT_ALLELE but that seems dangerous
                        } else {
                            byte[] insertionBases = {};
                            // add the padding base
                            insertionBases = ArrayUtils.add(insertionBases, ref[refPos - 1]);
                            insertionBases = ArrayUtils.addAll(insertionBases, Arrays.copyOfRange(alignment, alignmentPos, alignmentPos + elementLength));
                            if (BaseUtils.isAllRegularBases(insertionBases)) {
                                insertionAlleles.add(Allele.create(insertionBases, false));
                            }
                        }
                        if (insertionAlleles.size() == 2) {
                            // found a proper ref and alt allele
                            proposedEvents.add(new VariantContextBuilder(sourceNameToAdd, refLoc.getContig(), insertionStart, insertionStart, insertionAlleles).make());
                        }
                    }
                    alignmentPos += elementLength;
                    break;
                }
            case S:
                {
                    alignmentPos += elementLength;
                    break;
                }
            case D:
                {
                    if (refPos > 0) {
                        // protect against trying to create insertions/deletions at the beginning of a contig
                        // add padding base
                        final byte[] deletionBases = Arrays.copyOfRange(ref, refPos - 1, refPos + elementLength);
                        final List<Allele> deletionAlleles = new ArrayList<>();
                        final int deletionStart = refLoc.getStart() + refPos - 1;
                        final byte refByte = ref[refPos - 1];
                        if (BaseUtils.isRegularBase(refByte) && BaseUtils.isAllRegularBases(deletionBases)) {
                            deletionAlleles.add(Allele.create(deletionBases, true));
                            deletionAlleles.add(Allele.create(refByte, false));
                            proposedEvents.add(new VariantContextBuilder(sourceNameToAdd, refLoc.getContig(), deletionStart, deletionStart + elementLength, deletionAlleles).make());
                        }
                    }
                    refPos += elementLength;
                    break;
                }
            case M:
            case EQ:
            case X:
                {
                    for (int iii = 0; iii < elementLength; iii++) {
                        final byte refByte = ref[refPos];
                        final byte altByte = alignment[alignmentPos];
                        if (refByte != altByte) {
                            // SNP!
                            if (BaseUtils.isRegularBase(refByte) && BaseUtils.isRegularBase(altByte)) {
                                final List<Allele> snpAlleles = new ArrayList<>();
                                snpAlleles.add(Allele.create(refByte, true));
                                snpAlleles.add(Allele.create(altByte, false));
                                proposedEvents.add(new VariantContextBuilder(sourceNameToAdd, refLoc.getContig(), refLoc.getStart() + refPos, refLoc.getStart() + refPos, snpAlleles).make());
                            }
                        }
                        refPos++;
                        alignmentPos++;
                    }
                    break;
                }
            case N:
            case H:
            case P:
            default:
                throw new GATKException("Unsupported cigar operator created during SW alignment: " + ce.getOperator());
        }
    }
    for (final VariantContext proposedEvent : proposedEvents) addVC(proposedEvent, true);
}
Also used : Cigar(htsjdk.samtools.Cigar) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) VariantContext(htsjdk.variant.variantcontext.VariantContext) GATKException(org.broadinstitute.hellbender.exceptions.GATKException) CigarElement(htsjdk.samtools.CigarElement)

Example 29 with VariantContextBuilder

use of htsjdk.variant.variantcontext.VariantContextBuilder in project gatk by broadinstitute.

the class EventMap method updateToBlockSubstitutionIfBetter.

protected boolean updateToBlockSubstitutionIfBetter(final List<VariantContext> neighbors) {
    if (neighbors.size() < MIN_NUMBER_OF_EVENTS_TO_COMBINE_INTO_BLOCK_SUBSTITUTION)
        return false;
    // TODO -- need more tests to decide if this is really so good
    final VariantContext first = neighbors.get(0);
    final int refStartOffset = first.getStart() - refLoc.getStart();
    final int refEndOffset = neighbors.get(neighbors.size() - 1).getEnd() - refLoc.getStart();
    final byte[] refBases = Arrays.copyOfRange(ref, refStartOffset, refEndOffset + 1);
    final byte[] hapBases = AlignmentUtils.getBasesCoveringRefInterval(refStartOffset, refEndOffset, haplotype.getBases(), haplotype.getAlignmentStartHapwrtRef(), haplotype.getCigar());
    final VariantContextBuilder builder = new VariantContextBuilder(first);
    builder.stop(first.getStart() + refBases.length - 1);
    builder.alleles(Arrays.asList(Allele.create(refBases, true), Allele.create(hapBases)));
    final VariantContext block = builder.make();
    // remove all merged events
    for (final VariantContext merged : neighbors) {
        if (remove(merged.getStart()) == null)
            throw new IllegalArgumentException("Expected to remove variant context from the event map but remove said there wasn't any element there: " + merged);
    }
    // note must be after we remove the previous events as the treeset only allows one key per start
    logger.info("Transforming into block substitution at " + block);
    addVC(block, false);
    return true;
}
Also used : VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) VariantContext(htsjdk.variant.variantcontext.VariantContext)

Example 30 with VariantContextBuilder

use of htsjdk.variant.variantcontext.VariantContextBuilder in project gatk by broadinstitute.

the class SVVCFWriterUnitTest method testSortVariantsByCoordinate.

@Test
public void testSortVariantsByCoordinate() {
    final String insOne = "AAA";
    new String(SVDiscoveryTestDataProvider.makeDummySequence(100, (byte) 'A'));
    final String insTwo = "AAC";
    new String(SVDiscoveryTestDataProvider.makeDummySequence(100, (byte) 'C'));
    final String contig = "21";
    final int pos = 100001;
    final int end = 100501;
    final VariantContext inversionOne = new VariantContextBuilder().chr(contig).start(pos).stop(end).alleles("G", SVConstants.DiscoveryStepConstants.VCF_ALT_ALLELE_STRING_INV).attribute(GATKSVVCFHeaderLines.INSERTED_SEQUENCE, insOne).make();
    final VariantContext inversionTwo = new VariantContextBuilder().chr(contig).start(pos).stop(end).alleles("G", SVConstants.DiscoveryStepConstants.VCF_ALT_ALLELE_STRING_INV).attribute(GATKSVVCFHeaderLines.INSERTED_SEQUENCE, insTwo).make();
    final VariantContext upstreamVariant = new VariantContextBuilder().chr(contig).start(pos - 50).stop(end).alleles("T", SVConstants.DiscoveryStepConstants.VCF_ALT_ALLELE_STRING_DUP).attribute(GATKSVVCFHeaderLines.INSERTED_SEQUENCE, insOne).make();
    final VariantContext downstreamVariant = new VariantContextBuilder().chr(contig).start(pos + 20).stop(end + 20).alleles("C", SVConstants.DiscoveryStepConstants.VCF_ALT_ALLELE_STRING_INS).attribute(GATKSVVCFHeaderLines.INSERTED_SEQUENCE, insOne).make();
    final File referenceDictionaryFile = new File(ReferenceUtils.getFastaDictionaryFileName(b37_reference_20_21));
    final SAMSequenceDictionary dictionary = ReferenceUtils.loadFastaDictionary(referenceDictionaryFile);
    final List<VariantContext> sortedVariants = SVVCFWriter.sortVariantsByCoordinate(Arrays.asList(downstreamVariant, inversionTwo, inversionOne, upstreamVariant), dictionary);
    Assert.assertEquals(sortedVariants.size(), 4);
    Assert.assertEquals(sortedVariants.get(0), upstreamVariant);
    Assert.assertEquals(sortedVariants.get(1), inversionOne);
    Assert.assertEquals(sortedVariants.get(2), inversionTwo);
    Assert.assertEquals(sortedVariants.get(3), downstreamVariant);
}
Also used : VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) VariantContext(htsjdk.variant.variantcontext.VariantContext) File(java.io.File) SAMSequenceDictionary(htsjdk.samtools.SAMSequenceDictionary) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Aggregations

VariantContextBuilder (htsjdk.variant.variantcontext.VariantContextBuilder)57 VariantContext (htsjdk.variant.variantcontext.VariantContext)40 Test (org.testng.annotations.Test)26 Allele (htsjdk.variant.variantcontext.Allele)25 BaseTest (org.broadinstitute.hellbender.utils.test.BaseTest)21 Genotype (htsjdk.variant.variantcontext.Genotype)12 ReferenceContext (org.broadinstitute.hellbender.engine.ReferenceContext)8 SimpleInterval (org.broadinstitute.hellbender.utils.SimpleInterval)7 GenotypeBuilder (htsjdk.variant.variantcontext.GenotypeBuilder)6 GenotypesContext (htsjdk.variant.variantcontext.GenotypesContext)6 File (java.io.File)6 SAMSequenceDictionary (htsjdk.samtools.SAMSequenceDictionary)5 VCFHeader (htsjdk.variant.vcf.VCFHeader)5 CommandLineProgramTest (org.broadinstitute.hellbender.CommandLineProgramTest)5 BeforeTest (org.testng.annotations.BeforeTest)5 SAMSequenceRecord (htsjdk.samtools.SAMSequenceRecord)4 VariantContextWriter (htsjdk.variant.variantcontext.writer.VariantContextWriter)4 VCFHeaderLine (htsjdk.variant.vcf.VCFHeaderLine)3 ArrayList (java.util.ArrayList)3 Collectors (java.util.stream.Collectors)3