Java Code Examples for htsjdk.variant.variantcontext.VariantContext

The following examples show how to use htsjdk.variant.variantcontext.VariantContext. These examples are extracted from open source projects. You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example 1
public void add(final VariantContext vc) {
    if (lastIndel == null && vc.isIndel()) {
        buffer.add(vc);
        lastIndel = vc;
    } else if (nearby(lastIndel, vc)) {
        if (vc.isIndel()) {
            emitAllNonIndels(); // throw out {@code lastIndel} and {@code vc}
        } else {
            buffer.add(vc);
        }
    } else {
        emitAllVariants();
        buffer.add(vc);
    }

    lastIndel = vc.isIndel() ? vc : lastIndel;
}
 
Example 2
Source Project: Drop-seq   Source File: VariantContextSingletonFilterTest.java    License: MIT License 6 votes vote down vote up
@Test
public void filterOutHetSinglto32() {
 VariantContextBuilder b = new VariantContextBuilder();
 Allele a1 = Allele.create("A", true);
 Allele a2 = Allele.create("T", false);

 final List<Allele> allelesHet = new ArrayList<>(Arrays.asList(a1,a2));
 final List<Allele> allelesRef = new ArrayList<>(Arrays.asList(a1,a1));
 final List<Allele> allelesAlt = new ArrayList<>(Arrays.asList(a2,a2));

 // this does not have a het singleton because it has an alt.
 Collection<Genotype> genotypes = new ArrayList<>();
 genotypes.add(GenotypeBuilder.create("donor1", allelesRef));
 genotypes.add(GenotypeBuilder.create("donor2", allelesRef));
 genotypes.add(GenotypeBuilder.create("donor3", allelesAlt));

 VariantContext vc = b.alleles(allelesHet).start(1).stop(1).chr("1").genotypes(genotypes).make();
 Assert.assertNotNull(vc);

 Iterator<VariantContext> underlyingIterator = Collections.emptyIterator();

 VariantContextSingletonFilter f = new VariantContextSingletonFilter(underlyingIterator, true);
 boolean t1 = f.filterOut(vc);
 Assert.assertTrue(t1);

}
 
Example 3
private static void experimentalInterpretation(final JavaSparkContext ctx,
                                               final FindBreakpointEvidenceSpark.AssembledEvidenceResults assembledEvidenceResults,
                                               final SvDiscoveryInputMetaData svDiscoveryInputMetaData) {

    final JavaRDD<GATKRead> assemblyRawAlignments = getContigRawAlignments(ctx, assembledEvidenceResults, svDiscoveryInputMetaData);

    final String updatedOutputPath = svDiscoveryInputMetaData.getOutputPath() + "experimentalInterpretation_";
    svDiscoveryInputMetaData.updateOutputPath(updatedOutputPath);

    SvDiscoverFromLocalAssemblyContigAlignmentsSpark.AssemblyContigsClassifiedByAlignmentSignatures contigsByPossibleRawTypes
            = SvDiscoverFromLocalAssemblyContigAlignmentsSpark.preprocess(svDiscoveryInputMetaData, assemblyRawAlignments);

    final List<VariantContext> variants = SvDiscoverFromLocalAssemblyContigAlignmentsSpark
            .dispatchJobs(ctx, contigsByPossibleRawTypes, svDiscoveryInputMetaData, assemblyRawAlignments, true);
    contigsByPossibleRawTypes.unpersist();

    SvDiscoverFromLocalAssemblyContigAlignmentsSpark.filterAndWriteMergedVCF(updatedOutputPath, variants, svDiscoveryInputMetaData);
}
 
Example 4
@SuppressWarnings({"unchecked"})
private void assertGatherProducesCorrectVariants(GatherVcfsCloud.GatherType gatherType, File expected, List<File> inputs) throws IOException {
    final File output = createTempFile("gathered", ".vcf.gz");
    final ArgumentsBuilder args = new ArgumentsBuilder();
    inputs.forEach(args::addInput);
    args.addOutput(output)
            .add(GatherVcfsCloud.GATHER_TYPE_LONG_NAME, gatherType.toString());
    runCommandLine(args);

    try (final AbstractFeatureReader<VariantContext, ?> actualReader = AbstractFeatureReader.getFeatureReader(output.getAbsolutePath(), null, new VCFCodec(), false, Function.identity(), Function.identity());
         final AbstractFeatureReader<VariantContext, ?> expectedReader = AbstractFeatureReader.getFeatureReader(expected.getAbsolutePath(), null, new VCFCodec(), false, Function.identity(), Function.identity())) {

        final List<VariantContext> actualVariants = StreamSupport.stream(Spliterators.spliteratorUnknownSize(actualReader.iterator(),Spliterator.ORDERED), false).collect(Collectors.toList());
        final List<VariantContext> expectedVariants = StreamSupport.stream(Spliterators.spliteratorUnknownSize(expectedReader.iterator(),Spliterator.ORDERED), false).collect(Collectors.toList());
        VariantContextTestUtils.assertEqualVariants(actualVariants, expectedVariants);

        Assert.assertEquals(((VCFHeader) actualReader.getHeader()).getMetaDataInInputOrder(),
                ((VCFHeader) expectedReader.getHeader()).getMetaDataInInputOrder());
    }
}
 
Example 5
@Test
public void testMQHeadersAreUpdated() throws Exception {
    final File output = createTempFile("reblockedgvcf", ".vcf");
    final ArgumentsBuilder args = new ArgumentsBuilder();
    args.add("V", getToolTestDataDir() + "justHeader.g.vcf")
            .addOutput(output);
    runCommandLine(args);

    Pair<VCFHeader, List<VariantContext>> actual = VariantContextTestUtils.readEntireVCFIntoMemory(output.getAbsolutePath());
    VCFHeader header = actual.getLeft();
    List<VCFInfoHeaderLine> infoLines = new ArrayList<>(header.getInfoHeaderLines());
    //check all the headers in case there's one old and one updated
    for (final VCFInfoHeaderLine line : infoLines) {
        if (line.getID().equals(GATKVCFConstants.RAW_RMS_MAPPING_QUALITY_DEPRECATED)) {
            Assert.assertTrue(line.getType().equals(VCFHeaderLineType.Float));
            Assert.assertTrue(line.getDescription().contains("deprecated"));
        } else if (line.getID().equals(GATKVCFConstants.MAPPING_QUALITY_DEPTH_DEPRECATED)) {
            Assert.assertTrue(line.getDescription().contains("deprecated"));
        }
    }
}
 
Example 6
Source Project: gatk   Source File: AS_RankSumTest.java    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
private String makeReducedAnnotationString(final VariantContext vc, final Map<Allele,Double> perAltRankSumResults) {
    final StringBuilder annotationString = new StringBuilder();
    for (final Allele a : vc.getAlternateAlleles()) {
        if (annotationString.length() != 0) {
            annotationString.append(AnnotationUtils.ALLELE_SPECIFIC_REDUCED_DELIM);
        }
        if (!perAltRankSumResults.containsKey(a)) {
            logger.warn("VC allele not found in annotation alleles -- maybe there was trimming?  Allele " + a.getDisplayString() + " will be marked as missing.");
            annotationString.append(VCFConstants.MISSING_VALUE_v4); //add the missing value or the array size and indexes won't check out
        } else {
            final Double numericAlleleValue = perAltRankSumResults.get(a);
            final String perAlleleValue = numericAlleleValue != null ? String.format("%.3f", numericAlleleValue) : VCFConstants.MISSING_VALUE_v4;
            annotationString.append(perAlleleValue);
        }
    }
    return annotationString.toString();
}
 
Example 7
@VisibleForTesting
double artifactProbability(final ReferenceContext referenceContext, final VariantContext vc, final Genotype g) {
    // As of June 2018, genotype is hom ref iff we have the normal sample, but this may change in the future
    // TODO: handle MNVs
    if (g.isHomRef() || (!vc.isSNP() && !vc.isMNP()) ){
        return 0;
    } else if (!artifactPriorCollections.containsKey(g.getSampleName())) {
        return 0;
    }

    final double[] tumorLods = VariantContextGetters.getAttributeAsDoubleArray(vc, GATKVCFConstants.TUMOR_LOG_10_ODDS_KEY, () -> null, -1);
    final int indexOfMaxTumorLod = MathUtils.maxElementIndex(tumorLods);
    final Allele altAllele = vc.getAlternateAllele(indexOfMaxTumorLod);
    final byte[] altBases = altAllele.getBases();

    // for MNVs, treat each base as an independent substitution and take the maximum of all error probabilities
    return IntStream.range(0, altBases.length).mapToDouble(n -> {
        final Nucleotide altBase = Nucleotide.valueOf(new String(new byte[] {altBases[n]}));

        return artifactProbability(referenceContext, vc.getStart() + n, g, indexOfMaxTumorLod, altBase);
    }).max().orElse(0.0);

}
 
Example 8
@Test()
public void testTetraploidRun() throws IOException {
    final File output = createTempFile("combinegvcfs", ".vcf");

    final ArgumentsBuilder args = new ArgumentsBuilder();
    args.addReference(new File(b37_reference_20_21))
            .addOutput(output);
    args.add("variant",getToolTestDataDir()+"tetraploid-gvcf-1.vcf");
    args.add("variant",getToolTestDataDir()+"tetraploid-gvcf-2.vcf");
    args.add("variant",getToolTestDataDir()+"tetraploid-gvcf-3.vcf");
    args.add("intervals", getToolTestDataDir() + "tetraploid-gvcfs.interval_list");

    runCommandLine(args);

    final List<VariantContext> expectedVC = getVariantContexts(getTestFile("tetraploidRun.GATK3.g.vcf"));
    final List<VariantContext> actualVC = getVariantContexts(output);
    final VCFHeader header = getHeaderFromFile(output);
    assertForEachElementInLists(actualVC, expectedVC, (a, e) -> VariantContextTestUtils.assertVariantContextsAreEqualAlleleOrderIndependent(a, e, Arrays.asList(), Collections.emptyList(), header));

}
 
Example 9
@Test(dataProvider = "MaxMinIndelSize")
public void maxIndelSizeTest(final int size, final int otherSize, final int max, final int min, final String op) {

    final byte[] largerAllele = Utils.dupBytes((byte) 'A', size+1);
    final byte[] smallerAllele = Utils.dupBytes((byte) 'A', 1);

    final List<Allele> alleles = new ArrayList<>(2);
    final Allele ref = Allele.create(op.equals("I") ? smallerAllele : largerAllele, true);
    final Allele alt = Allele.create(op.equals("D") ? smallerAllele : largerAllele, false);
    alleles.add(ref);
    alleles.add(alt);

    final VariantContext vc = new VariantContextBuilder("test", "1", 10, 10 + ref.length() - 1, alleles).make();

    final boolean hasIndelTooLargeOrSmall = SelectVariants.containsIndelLargerOrSmallerThan(vc, max, min);
    Assert.assertEquals(hasIndelTooLargeOrSmall, size > max || size < min);
}
 
Example 10
Source Project: gatk   Source File: IndelSize.java    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
public List<Object> getRelevantStates(ReferenceContext referenceContext, ReadsContext readsContext, FeatureContext featureContext, VariantContext comp, String compName, VariantContext eval, String evalName, String sampleName, String FamilyName) {
    if (eval != null && eval.isIndel() && eval.isBiallelic()) {
        try {
            int eventLength = 0;
            if ( eval.isSimpleInsertion() ) {
                eventLength = eval.getAlternateAllele(0).length();
            } else if ( eval.isSimpleDeletion() ) {
                eventLength = -eval.getReference().length();
            }

            if (eventLength > MAX_INDEL_SIZE)
                eventLength = MAX_INDEL_SIZE;
            else if (eventLength < -MAX_INDEL_SIZE)
                eventLength = -MAX_INDEL_SIZE;

            return Collections.singletonList((Object)eventLength);
        } catch (Exception e) {
            return Collections.emptyList();
        }
    }

    return Collections.emptyList();
}
 
Example 11
Source Project: hmftools   Source File: MNVMerger.java    License: GNU General Public License v3.0 6 votes vote down vote up
@NotNull
private static List<Allele> createMnvAlleles(@NotNull final List<VariantContext> variants,
        @NotNull final Map<Integer, Character> gapReads) {
    final StringBuilder refBases = new StringBuilder();
    final StringBuilder altBases = new StringBuilder();
    int currentPosition = variants.get(0).getStart();
    for (final VariantContext variant : variants) {
        while (currentPosition != variant.getStart()) {
            final Character gapRead = gapReads.get(currentPosition);
            refBases.append(gapRead);
            altBases.append(gapRead);
            currentPosition++;
        }
        refBases.append(variant.getReference().getBaseString());
        altBases.append(variant.getAlternateAllele(0).getBaseString());
        currentPosition = variant.getStart() + variant.getReference().length();
    }
    final Allele ref = Allele.create(refBases.toString(), true);
    final Allele alt = Allele.create(altBases.toString(), false);
    return Lists.newArrayList(ref, alt);
}
 
Example 12
@NotNull
static Pair<Integer, String> relativePositionAndRef(@NotNull final IndexedFastaSequenceFile reference, @NotNull final VariantContext variant) {
    final int refLength = variant.getReference().getBaseString().length();
    @Nullable
    final SAMSequenceRecord samSequenceRecord = reference.getSequenceDictionary().getSequence(variant.getContig());
    if (samSequenceRecord == null) {
        LOGGER.warn("Unable to locate contig {} in ref genome", variant.getContig());
        return new Pair<>(-1, Strings.EMPTY);
    }

    final int chromosomeLength = samSequenceRecord.getSequenceLength();
    long positionBeforeEvent = variant.getStart();

    long start = Math.max(positionBeforeEvent - 100, 1);
    long end = Math.min(positionBeforeEvent + refLength + 100 - 1, chromosomeLength - 1);
    int relativePosition = (int) (positionBeforeEvent - start);
    final String sequence;
    if (start < chromosomeLength && end < chromosomeLength) {
        sequence = reference.getSubsequenceAt(variant.getContig(), start, end).getBaseString();
    } else {
        sequence = Strings.EMPTY;
        LOGGER.warn("Requested base sequence outside of chromosome region!");
    }
    return new Pair<>(relativePosition, sequence);
}
 
Example 13
Source Project: gatk   Source File: CNVInputReader.java    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
/**
 * Loads an external cnv call list and returns the results in an SVIntervalTree. NB: the contig indices in the SVIntervalTree
 * are based on the sequence indices in the SAM header, _NOT_ the ReadMetadata (which we might not have access to at this
 * time).
 */
public static SVIntervalTree<VariantContext> loadCNVCalls(final String cnvCallsFile,
                                                           final SAMFileHeader headerForReads) {
    Utils.validate(cnvCallsFile != null, "Can't load null CNV calls file");
    try ( final FeatureDataSource<VariantContext> dataSource = new FeatureDataSource<>(cnvCallsFile, null, 0, null) ) {

        final String sampleId = SVUtils.getSampleId(headerForReads);
        validateCNVcallDataSource(headerForReads, sampleId, dataSource);

        final SVIntervalTree<VariantContext> cnvCallTree = new SVIntervalTree<>();
        Utils.stream(dataSource.iterator())
                .map(vc -> new VariantContextBuilder(vc).genotypes(vc.getGenotype(sampleId)).make()) // forces a decode of the genotype for serialization purposes
                .map(vc -> new Tuple2<>(new SVInterval(headerForReads.getSequenceIndex(vc.getContig()), vc.getStart(), vc.getEnd()),vc))
                .forEach(pv -> cnvCallTree.put(pv._1(), pv._2()));
        return cnvCallTree;
    }
}
 
Example 14
@VisibleForTesting
static VariantContextBuilder annotateWithExternalCNVCalls(final String recordContig, final int pos, final int end,
                                                          final VariantContextBuilder inputBuilder,
                                                          final Broadcast<SAMSequenceDictionary> broadcastSequenceDictionary,
                                                          final Broadcast<SVIntervalTree<VariantContext>> broadcastCNVCalls,
                                                          final String sampleId) {
    if (broadcastCNVCalls == null)
        return inputBuilder;
    final SVInterval variantInterval = new SVInterval(broadcastSequenceDictionary.getValue().getSequenceIndex(recordContig), pos, end);
    final SVIntervalTree<VariantContext> cnvCallTree = broadcastCNVCalls.getValue();
    final String cnvCallAnnotation =
            Utils.stream(cnvCallTree.overlappers(variantInterval))
                    .map(overlapper -> formatExternalCNVCallAnnotation(overlapper.getValue(), sampleId))
                    .collect(Collectors.joining(VCFConstants.INFO_FIELD_ARRAY_SEPARATOR));
    if (!cnvCallAnnotation.isEmpty()) {
        return inputBuilder.attribute(GATKSVVCFConstants.EXTERNAL_CNV_CALLS, cnvCallAnnotation);
    } else
        return inputBuilder;
}
 
Example 15
Source Project: picard   Source File: MakeVcfSampleNameMap.java    License: MIT License 5 votes vote down vote up
private static AbstractFeatureReader<VariantContext, LineIterator> getReaderFromPath(final Path variantPath) {
    final String variantURI = variantPath.toAbsolutePath().toUri().toString();
    try {
        return AbstractFeatureReader.getFeatureReader(variantURI, null, new VCFCodec(),
                false, Function.identity(), Function.identity());
    } catch (final TribbleException e) {
        throw new PicardException("Failed to create reader from " + variantURI, e);
    }
}
 
Example 16
Source Project: gatk   Source File: ApplyVQSR.java    License: BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
private VariantContext getMatchingRecalVC(final VariantContext target, final List<VariantContext> recalVCs, final Allele allele) {
    for( final VariantContext recalVC : recalVCs ) {
        if ( target.getEnd() == recalVC.getEnd() ) {
            if (!useASannotations)
                return recalVC;
            else if (allele.equals(recalVC.getAlternateAllele(0)))
                return recalVC;
        }
    }
    return null;
}
 
Example 17
/**
 * Get the {@link org.broadinstitute.hellbender.tools.funcotator.dataSources.gencode.GencodeFuncotation.VariantClassification} for a given {@code variant}/{@code allele} in a coding region.
 * @param variant The {@link VariantContext} to classify.
 * @param altAllele The {@link Allele} of the given {@code variant} to classify.
 * @param variantType The {@link org.broadinstitute.hellbender.tools.funcotator.dataSources.gencode.GencodeFuncotation.VariantType} of the given {@code variant}.
 * @param sequenceComparison The {@link org.broadinstitute.hellbender.tools.funcotator.SequenceComparison} for the given {@code variant}.  Must have a non-null proteinChangeInfo Object.
 * @return A {@link org.broadinstitute.hellbender.tools.funcotator.dataSources.gencode.GencodeFuncotation.VariantClassification} based on the given {@code allele}, {@code variant}, {@code exon}, and {@code sequenceComparison}.
 */
private static GencodeFuncotation.VariantClassification getVariantClassificationForCodingRegion(final VariantContext variant,
                                                                                         final Allele altAllele,
                                                                                         final GencodeFuncotation.VariantType variantType,
                                                                                         final SequenceComparison sequenceComparison) {
    final GencodeFuncotation.VariantClassification varClass;

    if (variantType == GencodeFuncotation.VariantType.INS) {
        if ( GATKVariantContextUtils.isFrameshift(variant.getReference(), altAllele)) {
            varClass = GencodeFuncotation.VariantClassification.FRAME_SHIFT_INS;
        }
        else {
            varClass = GencodeFuncotation.VariantClassification.IN_FRAME_INS;
        }
    }
    else if (variantType == GencodeFuncotation.VariantType.DEL) {
        if (GATKVariantContextUtils.isFrameshift(variant.getReference(), altAllele)) {
            varClass = GencodeFuncotation.VariantClassification.FRAME_SHIFT_DEL;
        }
        else {
            varClass = GencodeFuncotation.VariantClassification.IN_FRAME_DEL;
        }
    }
    else {
        // This is a SNP/MNP
        // We just check to see what the protein change is to check for MISSENSE/NONSENSE/SILENT:
        varClass = getVarClassFromEqualLengthCodingRegions( sequenceComparison );
    }

    return varClass;
}
 
Example 18
@Override
public void apply(final VariantContext variant,
                  final ReadsContext readsContext,
                  final ReferenceContext referenceContext,
                  final FeatureContext featureContext) {

    final Collection<VariantContext> otherVCs = featureContext.getValues(supportVariants);

    //If no resource contains a matching variant, then add numRefIfMissing as a pseudocount to the priors
    List<VariantContext> resourcesWithMatchingStarts = otherVCs.stream()
            .filter(vc -> variant.getStart() == vc.getStart()).collect(Collectors.toList());

    //do family priors first (if applicable)
    final VariantContextBuilder builder = new VariantContextBuilder(variant);
    //only compute family priors for biallelelic sites
    if (!skipFamilyPriors && variant.isBiallelic()){
        final GenotypesContext gc = famUtils.calculatePosteriorGLs(variant);
        builder.genotypes(gc);
    }
    VariantContextUtils.calculateChromosomeCounts(builder, false);
    final VariantContext vc_familyPriors = builder.make();

    final VariantContext vc_bothPriors;
    if (!skipPopulationPriors) {
        vc_bothPriors = PosteriorProbabilitiesUtils.calculatePosteriorProbs(vc_familyPriors, resourcesWithMatchingStarts,
                resourcesWithMatchingStarts.isEmpty() ? numRefIfMissing : 0, options);
    } else {
        vc_bothPriors = vc_familyPriors;
    }
    vcfWriter.add(vc_bothPriors);
}
 
Example 19
@DataProvider(name = "forTestDeletionConsistencyCheck")
private Object[][] forTestDeletionConsistencyCheck() {
    final List<Object[]> data = new ArrayList<>(20);

    VariantContext del = makeDeletion(new SimpleInterval("chr6:857169-857562"), Allele.create("G", true)).make();
    HashSet<SimpleInterval> missingSegments = Sets.newHashSet(new SimpleInterval("chr6:857170-857564"), new SimpleInterval("chr6:857767-857852"));

    data.add(new Object[]{del, Collections.emptySet(), false});
    data.add(new Object[]{del, missingSegments, true});

    data.add(new Object[]{del, Collections.singleton(new SimpleInterval("chr6:857767-857852")), false});

    return data.toArray(new Object[data.size()][]);
}
 
Example 20
private static Tuple2<GATKRead, Iterable<GATKVariant>> getVariantsOverlappingRead(final GATKRead read, final List<FeatureDataSource<VariantContext>> variantSources) {
    if (SimpleInterval.isValid(read.getContig(), read.getStart(), read.getEnd())) {
        return new Tuple2<>(read, getVariantsOverlappingInterval(variantSources, new SimpleInterval(read)));
    } else {
        //Sometimes we have reads that do not form valid intervals (reads that do not consume any ref bases, eg CIGAR 61S90I
        //In those cases, we'll just say that nothing overlaps the read
        return new Tuple2<>(read, Collections.emptyList());
    }
}
 
Example 21
Source Project: Hadoop-BAM   Source File: BCFRecordWriter.java    License: MIT License 5 votes vote down vote up
protected void writeRecord(VariantContext vc) {
	final GenotypesContext gc = vc.getGenotypes();
	if (gc instanceof LazyParsingGenotypesContext)
		((LazyParsingGenotypesContext)gc).getParser().setHeaderDataCache(
			gc instanceof LazyVCFGenotypesContext ? vcfHeaderDataCache
			                                      : bcfHeaderDataCache);

	writer.add(vc);
}
 
Example 22
@Override
protected void nthPassApply(final VariantContext variant,
                            final ReadsContext readsContext,
                            final ReferenceContext referenceContext,
                            final FeatureContext featureContext,
                            final int n) {
    ParamUtils.isPositiveOrZero(n, "Passes must start at the 0th pass.");
    if (n <= NUMBER_OF_LEARNING_PASSES) {
        filteringEngine.accumulateData(variant, referenceContext);
    } else if (n == NUMBER_OF_LEARNING_PASSES + 1) {
        vcfWriter.add(filteringEngine.applyFiltersAndAccumulateOutputStats(variant, referenceContext));
    } else {
        throw new GATKException.ShouldNeverReachHereException("This walker should never reach (zero-indexed) pass " + n);
    }
}
 
Example 23
@DataProvider(name = "DetectFeatureTypeFromFeatureInputFieldTestData")
public Object[][] getDetectFeatureTypeFromFeatureInputFieldTestData() {
    return new Object[][] {
            { "variantContextFeatureInput", VariantContext.class },
            { "variantContextListFeatureInput", VariantContext.class },
            { "bedFeatureInput", BEDFeature.class },
            { "bedListFeatureInput", BEDFeature.class }
    };
}
 
Example 24
/**
 * Creates a result indicating that there was no trimming to be done.
 */
protected static Result noTrimming(final boolean emitReferenceConfidence,
                                   final AssemblyRegion targetRegion, final int padding,
                                   final int usableExtension,final List<VariantContext> events) {
    final SimpleInterval targetRegionLoc = targetRegion.getSpan();
    final Result result = new Result(emitReferenceConfidence,false,targetRegion,padding,usableExtension,events,Pair.of(null, null),
            targetRegionLoc,targetRegionLoc,targetRegionLoc,targetRegionLoc);
    result.callableRegion = targetRegion;
    return result;
}
 
Example 25
Source Project: genomewarp   Source File: VcfToVariantTest.java    License: Apache License 2.0 5 votes vote down vote up
@Test
public void testGetFilter() throws Exception {
  File vcfFile =
      new File(VcfToVariant.class.getClassLoader().getResource(VALID_VCF_4_1).getFile());
  VCFFileReader vcfReader = new VCFFileReader(vcfFile, false);
  int currVariant = 0;

  for (final VariantContext vc : vcfReader) {
    List<String> filters = VcfToVariant.getFilter(vc);
    assertEquals(filters, TRUTH.get(currVariant).getFilterList());
    currVariant++;
  }
}
 
Example 26
/**
 * Test that we can annotate a b37 (here it is clinvar) datasource when GENCODE is using "chr*" and the datasource is not.
 */
@Test
public void testCanAnnotateMixedContigHg19Clinvar() {
    final FuncotatorArgumentDefinitions.OutputFormatType outputFormatType = FuncotatorArgumentDefinitions.OutputFormatType.VCF;
    final File outputFile = getOutputFile(outputFormatType);

    final ArgumentsBuilder arguments = createBaselineArgumentsForFuncotator(
            PIK3CA_VCF_HG19,
            outputFile,
            b37Chr3Ref,
            DS_PIK3CA_DIR,
            FuncotatorTestConstants.REFERENCE_VERSION_HG19,
            outputFormatType,
            false);

    // We need this argument since we are testing on a subset of b37
    arguments.add(FuncotatorArgumentDefinitions.FORCE_B37_TO_HG19_REFERENCE_CONTIG_CONVERSION, true);

    runCommandLine(arguments);

    final List<VariantContext> variantContexts = new ArrayList<>();
    final FeatureDataSource<VariantContext> featureDataSource = new FeatureDataSource<>(outputFile);
    for (final VariantContext vc : featureDataSource) {
        variantContexts.add(vc);
    }

    final int NUM_VARIANTS = 21;
    final int NUM_CLINVAR_HITS = 4;
    Assert.assertEquals(variantContexts.size(), NUM_VARIANTS, "Found unexpected number of variants!");

    // Look for "MedGen" to know that we have a clinvar hit.
    Assert.assertEquals(variantContexts.stream()
            .filter(vc -> StringUtils.contains(vc.getAttributeAsString(VcfOutputRenderer.FUNCOTATOR_VCF_FIELD_NAME, ""), "MedGen"))
            .count(), NUM_CLINVAR_HITS, "Found unexpected number of ClinVar hits!");
}
 
Example 27
Source Project: hmftools   Source File: GermlineFilters.java    License: GNU General Public License v3.0 5 votes vote down vote up
public static boolean hasStrandBias(final StructuralVariant sv, final VariantContext variant)
{
    // Filter deletion or duplication breakpoints under 1000bp with a split read strand bias of 0.95 or greater
    if(!isShortDelDup(sv))
        return false;

    return variant.getCommonInfo().getAttributeAsDouble(SB, 0.0) >= MAX_ALLOWABLE_SHORT_EVENT_STRAND_BIAS;
}
 
Example 28
/**
 * 40-'C' + 'ATCG' + 34 bases of unique sequence + 'ATCG' + 40-'T' is shrunk to 40-'C' + 'ATCG' + 40-'T' (forward strand representation)
 * Return a list of two entries for positive and reverse strand representations.
 */
private static List<TestDataForSimpleSV>
forDeletionWithHomology() {

    final List<TestDataForSimpleSV> result = new ArrayList<>();

    // simple deletion with homology '+' strand representation
    final String leftRefFlank = TestUtilsForAssemblyBasedSVDiscovery.makeDummySequence('C', 40);
    final String rightRefFlank = TestUtilsForAssemblyBasedSVDiscovery.makeDummySequence('T', 40);
    final String homology = "ATCG";
    byte[] contigSeq = (leftRefFlank + homology + rightRefFlank).getBytes();
    String contigName = "simple_del_with_hom_+";

    final SimpleInterval expectedLeftBreakpoint = new SimpleInterval("21:17000040-17000040");
    final SimpleInterval expectedRightBreakpoint = new SimpleInterval("21:17000078-17000078");
    final BreakpointComplications expectedBreakpointComplications = new BreakpointComplications.SimpleInsDelOrReplacementBreakpointComplications(homology, "");
    final byte[] expectedAltSeq = EMPTY_BYTE_ARRAY;
    final NovelAdjacencyAndAltHaplotype expectedNovelAdjacencyAndAltHaplotype = new NovelAdjacencyAndAltHaplotype(expectedLeftBreakpoint, expectedRightBreakpoint, NO_SWITCH, expectedBreakpointComplications, SIMPLE_DEL, expectedAltSeq);
    AlignmentInterval region1 = new AlignmentInterval(new SimpleInterval("21", 17000001, 17000044), 1 ,44, TextCigarCodec.decode("44M40S"), true, 60, 0, 100, ContigAlignmentsModifier.AlnModType.NONE);
    AlignmentInterval region2 = new AlignmentInterval(new SimpleInterval("21", 17000079, 17000122), 41 ,84, TextCigarCodec.decode("40S44M"), true, 60, 0, 100, ContigAlignmentsModifier.AlnModType.NONE);
    SimpleChimera expectedSimpleChimera = new SimpleChimera(contigName, region1, region2, NO_SWITCH, true, Collections.emptyList(), NO_GOOD_MAPPING_TO_NON_CANONICAL_CHROMOSOME);
    DistancesBetweenAlignmentsOnRefAndOnRead expectedDistances = new DistancesBetweenAlignmentsOnRefAndOnRead(34, -4, 17000044, 17000079, 44, 41);
    final List<SvType> expectedSVTypes = Collections.singletonList(makeDeletionType(new SimpleInterval("21:17000040-17000078"), Allele.create("G", true),false));
    final List<VariantContext> expectedVariants = Collections.singletonList(
            addStandardAttributes(makeDeletion(new SimpleInterval("21:17000040-17000078"), Allele.create("G", true), false),
                    40, contigName, SimpleSVType.SupportedType.DEL.name(), 17000078, -38, "", homology, "").make());
    result.add(new TestDataForSimpleSV(region1, region2, contigName, contigSeq, false, expectedSimpleChimera, expectedNovelAdjacencyAndAltHaplotype, expectedSVTypes, expectedVariants, expectedDistances, BreakpointsInference.SimpleInsertionDeletionBreakpointsInference.class));

    // simple deletion with homology '-' strand representation
    contigSeq = (SequenceUtil.reverseComplement(rightRefFlank) + SequenceUtil.reverseComplement(homology) + SequenceUtil.reverseComplement(leftRefFlank)).getBytes();
    contigName = "simple_del_with_hom_-";
    region1 = new AlignmentInterval(new SimpleInterval("21", 17000079, 17000122), 1 ,44, TextCigarCodec.decode("44M40S"), false, 60, 0, 100, ContigAlignmentsModifier.AlnModType.NONE);
    region2 = new AlignmentInterval(new SimpleInterval("21", 17000001, 17000044), 41 ,84, TextCigarCodec.decode("40S44M"), false, 60, 0, 100, ContigAlignmentsModifier.AlnModType.NONE);
    expectedSimpleChimera = new SimpleChimera(contigName, region1, region2, NO_SWITCH, false, Collections.emptyList(), NO_GOOD_MAPPING_TO_NON_CANONICAL_CHROMOSOME);
    expectedDistances = new DistancesBetweenAlignmentsOnRefAndOnRead(34, -4, 17000044, 17000079, 44, 41);
    result.add(new TestDataForSimpleSV(region1, region2, contigName, contigSeq, true, expectedSimpleChimera, expectedNovelAdjacencyAndAltHaplotype, expectedSVTypes, expectedVariants, expectedDistances,  BreakpointsInference.SimpleInsertionDeletionBreakpointsInference.class));

    return result;
}
 
Example 29
@Override
public boolean test(final VariantContext variant) {
    final boolean accept = delegateFilter.test(variant);
    if (!accept) {
        filteredCount++;
    }
    return accept;
}
 
Example 30
private void assertInfoFieldsAreClose(File actualVcf, File expectedVcf, String infoKey){
    Iterator<VariantContext> expectedVi = VariantContextTestUtils.streamVcf(expectedVcf).collect(Collectors.toList()).iterator();
    Iterator<VariantContext> actualVi = VariantContextTestUtils.streamVcf(actualVcf).collect(Collectors.toList()).iterator();
    while (expectedVi.hasNext() && actualVi.hasNext()) {
        VariantContext expectedVc = expectedVi.next();
        VariantContext actualVc = actualVi.next();
        double expectedScore = expectedVc.getAttributeAsDouble(infoKey, 0.0); // Different defaults trigger failures on missing scores
        double actualScore = actualVc.getAttributeAsDouble(infoKey, EPSILON+1.0);
        double diff = Math.abs(expectedScore-actualScore);
        Assert.assertTrue(diff < EPSILON);
        VariantContextTestUtils.assertVariantContextsAreEqual(actualVc, expectedVc, Collections.singletonList(infoKey), Collections.emptyList());
    }
    Assert.assertTrue(!expectedVi.hasNext() && !actualVi.hasNext());
}