Java Code Examples for htsjdk.variant.variantcontext.VariantContext#isFiltered()

private SiteStatus calcSiteStatus(VariantContext vc) {
    if ( vc == null ) return SiteStatus.NO_CALL;
    if ( vc.isFiltered() ) return SiteStatus.FILTERED;
    if ( vc.isMonomorphicInSamples() ) return SiteStatus.MONO;
    if ( vc.hasGenotypes() ) return SiteStatus.POLY;  // must be polymorphic if isMonomorphicInSamples was false and there are genotypes

    if ( vc.hasAttribute(VCFConstants.ALLELE_COUNT_KEY) ) {
        int ac = 0;
        if ( vc.getNAlleles() > 2 ) {
            return SiteStatus.POLY;
        }
        else
            ac = vc.getAttributeAsInt(VCFConstants.ALLELE_COUNT_KEY, 0);
        return ac > 0 ? SiteStatus.POLY : SiteStatus.MONO;
    } else {
        return TREAT_ALL_SITES_IN_EVAL_VCF_AS_CALLED ? SiteStatus.POLY : SiteStatus.NO_CALL; // we can't figure out what to do
    }
}
 

@Override
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){
        return Collections.emptyList();
    }

    ArrayList<Object> relevantStates = new ArrayList<Object>();
    if (eval.isFiltered()){
        relevantStates.addAll(eval.getFilters());
    }
    else {
        relevantStates.add("PASS");
    }

    return relevantStates;
}
 

@Test(dataProvider = "goodInputVcfs")
public void testJavaScript(final File input) throws Exception {
    final File out = VcfTestUtils.createTemporaryIndexedFile("filterVcfTestJS.", ".vcf");
    final FilterVcf filterer = new FilterVcf();
    filterer.INPUT = input;
    filterer.OUTPUT = out;
    filterer.JAVASCRIPT_FILE = quickJavascriptFilter("variant.getStart()%5 != 0");

    final int retval = filterer.doWork();
    Assert.assertEquals(retval, 0);

    //count the number of reads
    final int expectedNumber = 4;
    int count = 0;
    VCFFileReader in = new VCFFileReader(filterer.OUTPUT, false);
    CloseableIterator<VariantContext> iter = in.iterator();
    while (iter.hasNext()) {
        final VariantContext ctx = iter.next();
        count += (ctx.isFiltered() ? 1 : 0);
    }
    iter.close();
    in.close();
    Assert.assertEquals(count, expectedNumber);
}
 

private void updateDetailMetric(final VariantCallingDetailMetrics metric,
                                final Genotype genotype,
                                final VariantContext vc,
                                final boolean hasSingletonSample) {
    updateSummaryMetric(metric, genotype, vc, hasSingletonSample);

    if (genotype != null && !vc.isFiltered()) {
        if (genotype.getGQ() == 0) {
            ++metric.TOTAL_GQ0_VARIANTS;
        }
        if (genotype.isHet()) {
            ++metric.numHets;
        } else if (genotype.isHomVar()) {
            ++metric.numHomVar;
        }
    }
}
 

/**
 * Tests that all records get PASS set as their filter when extreme values are used for filtering.
 */
@Test(dataProvider = "goodInputVcfs")
public void testNoFiltering(final File input) throws Exception {
    final File out = testFiltering(input, ".vcf.gz", 0, 0, 0, Double.MAX_VALUE);
    final VCFFileReader in = new VCFFileReader(out, false);
    for (final VariantContext ctx : in) {
        if (!ctx.filtersWereApplied() || ctx.isFiltered()) {
            Assert.fail("Context should not have been filtered: " + ctx.toString());
        }
    }
    in.close();
}
 

/**
 * Returns the ID fields (separated by semicolons) of all VariantContexts in rsIDSourceVCs that have the same reference and alternate alleles
 * as vcToAnnotate, after splitting to bi-allelics and trimming both the rsIDSourceVCs and vcToAnnotate
 *
 * Doesn't require vcToAnnotate to be a complete match, so
 *
 * A/C/G in VC in rsIDSourceVCs
 *
 * would match the a VC with A/C but not A/T.  Also we don't require all alleles to match
 * so we would also match A/C/T to A/C/G.
 *
 * Will only match rsIDSourceVCs that aren't failing filters.
 *
 * @param rsIDSourceVCs a non-null list of potential overlaps that start at vcToAnnotate
 * @param vcToAnnotate a non-null VariantContext to annotate
 * @return a String to use for the rsID from rsIDSourceVCs if any match, or null if none match
 */
private static String getRsID(final List<VariantContext> rsIDSourceVCs, final VariantContext vcToAnnotate) {
    Utils.nonNull(rsIDSourceVCs, "rsIDSourceVCs cannot be null");
    Utils.nonNull(vcToAnnotate, "vcToAnnotate cannot be null");
    final List<String> rsids = new ArrayList<>();

    final List<VariantContext> vcAnnotateList = GATKVariantContextUtils.splitVariantContextToBiallelics(vcToAnnotate, true,
            GenotypeAssignmentMethod.SET_TO_NO_CALL_NO_ANNOTATIONS, true);

    for ( final VariantContext vcCompSource : rsIDSourceVCs ) {
        if ( vcCompSource.isFiltered() ) {
            continue; // don't process any failed VCs
        }

        if (!vcCompSource.getContig().equals(vcToAnnotate.getContig())) {
            throw new IllegalArgumentException("source rsID VariantContext " + vcCompSource + " is not on same chromosome as vcToAnnotate " + vcToAnnotate);
        }

        final List<VariantContext> vcCompList = GATKVariantContextUtils.splitVariantContextToBiallelics(vcCompSource, true,
                GenotypeAssignmentMethod.SET_TO_NO_CALL_NO_ANNOTATIONS, true);
        boolean addThisID = false;
        for (final VariantContext vcComp : vcCompList) {
            for (final VariantContext vcToAnnotateBi : vcAnnotateList) {
                if (vcComp.getStart() == vcToAnnotateBi.getStart() && vcToAnnotateBi.getReference().equals(vcComp.getReference()) && vcComp.getAlternateAlleles().equals(vcToAnnotateBi.getAlternateAlleles())) {
                    addThisID = true;
                    break;
                }
            }

            if (addThisID) {
                rsids.add(vcCompSource.getID());
                break;
            }
        }
    }

    return rsids.isEmpty()? null : String.join(";", rsids);
}
 

@Override
public void apply(final VariantContext variant, final ReadsContext readsContext, final ReferenceContext referenceContext, final FeatureContext featureContext) {
    if (variant.isFiltered()) {
        return;
    }

    if (variant.isIndel()) {
        indelFalsePositiveCount++;
    } else {
        snpFalsePositiveCount++;
    }
}
 

@Override
public void apply(final VariantContext variant, final ReadsContext readsContext, final ReferenceContext referenceContext, final FeatureContext featureContext) {
    if (variant.isFiltered()) {
        return;
    }

    if (variant.isIndel()) {
        indelFalsePositiveCount++;
    } else {
        snpFalsePositiveCount++;
    }
}
 

@VisibleForTesting
static List<String> getFilter(VariantContext vc) {
  if (vc.isFiltered()) {
    return (List<String>) ParsingUtils.sortList(vc.getFilters());
  }
  if (vc.filtersWereApplied()) {
    return Arrays.asList(VCFConstants.PASSES_FILTERS_v4);
  }

  return Arrays.asList(VCFConstants.UNFILTERED);
}
 

private static void attachFilter(@NotNull final ImmutableSomaticVariantImpl.Builder builder, @NotNull VariantContext context) {
    if (context.isFiltered()) {
        StringJoiner joiner = new StringJoiner(";");
        context.getFilters().forEach(joiner::add);
        builder.filter(joiner.toString());
    } else {
        builder.filter(PASS_FILTER);
    }
}
 

private void updateDetailMetric(final CollectArraysVariantCallingMetrics.ArraysVariantCallingDetailMetrics metric,
                                final Genotype genotype,
                                final VariantContext vc,
                                final boolean hasSingletonSample) {
    metric.NUM_ASSAYS++;
    if (!vc.isFiltered() || vc.getCommonInfo().getFilters().contains(InfiniumVcfFields.DUPE)) {
        metric.NUM_NON_FILTERED_ASSAYS++;
        if (genotype.isCalled()) {
            metric.NUM_CALLS++;
            String gtA = (String) genotype.getExtendedAttribute(InfiniumVcfFields.GTA, genotype.getGenotypeString());
            if (!gtA.equals(VCFConstants.EMPTY_GENOTYPE)) {
                metric.NUM_AUTOCALL_CALLS++;
            }
        } else {
            metric.NUM_NO_CALLS++;
        }

        if (vc.isSNP()) {
            // Biallelic SNPs
            final boolean isInDbSnp = dbsnp.snps.isDbSnpSite(vc.getContig(), vc.getStart());

            metric.NUM_SNPS++;

            if (isInDbSnp) {
                metric.NUM_IN_DB_SNP++;
            }
        } else if (vc.isIndel()) {
            metric.NUM_INDELS++;
        }

        if (hasSingletonSample) {
            metric.NUM_SINGLETONS++;
        }

        if (genotype.isHet()) {
            metric.numHets++;
        } else if (genotype.isHomVar()) {
            metric.numHomVar++;
        }
    } else {
        metric.NUM_FILTERED_ASSAYS++;
        if (vc.getCommonInfo().getFilters().contains(InfiniumVcfFields.ZEROED_OUT_ASSAY)) {
            // A "zeroed-out SNP".  Marked as unusable/uncallable
            metric.NUM_ZEROED_OUT_ASSAYS++;
        }
    }
}
 

public IntervalList processData(final File vcfFile, final File sdFile, final Set<String> sample, int GQThreshold, final boolean hetSNPsOnly) {

		final VCFFileReader reader = new VCFFileReader(vcfFile, false);
		if (!VCFUtils.GQInHeader(reader)) {
			GQThreshold=-1;
			log.info("Genotype Quality [GQ] not found in header.  Disabling GQ_THRESHOLD parameter");
		}
		
		final VCFHeader inputVcfHeader = new VCFHeader(reader.getFileHeader().getMetaDataInInputOrder());
		SAMSequenceDictionary sequenceDictionary = inputVcfHeader.getSequenceDictionary();
		Set<String> sampleListFinal = sample;
		if (sample==null || sample.isEmpty()) {
			ArrayList<String> s = reader.getFileHeader().getSampleNamesInOrder();
			sampleListFinal=new TreeSet<String>(s);
		}

		if (sdFile != null)
			sequenceDictionary = getSequenceDictionary(sdFile);

		final ProgressLogger progress = new ProgressLogger(this.log, 500000);

		final SAMFileHeader samHeader = new SAMFileHeader();
		samHeader.setSequenceDictionary(sequenceDictionary);
		IntervalList result = new IntervalList(samHeader);

		// Go through the input, find sites we want to keep.
		final PeekableIterator<VariantContext> iterator = new PeekableIterator<>(reader.iterator());

		validateRequestedSamples (iterator, sampleListFinal);

		while (iterator.hasNext()) {
			final VariantContext site = iterator.next();
			progress.record(site.getContig(), site.getStart());
			// for now drop any filtered site.
			if (site.isFiltered())
				continue;
			// move onto the next record if the site is not a SNP or the samples aren't all heterozygous.
			if (!site.isSNP())
				continue;
			if (!sitePassesFilters(site, sampleListFinal, GQThreshold, hetSNPsOnly))
				continue;
			Interval varInt = new Interval(site.getContig(), site.getStart(),
					site.getEnd(), true, site.getID());

			// final Interval site = findHeterozygousSites(full, SAMPLE);
			result.add(varInt);

		}

		CloserUtil.close(iterator);
		CloserUtil.close(reader);
		return (result);
	}
 

private TruthVersusEval nextEvalOnlyVariant() {
    final VariantContext evalVariant = evalIterator.next();
    return evalVariant.isFiltered() ? TruthVersusEval.filteredTrueNegative(evalVariant) : TruthVersusEval.falsePositive(evalVariant);
}
 

@Override public boolean test(final VariantContext variant) {
    return !variant.isFiltered();
}
 

@Override
public void apply(final VariantContext vc, final ReadsContext readsContext, final ReferenceContext ref, final FeatureContext featureContext) {
    if (DO_NOT_VALIDATE_FILTERED && vc.isFiltered()) {
        return;
    }
    // get the true reference allele
    final Allele reportedRefAllele = vc.getReference();
    final int refLength = reportedRefAllele.length();

    final Allele observedRefAllele = hasReference() ? Allele.create(Arrays.copyOf(ref.getBases(), refLength)) : null;

    final Set<String> rsIDs = getRSIDs(featureContext);

    if (VALIDATE_GVCF) {
        final SimpleInterval refInterval = ref.getInterval();

        validateVariantsOrder(vc);

        // GenomeLocSortedSet will automatically merge intervals that are overlapping when setting `mergeIfIntervalOverlaps`
        // to true.  In a GVCF most blocks are adjacent to each other so they wouldn't normally get merged.  We check
        // if the current record is adjacent to the previous record and "overlap" them if they are so our set is as
        // small as possible while still containing the same bases.
        final int start = (previousInterval != null && previousInterval.overlapsWithMargin(refInterval, 1)) ?
                previousInterval.getStart() : refInterval.getStart();
        final int end = (previousInterval != null && previousInterval.overlapsWithMargin(refInterval, 1)) ?
                Math.max(previousInterval.getEnd(), vc.getEnd()) : vc.getEnd();
        final GenomeLoc possiblyMergedGenomeLoc = genomeLocSortedSet.getGenomeLocParser().createGenomeLoc(refInterval.getContig(), start, end);
        genomeLocSortedSet.add(possiblyMergedGenomeLoc, true);

        previousInterval = new SimpleInterval(possiblyMergedGenomeLoc);
        previousStart = vc.getStart();
        validateGVCFVariant(vc);
    }

    for (final ValidationType t : validationTypes) {
        try{
            applyValidationType(vc, reportedRefAllele, observedRefAllele, rsIDs, t);
        } catch (TribbleException e) {
            throwOrWarn(new UserException.FailsStrictValidation(drivingVariantFile, t, e.getMessage()));
        }
    }
}
 

private byte[] handlePosition(SimpleInterval interval, byte base, FeatureContext features) {
    if (deletionBasesRemaining > 0) {
        deletionBasesRemaining--;
        return NO_BASES;
    }

    // If we have a mask at this site, use it
    if ( snpmaskPriority ){
        if (isMasked(features) )
            return N_BYTES;
    }

    // Check to see if we have a called snp
    for ( final VariantContext vc : features.getValues(variants) ) {
        if ( vc.isFiltered() || vc.getStart() != interval.getStart()  )
            continue;

        if ( vc.isSimpleDeletion()) {
            deletionBasesRemaining = vc.getReference().length() - 1;
            // delete the next n bases, not this one
            return baseToByteArray(base);
        } else if ( vc.isSimpleInsertion() || vc.isSNP() ) {
            // Get the first alt allele that is not a spanning deletion. If none present, use the empty allele
            final Optional<Allele> optionalAllele = getFirstConcreteAltAllele(vc.getAlternateAlleles());
            final Allele allele = optionalAllele.orElseGet(() -> Allele.create(EMPTY_BASE, false));
            if ( vc.isSimpleInsertion() ) {
                return allele.getBases();
            } else {
                final String iupacBase = (iupacSample != null) ? getIUPACBase(vc.getGenotype(iupacSample)) : allele.toString();
                return iupacBase.getBytes();
            }
        }
    }

    if ( !snpmaskPriority ){
        if ( isMasked(features)) {
            return N_BYTES;
        }
    }

    // if we got here then we're just ref
    return baseToByteArray(base);
}
 

private TruthVersusEval nextEvalOnlyVariant() {
    final VariantContext evalVariant = evalIterator.next();
    return evalVariant.isFiltered() ? TruthVersusEval.filteredTrueNegative(evalVariant) : TruthVersusEval.falsePositive(evalVariant);
}