Java Allele.create方法代码示例

import htsjdk.variant.variantcontext.Allele; //导入方法依赖的package包/类
/**
 * Trim up alleles in inputVC, cutting out all bases up to fwdTrimEnd inclusive and
 * the last revTrim bases from the end
 *
 * @param inputVC a non-null input VC
 * @param fwdTrimEnd bases up to this index (can be -1) will be removed from the start of all alleles
 * @param revTrim the last revTrim bases of each allele will be clipped off as well
 * @return a non-null VariantContext (may be == to inputVC) with trimmed up alleles
 */
protected static VariantContext trimAlleles(final VariantContext inputVC,
                                            final int fwdTrimEnd,
                                            final int revTrim) {
    if( fwdTrimEnd == -1 && revTrim == 0 ) // nothing to do, so just return inputVC unmodified
        return inputVC;

    final List<Allele> alleles = new LinkedList<>();
    final Map<Allele, Allele> originalToTrimmedAlleleMap = new HashMap<>();

    for (final Allele a : inputVC.getAlleles()) {
        if (a.isSymbolic()) {
            alleles.add(a);
            originalToTrimmedAlleleMap.put(a, a);
        } else {
            // get bases for current allele and create a new one with trimmed bases
            final byte[] newBases = Arrays.copyOfRange(a.getBases(), fwdTrimEnd+1, a.length()-revTrim);
            final Allele trimmedAllele = Allele.create(newBases, a.isReference());
            alleles.add(trimmedAllele);
            originalToTrimmedAlleleMap.put(a, trimmedAllele);
        }
    }

    // now we can recreate new genotypes with trimmed alleles
    final AlleleMapper alleleMapper = new AlleleMapper(originalToTrimmedAlleleMap);
    final GenotypesContext genotypes = updateGenotypesWithMappedAlleles(inputVC.getGenotypes(), alleleMapper);

    final int start = inputVC.getStart() + (fwdTrimEnd + 1);
    final VariantContextBuilder builder = new VariantContextBuilder(inputVC);
    builder.start(start);
    builder.stop(start + alleles.get(0).length() - 1);
    builder.alleles(alleles);
    builder.genotypes(genotypes);
    return builder.make();
}
 

import htsjdk.variant.variantcontext.Allele; //导入方法依赖的package包/类
/**
 * Create a block substitution out of two variant contexts that start at the same position
 * <p>
 * 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) {
    if (vc1.getStart() != vc2.getStart())
        throw new IllegalArgumentException("vc1 and 2 must have the same start but got " + vc1 + " and " + vc2);
    if (!vc1.isBiallelic()) throw new IllegalArgumentException("vc1 must be biallelic");
    if (!vc1.isSNP()) {
        if (!((vc1.isSimpleDeletion() && vc2.isSimpleInsertion()) || (vc1.isSimpleInsertion() && vc2.isSimpleDeletion())))
            throw new IllegalArgumentException("Can only merge single insertion with deletion (or vice versa) but got " + vc1 + " merging with " + vc2);
    } else if (vc2.isSNP()) {
        throw new IllegalArgumentException("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();
}
 

import htsjdk.variant.variantcontext.Allele; //导入方法依赖的package包/类
/**
 * This method does a couple of things:
 * <ul><li>
 *     remaps the vc alleles considering the differences between the final reference allele and its own reference,</li>
 * <li>
 *     collects alternative alleles present in variant context and add them to the {@code finalAlleles} set.
 * </li></ul>
 *
 * @param vcAlleles the variant context allele list.
 * @param refAllele final reference allele.
 * @param finalAlleles where to add the final set of non-ref called alleles.
 * @return never {@code null}
 */
//TODO as part of a larger refactoring effort {@link #remapAlleles} can be merged with {@link GATKVariantContextUtils#remapAlleles}.
private static List<Allele> remapAlleles(final List<Allele> vcAlleles, final Allele refAllele, final LinkedHashSet<Allele> finalAlleles) {
    final Allele vcRef = vcAlleles.get(0);
    if (!vcRef.isReference()) throw new IllegalStateException("the first allele of the vc allele list must be reference");
    final byte[] refBases = refAllele.getBases();
    final int extraBaseCount = refBases.length - vcRef.getBases().length;
    if (extraBaseCount < 0) throw new IllegalStateException("the wrong reference was selected");
    final List<Allele> result = new ArrayList<>(vcAlleles.size());

    for (final Allele a : vcAlleles) {
        if (a.isReference()) {
            result.add(refAllele);
        } else if (a.isSymbolic()) {
            result.add(a);
            // we always skip <NON_REF> when adding to finalAlleles this is done outside if applies.
            if (!a.equals(GATKVariantContextUtils.NON_REF_SYMBOLIC_ALLELE))
                finalAlleles.add(a);
        } else if (a.isCalled()) {
            final Allele newAllele;
            if (extraBaseCount > 0) {
                final byte[] oldBases = a.getBases();
                final byte[] newBases = Arrays.copyOf(oldBases,oldBases.length + extraBaseCount);
                System.arraycopy(refBases,refBases.length - extraBaseCount,newBases,oldBases.length,extraBaseCount);
                newAllele = Allele.create(newBases,false);
            } else
                newAllele = a;
            result.add(newAllele);
            finalAlleles.add(newAllele);
        } else { // NO_CALL and strange miscellanea
            result.add(a);
        }
    }
    return result;
}
 

import htsjdk.variant.variantcontext.Allele; //导入方法依赖的package包/类
protected VariantContext createMergedVariantContext( final VariantContext thisVC, final VariantContext nextVC, final byte[] ref, final GenomeLoc refLoc ) {
    final int thisStart = thisVC.getStart();
    final int nextStart = nextVC.getStart();
    byte[] refBases = new byte[]{};
    byte[] altBases = new byte[]{};
    refBases = ArrayUtils.addAll(refBases, thisVC.getReference().getBases());
    altBases = ArrayUtils.addAll(altBases, thisVC.getAlternateAllele(0).getBases());
    int locus;
    for( locus = thisStart + refBases.length; locus < nextStart; locus++ ) {
        final byte refByte = ref[locus - refLoc.getStart()];
        refBases = ArrayUtils.add(refBases, refByte);
        altBases = ArrayUtils.add(altBases, refByte);
    }
    refBases = ArrayUtils.addAll(refBases, ArrayUtils.subarray(nextVC.getReference().getBases(), locus > nextStart ? 1 : 0, nextVC.getReference().getBases().length)); // special case of deletion including the padding base of consecutive indel
    altBases = ArrayUtils.addAll(altBases, nextVC.getAlternateAllele(0).getBases());

    int iii = 0;
    if( refBases.length == altBases.length ) { // insertion + deletion of same length creates an MNP --> trim common prefix bases off the beginning of the allele
        while( iii < refBases.length && refBases[iii] == altBases[iii] ) { iii++; }
        if ( iii == refBases.length ) {
            // we've become a null allele, such as with CA/C + A/AA -> CA/CA => after trimming there's nothing left
            // so return a null variant context so we can eliminate the variants from consideration
            return null;
        }
    }


    final Allele refAllele = Allele.create( ArrayUtils.subarray(refBases, iii, refBases.length), true );
    final Allele altAllele =  Allele.create( ArrayUtils.subarray(altBases, iii, altBases.length), false );
    return new VariantContextBuilder("merged", thisVC.getChr(), thisVC.getStart() + iii, nextVC.getEnd(), Arrays.asList(refAllele, altAllele)).make();
}
 

import htsjdk.variant.variantcontext.Allele; //导入方法依赖的package包/类
/**
 * Splits the alleles for the provided variant context into its primitive parts.
 * Requires that the input VC be bi-allelic, so calling methods should first call splitVariantContextToBiallelics() if needed.
 * Currently works only for MNPs.
 *
 * @param vc  the non-null VC to split
 * @return a non-empty list of VCs split into primitive parts or the original VC otherwise
 */
public static List<VariantContext> splitIntoPrimitiveAlleles(final VariantContext vc) {
    if ( vc == null )
        throw new IllegalArgumentException("Trying to break a null Variant Context into primitive parts");

    if ( !vc.isBiallelic() )
        throw new IllegalArgumentException("Trying to break a multi-allelic Variant Context into primitive parts");

    // currently only works for MNPs
    if ( !vc.isMNP() )
        return Arrays.asList(vc);

    final byte[] ref = vc.getReference().getBases();
    final byte[] alt = vc.getAlternateAllele(0).getBases();

    if ( ref.length != alt.length )
        throw new IllegalStateException("ref and alt alleles for MNP have different lengths");

    final List<VariantContext> result = new ArrayList<>(ref.length);

    for ( int i = 0; i < ref.length; i++ ) {

        // if the ref and alt bases are different at a given position, create a new SNP record (otherwise do nothing)
        if ( ref[i] != alt[i] ) {

            // create the ref and alt SNP alleles
            final Allele newRefAllele = Allele.create(ref[i], true);
            final Allele newAltAllele = Allele.create(alt[i], false);

            // create a new VariantContext with the new SNP alleles
            final VariantContextBuilder newVC = new VariantContextBuilder(vc).start(vc.getStart() + i).stop(vc.getStart() + i).alleles(Arrays.asList(newRefAllele, newAltAllele));

            // create new genotypes with updated alleles
            final Map<Allele, Allele> alleleMap = new HashMap<>();
            alleleMap.put(vc.getReference(), newRefAllele);
            alleleMap.put(vc.getAlternateAllele(0), newAltAllele);
            final GenotypesContext newGenotypes = updateGenotypesWithMappedAlleles(vc.getGenotypes(), new AlleleMapper(alleleMap));

            result.add(newVC.genotypes(newGenotypes).make());
        }
    }

    if ( result.isEmpty() )
        result.add(vc);

    return result;
}
 

import htsjdk.variant.variantcontext.Allele; //导入方法依赖的package包/类
/**
 * Calculate the genotype likelihoods for the sample in pileup for being hom-ref contrasted with being ref vs. alt
 *
 * @param sampleName target sample name.
 * @param ploidy target sample ploidy.
 * @param genotypingModel model to calculate likelihoods and genotypes.
 * @param pileup the read backed pileup containing the data we want to evaluate
 * @param refBase the reference base at this pileup position
 * @param minBaseQual the min base quality for a read in the pileup at the pileup position to be included in the calculation
 * @param hqSoftClips running average data structure (can be null) to collect information about the number of high quality soft clips
 * @return a RefVsAnyResult genotype call.
 */
public RefVsAnyResult calcGenotypeLikelihoodsOfRefVsAny(final String sampleName, final int ploidy,
                                                        final GenotypingModel genotypingModel,
                                                        final ReadBackedPileup pileup, final byte refBase, final byte minBaseQual, final MathUtils.RunningAverage hqSoftClips) {
    final AlleleList<Allele> alleleList = new IndexedAlleleList<>(Allele.create(refBase,true),GATKVariantContextUtils.NON_REF_SYMBOLIC_ALLELE);
    // Notice that the sample name is rather irrelevant as this information is never used, just need to be the same in both lines bellow.

    final int maximumReadCount = pileup.getReads().size();

    final List<GATKSAMRecord> reads = new ArrayList<>(maximumReadCount);
    final double[][] likelihoods = new double[2][maximumReadCount];
    final int[] adCounts = new int[2];
    int nextIndex = 0;
    for (final PileupElement p : pileup) {
        final byte qual = p.isDeletion() ? REF_MODEL_DELETION_QUAL : p.getQual();
        if (!p.isDeletion() && qual <= minBaseQual)
            continue;
        final GATKSAMRecord read = p.getRead();
        reads.add(read);
        final boolean isAlt = p.getBase() != refBase || p.isDeletion() || p.isBeforeDeletionStart()
                || p.isAfterDeletionEnd() || p.isBeforeInsertion() || p.isAfterInsertion() || p.isNextToSoftClip();
        final int bestAllele;
        final int worstAllele;
        if (isAlt) {
            bestAllele = 1;
            worstAllele = 0;
        } else {
            bestAllele = 0;
            worstAllele = 1;
        }

        likelihoods[bestAllele][nextIndex] = QualityUtils.qualToProbLog10(qual);
        likelihoods[worstAllele][nextIndex++] = QualityUtils.qualToErrorProbLog10(qual) + MathUtils.LOG_ONE_THIRD;
        adCounts[bestAllele]++;
        if (isAlt && hqSoftClips != null && p.isNextToSoftClip())
            hqSoftClips.add(AlignmentUtils.calcNumHighQualitySoftClips(read, (byte) 28));
    }

    final Map<String,List<GATKSAMRecord>> sampleToReads = Collections.singletonMap(sampleName,reads);
    final ReadLikelihoods<Allele> readLikelihoods = new ReadLikelihoods<>(new IndexedSampleList(sampleName),alleleList,sampleToReads, true);
    final ReadLikelihoods.Matrix<Allele> sampleLikelihoods = readLikelihoods.sampleMatrix(0);
    final int readCount = sampleLikelihoods.readCount();
    for (int i = 0; i < readCount; i++) {
        sampleLikelihoods.set(0,i,likelihoods[0][i]);
        sampleLikelihoods.set(1,i,likelihoods[1][i]);
    }

    final PloidyModel ploidyModel = new HomogeneousPloidyModel(new IndexedSampleList(sampleName),ploidy);
    final GenotypingLikelihoods<Allele> genotypingLikelihoods = genotypingModel.calculateLikelihoods(alleleList, new GenotypingData<>(ploidyModel, readLikelihoods));
    final double[] genotypeLikelihoodArray = genotypingLikelihoods.sampleLikelihoods(0).getAsVector();
    final RefVsAnyResult result = new RefVsAnyResult(genotypeLikelihoodArray.length);
    System.arraycopy(genotypeLikelihoodArray,0,result.genotypeLikelihoods,0,genotypeLikelihoodArray.length);
    System.arraycopy(adCounts,0,result.AD_Ref_Any,0,2);
    return result;
}
 

import htsjdk.variant.variantcontext.Allele; //导入方法依赖的package包/类
private List<Allele> consensusCountsToAlleles(final ReferenceContext ref,
                                              final Map<String, Integer> consensusIndelStrings) {
    final GenomeLoc loc = ref.getLocus();
    final Collection<VariantContext> vcs = new ArrayList<VariantContext>();
    int maxAlleleCnt = 0;
    Allele refAllele, altAllele;

    for (final Map.Entry<String, Integer> elt : consensusIndelStrings.entrySet()) {
        final String s = elt.getKey();
        final int curCnt = elt.getValue();
        int stop = 0;

        // if observed count if above minimum threshold, we will genotype this allele
        if (curCnt < minIndelCountForGenotyping)
            continue;

        if (s.startsWith("D")) {
            // get deletion length
            final int dLen = Integer.valueOf(s.substring(1));
            // get ref bases of accurate deletion
            final int startIdxInReference = 1 + loc.getStart() - ref.getWindow().getStart();
            stop = loc.getStart() + dLen;
            final byte[] refBases = Arrays.copyOfRange(ref.getBases(), startIdxInReference - 1, startIdxInReference + dLen);   // add reference padding

            if (Allele.acceptableAlleleBases(refBases, false)) {
                refAllele = Allele.create(refBases, true);
                altAllele = Allele.create(ref.getBase(), false);
            }
            else continue; // don't go on with this allele if refBases are non-standard
        } else {
            // insertion case
            final String insertionBases = (char)ref.getBase() + s;  // add reference padding
            if (Allele.acceptableAlleleBases(insertionBases, false)) { // don't allow N's in insertions
                refAllele = Allele.create(ref.getBase(), true);
                altAllele = Allele.create(insertionBases, false);
                stop = loc.getStart();
            }
            else continue; // go on to next allele if consensus insertion has any non-standard base.
        }


        final VariantContextBuilder builder = new VariantContextBuilder().source("");
        builder.loc(loc.getContig(), loc.getStart(), stop);
        builder.alleles(Arrays.asList(refAllele, altAllele));
        builder.noGenotypes();
        if (doMultiAllelicCalls) {
            vcs.add(builder.make());
            if (vcs.size() >= GenotypeLikelihoods.MAX_ALT_ALLELES_THAT_CAN_BE_GENOTYPED)
                break;
        } else if (curCnt > maxAlleleCnt) {
            maxAlleleCnt = curCnt;
            vcs.clear();
            vcs.add(builder.make());
        }
    }

    if (vcs.isEmpty())
        return Collections.emptyList(); // nothing else to do, no alleles passed minimum count criterion

    final VariantContext mergedVC = GATKVariantContextUtils.simpleMerge(vcs, null, GATKVariantContextUtils.FilteredRecordMergeType.KEEP_IF_ANY_UNFILTERED, GATKVariantContextUtils.GenotypeMergeType.UNSORTED, false, false, null, false, false);
    return mergedVC.getAlleles();
}
 

import htsjdk.variant.variantcontext.Allele; //导入方法依赖的package包/类
/**
 * Determines the ref allele given the provided reference base at this position
 *
 * @param VCs     collection of unsorted genomic VCs
 * @param loc     the current location
 * @param refBase the reference allele to use if all contexts in the VC are spanning
 * @return new Allele or null if no reference allele/base is available
 */
private static Allele determineReferenceAlleleGivenReferenceBase(final List<VariantContext> VCs, final GenomeLoc loc, final Byte refBase) {
    final Allele refAllele = GATKVariantContextUtils.determineReferenceAllele(VCs, loc);
    if ( refAllele == null )
        return ( refBase == null ? null : Allele.create(refBase, true) );
    return refAllele;
}