Java Code Examples for htsjdk.variant.variantcontext.Genotype#hasAD()

@NotNull
public static VariantContext enrich(@NotNull final VariantContext context) {
    if (!context.isIndel() && !context.isSNP() ) {
        return context;
    }

    final VariantContextBuilder contextBuilder = new VariantContextBuilder(context).noGenotypes();

    final List<Allele> alleles = context.getAlleles();
    final Function<Allele, String> alleleKey = alleleKey(context);

    final List<Genotype> updatedGenotypes = Lists.newArrayList();
    for (Genotype genotype : context.getGenotypes()) {
        if (!genotype.hasAD() && hasRequiredAttributes(genotype, alleles, alleleKey)) {
            final GenotypeBuilder genotypeBuilder = new GenotypeBuilder(genotype).AD(readAD(genotype, alleles, alleleKey));
            updatedGenotypes.add(genotypeBuilder.make());
        } else {
            updatedGenotypes.add(genotype);
        }
    }

    return contextBuilder.genotypes(updatedGenotypes).make();
}
 

/**
 * @param genotype The genotype to test whether we can attempt to validate.  Never {@code null}
 * @param referenceAllele The reference allele (from parent VariantContext).  Never {@code null}
 * @return whether this class can even attempt a validation of the genotype in question.
 */
public static boolean isAbleToValidateGenotype(final Genotype genotype, final Allele referenceAllele) {
    Utils.nonNull(genotype);
    Utils.nonNull(referenceAllele);

    // In order to proceed, we have some assumptions:
    //  - The genotype has a ploidy of 2
    //  - The genotype is a simple indel or a xNP
    //  - The first allele of the genotype is reference
    final boolean isDiploid = genotype.getAlleles().size() == 2;
    final boolean doesGenotypeHaveReference = genotype.getAllele(0).equals(referenceAllele);
    final boolean isReferenceNotSymbolic = !referenceAllele.isSymbolic();
    final VariantContext.Type variantType = GATKVariantContextUtils.typeOfVariant(genotype.getAllele(0), genotype.getAllele(1));
    final boolean isValidatableVariantType = VALIDATABLE_TYPES.contains(variantType)
                    && !GATKVariantContextUtils.isComplexIndel(genotype.getAllele(0), genotype.getAllele(1));
    final boolean hasKnownCoverage = genotype.hasAD() && (genotype.getAD().length == 2);
    final boolean isValidateable = (isDiploid && doesGenotypeHaveReference && isValidatableVariantType && hasKnownCoverage && isReferenceNotSymbolic);
    if (!isValidateable) {
        logger.info("Cannot validate genotype: " + genotype + "  ploidy2: " + isDiploid +
                "  genotypeHasReferenceAllele: " + doesGenotypeHaveReference + "   validatableVariant: " + isValidatableVariantType +
                "  hasCompleteAD field: " + hasKnownCoverage + "  isNonSymbolicReference: " + isReferenceNotSymbolic);
    }

    return isValidateable;
}
 

private Map<Allele, Integer> getADcounts(final VariantContext vc) {
    final GenotypesContext genotypes = vc.getGenotypes();
    if ( genotypes == null || genotypes.size() == 0 ) {
        genotype_logger.warn("VC does not have genotypes -- annotations were calculated in wrong order");
        return null;
    }

    final Map<Allele, Integer> variantADs = new HashMap<>();
    for(final Allele a : vc.getAlleles())
        variantADs.put(a,0);

    for (final Genotype gt : vc.getGenotypes()) {
        if(gt.hasAD()) {
            final int[] ADs = gt.getAD();
            for (int i = 1; i < vc.getNAlleles(); i++) {
                variantADs.put(vc.getAlternateAllele(i - 1), variantADs.get(vc.getAlternateAllele(i - 1)) + ADs[i]); //here -1 is to reconcile allele index with alt allele index
            }
        }
    }
    return variantADs;
}
 

private void run() {
    fileWriter.writeHeader(fileReader.getFileHeader());
    for (final VariantContext context : fileReader) {

        final VariantContextBuilder builder = new VariantContextBuilder(context);
        final List<Genotype> genotypes = Lists.newArrayList();

        for (int genotypeIndex = 0; genotypeIndex < context.getGenotypes().size(); genotypeIndex++) {
            Genotype genotype = context.getGenotype(genotypeIndex);

            if (genotype.hasAD() && genotype.hasDP()) {
                int[] ad = genotype.getAD();
                int total = 0;
                boolean updateRecord = false;
                for (int i = 0; i < ad.length; i++) {
                    int adValue = ad[i];
                    if (adValue < 0) {
                        updateRecord = true;
                        ad[i] = Math.abs(adValue);
                    }
                    total += Math.abs(adValue);
                }

                if (updateRecord) {
                    LOGGER.info("Updated entry at {}:{}", context.getContig(), context.getStart());
                    Genotype updated = new GenotypeBuilder(genotype).AD(ad).DP(total).make();
                    genotypes.add(updated);
                }
                else {
                    genotypes.add(genotype);
                }
            }
        }

        builder.genotypes(genotypes);
        fileWriter.add(builder.make());
    }
}
 

@NotNull
public VariantContext enrich(@NotNull final VariantContext variant) {
    final Genotype genotype = variant.getGenotype(tumorSample);
    if (genotype != null && genotype.hasAD() && HumanChromosome.contains(variant.getContig())) {
        final GenomePosition position = GenomePositions.create(variant.getContig(), variant.getStart());
        final AllelicDepth depth = AllelicDepth.fromGenotype(genotype);
        enrich(position, depth, PurityAdjustedSomaticVariantBuilder.fromVariantContex(variant));
    }
    return variant;
}
 

public static int getDepth(final GenotypesContext genotypes, final AlleleLikelihoods<GATKRead, Allele> likelihoods) {
    int depth = 0;
    int ADrestrictedDepth = 0;

    for ( final Genotype genotype : genotypes ) {
        // we care only about variant calls with likelihoods
        if ( !genotype.isHet() && !genotype.isHomVar() ) {
            continue;
        }

        // if we have the AD values for this sample, let's make sure that the variant depth is greater than 1!
        if ( genotype.hasAD() ) {
            final int[] AD = genotype.getAD();
            final int totalADdepth = (int) MathUtils.sum(AD);
            if ( totalADdepth != 0 ) {
                if (totalADdepth - AD[0] > 1) {
                    ADrestrictedDepth += totalADdepth;
                }
                depth += totalADdepth;
                continue;
            }
        }
        // if there is no AD value or it is a dummy value, we want to look to other means to get the depth
        if (likelihoods != null) {
            depth += likelihoods.sampleEvidenceCount(likelihoods.indexOfSample(genotype.getSampleName()));
        } else if ( genotype.hasDP() ) {
            depth += genotype.getDP();
        }
    }

    // if the AD-restricted depth is a usable value (i.e. not zero), then we should use that one going forward
    if ( ADrestrictedDepth > 0 ) {
        depth = ADrestrictedDepth;
    }
    return depth;
}
 

static boolean containsAllelicDepth(final Genotype genotype) {
    return genotype.hasAD() && genotype.getAD().length > 1;
}
 

/**
 * method to swap the reference and alt alleles of a bi-allelic, SNP
 *
 * @param vc                   the {@link VariantContext} (bi-allelic SNP) that needs to have it's REF and ALT alleles swapped.
 * @param annotationsToReverse INFO field annotations (of double value) that will be reversed (x->1-x)
 * @param annotationsToDrop    INFO field annotations that will be dropped from the result since they are invalid when REF and ALT are swapped
 * @return a new {@link VariantContext} with alleles swapped, INFO fields modified and in the genotypes, GT, AD and PL corrected appropriately
 */
public static VariantContext swapRefAlt(final VariantContext vc, final Collection<String> annotationsToReverse, final Collection<String> annotationsToDrop) {

    if (!vc.isBiallelic() || !vc.isSNP()) {
        throw new IllegalArgumentException("swapRefAlt can only process biallelic, SNPS, found " + vc.toString());
    }

    final VariantContextBuilder swappedBuilder = new VariantContextBuilder(vc);

    swappedBuilder.attribute(SWAPPED_ALLELES, true);

    // Use getBaseString() (rather than the Allele itself) in order to create new Alleles with swapped
    // reference and non-variant attributes
    swappedBuilder.alleles(Arrays.asList(vc.getAlleles().get(1).getBaseString(), vc.getAlleles().get(0).getBaseString()));

    final Map<Allele, Allele> alleleMap = new HashMap<>();

    // A mapping from the old allele to the new allele, to be used when fixing the genotypes
    alleleMap.put(vc.getAlleles().get(0), swappedBuilder.getAlleles().get(1));
    alleleMap.put(vc.getAlleles().get(1), swappedBuilder.getAlleles().get(0));

    final GenotypesContext swappedGenotypes = GenotypesContext.create(vc.getGenotypes().size());
    for (final Genotype genotype : vc.getGenotypes()) {
        final List<Allele> swappedAlleles = new ArrayList<>();
        for (final Allele allele : genotype.getAlleles()) {
            if (allele.isNoCall()) {
                swappedAlleles.add(allele);
            } else {
                swappedAlleles.add(alleleMap.get(allele));
            }
        }
        // Flip AD
        final GenotypeBuilder builder = new GenotypeBuilder(genotype).alleles(swappedAlleles);
        if (genotype.hasAD() && genotype.getAD().length == 2) {
            final int[] ad = ArrayUtils.clone(genotype.getAD());
            ArrayUtils.reverse(ad);
            builder.AD(ad);
        } else {
            builder.noAD();
        }

        //Flip PL
        if (genotype.hasPL() && genotype.getPL().length == 3) {
            final int[] pl = ArrayUtils.clone(genotype.getPL());
            ArrayUtils.reverse(pl);
            builder.PL(pl);
        } else {
            builder.noPL();
        }
        swappedGenotypes.add(builder.make());
    }
    swappedBuilder.genotypes(swappedGenotypes);

    for (final String key : vc.getAttributes().keySet()) {
        if (annotationsToDrop.contains(key)) {
            swappedBuilder.rmAttribute(key);
        } else if (annotationsToReverse.contains(key) && !vc.getAttributeAsString(key, "").equals(VCFConstants.MISSING_VALUE_v4)) {
            final double attributeToReverse = vc.getAttributeAsDouble(key, -1);

            if (attributeToReverse < 0 || attributeToReverse > 1) {
                log.warn("Trying to reverse attribute " + key +
                        " but found value that isn't between 0 and 1: (" + attributeToReverse + ") in variant " + vc + ". Results might be wrong.");
            }
            swappedBuilder.attribute(key, 1 - attributeToReverse);
        }
    }

    return swappedBuilder.make();
}
 

private static final int altCount(final Genotype g) {
    return g.hasAD() ? (int) MathUtils.sum(g.getAD()) - g.getAD()[0] : 0;
}
 

/**
 * Create the MAF count fields (See {@link CustomMafFuncotationCreator#COUNT_FIELD_NAMES}) for a single pair in the given variant.
 *
 * If no samples can be linked to a tumor and normal combination, then blank values are generated for the count fields.
 * @param variant Never {@code null}
 * @param tnPair Never {@code null}
 * @return List of funcotations for the pair.  This list will usually be of length one, unless multiallelics are involved.
 *  Never {@code null}  Field values can be blank:
 *  - if there are no AD or AF format fields
 *  - if the specified pair is not in the variant sample list.
 */
private static List<Funcotation> createCustomMafCountFields(final VariantContext variant, final TumorNormalPair tnPair) {
    Utils.nonNull(variant);
    Utils.nonNull(tnPair);
    final List<Funcotation> result = new ArrayList<>();

    for (int i = 0; i < variant.getAlternateAlleles().size(); i++) {
        final Allele allele = variant.getAlternateAllele(i);
        final String tumorSampleName = tnPair.getTumor();
        final String normalSampleName = tnPair.getNormal();

        // This code assumes that the AD field is of Type "R"
        final Genotype tumorGenotype = variant.getGenotype(tumorSampleName);
        final boolean hasTumorAD = (tumorGenotype != null) && (tumorGenotype.hasAD());
        final boolean hasTumorAF = (tumorGenotype != null) && (tumorGenotype.hasAnyAttribute(GATKVCFConstants.ALLELE_FRACTION_KEY));

        // Just make a quick check that the tumor genotype has allelic depth for both a ref and alt
        // This is REQUIRED, but can be missing and our validations seem to miss it:
        if ( hasTumorAD && tumorGenotype.getAD().length < 2 ) {
            throw new UserException("Allelic Depth (AD field) for Variant[" +
                    variant.getContig() + ":" + variant.getStart() + "_" +
                    variant.getReference().toString() + "->" + variant.getAlternateAlleles().get(0).toString() + ']' +
                    " Does not contain both a REF and an ALT value (only one value is present)!");
        }

        final String tAltCount = hasTumorAD ? Integer.toString(tumorGenotype.getAD()[i + 1]) : "";
        final String tRefCount = hasTumorAD ? Integer.toString(tumorGenotype.getAD()[0]) : "";
        final String tumorFAllAlleles = hasTumorAF ? tumorGenotype.getAnyAttribute(GATKVCFConstants.ALLELE_FRACTION_KEY).toString() : "";
        final String tumorF = hasTumorAF ? StringUtils.split(tumorFAllAlleles, ",")[i] : "";

        final Genotype normalGenotype = variant.getGenotype(normalSampleName);
        final boolean hasNormalAD = (normalGenotype != null) && (normalGenotype.hasAD());
        final String nAltCount = hasNormalAD ? Integer.toString(normalGenotype.getAD()[i + 1]) : "";
        final String nRefCount = hasNormalAD ? Integer.toString(normalGenotype.getAD()[0]) : "";

        final List<String> fieldValues = Arrays.asList(
                tAltCount,
                tRefCount,
                nAltCount,
                nRefCount,
                tumorF);

        result.add(TableFuncotation.create(COUNT_FIELD_NAMES, fieldValues, allele, MafOutputRendererConstants.MAF_COUNT_RENDERING_DATASOURCE_DUMMY_NAME, createCustomMafCountFieldsMetadata()));
    }

    return result;
}