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();
}
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();
}
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);
}
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;
}
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);
}
Aggregations