Java Code Examples for htsjdk.samtools.util.Locatable#getEnd()

/**
 * Given an interval and collections of targets and SNPs, returns a trimmed interval produced by removing the empty
 * portions at the start and the end of the original interval that do not overlap the targets and SNPs that overlap
 * with the original interval.  If this procedure would remove the entire interval, the original interval is
 * returned instead.  Note that this method will not expand an interval to the start of the first overlapping target
 * and the end of the last overlapping target; it will only shrink the interval or leave it alone.  This is to
 * avoid overlapping segments (which would occur if a SNP breakpoint fell within a target and the interval
 * were expanded, for example).
 */
public static SimpleInterval trimInterval(final Locatable interval,
                                          final TargetCollection<? extends Locatable> targets,
                                          final TargetCollection<? extends Locatable> snps) {
    Utils.nonNull(interval, "The interval cannot be null.");
    Utils.nonNull(targets, "The collection of targets cannot be null.");
    Utils.nonNull(snps, "The collection of SNPs cannot be null.");

    final IndexRange targetRange = targets.indexRange(interval);
    final IndexRange snpRange = snps.indexRange(interval);

    final int numTargetsInInterval = targetRange.size();
    final int numSNPsInInterval = snpRange.size();

    int start = interval.getStart();
    int end = interval.getEnd();

    if (numTargetsInInterval == 0 && numSNPsInInterval > 0) {
        //if there are no targets overlapping interval, use SNPs to determine trimmed interval
        start = snps.target(snpRange.from).getStart();
        end = snps.target(snpRange.to - 1).getEnd();
    } else if (numTargetsInInterval > 0) {
        //if interval start does not fall within first target, use start of first target as start of trimmed interval
        start = Math.max(start, targets.target(targetRange.from).getStart());
        //if interval end does not fall within last target, use end of last target as end of trimmed interval
        end = Math.min(end, targets.target(targetRange.to - 1).getEnd());

        if (numSNPsInInterval > 0) {
            //if there are also SNPs within interval, check to see if they give a larger trimmed interval
            start = Math.min(start, snps.target(snpRange.from).getStart());
            end = Math.max(end, snps.target(snpRange.to - 1).getEnd());
        }
    }
    if (start < interval.getStart() || end > interval.getEnd() || end < start) {
        throw new GATKException.ShouldNeverReachHereException("Something went wrong in trimming interval.");
    }
    return new SimpleInterval(interval.getContig(), start, end);
}
 

/**
 * @param transcript Never {@code null}
 * @param segment Genomic region to determine which exons are covered in the transcript.  Never {@code null}
 * @return Instance of {@link SegmentExonOverlaps}.  A typical value will be a number appended to a "-" or a "+".
 *   The exon numbers are 0-based.  "-" if it covers upstream exons (considering coding direction -- i.e. previous exons
 *   are higher genomic coordinates on negative strand transcripts) and "+" for downstream exons.
 *   For example:
 *   - "6+" would be sixth exon and above.  This would indicate downstream exons (higher genomic
 *     coordinates on positive coding directions and lower genomic coordinates when negative coding) are covered by
 *     the segment.
 *   - "" indicates that the endpoints of the segment do not cover any transcript exons.  So either the entire
 *      transcript is covered or none of the transcript is covered.
 *
 *   Note that this will return the first exon covered, as long as the breakpoint is in the transcript.  Even if the
 *    breakpoint itself does not cover an exon.
 *    For example, this will yield a non-blank value when a segment breakpoint is in an intron.
 */
public static SegmentExonOverlaps determineSegmentExonPosition(final GencodeGtfTranscriptFeature transcript, final Locatable segment) {

    Utils.nonNull(transcript);
    Utils.nonNull(segment);

    final String NOT_IN_TRANSCRIPT = "";

    // Internally, this method assumes that the exons are sorted in genomic order.  Which is NOT how these are
    //  stored in the GencodeGtfTranscriptFeature
    final List<GencodeGtfExonFeature> exons = transcript.getGenomicStrand() == Strand.NEGATIVE ?
            Lists.reverse(transcript.getExons()) : transcript.getExons();

    if (exons.size() == 0) {
        return new SegmentExonOverlaps(NOT_IN_TRANSCRIPT, NOT_IN_TRANSCRIPT);
    }

    final SimpleInterval segmentStart = new SimpleInterval(segment.getContig(), segment.getStart(), segment.getStart());
    final SimpleInterval segmentEnd = new SimpleInterval(segment.getContig(), segment.getEnd(), segment.getEnd());

    // Find the proper index for the start.
    // If the start of the segment does not overlap the first exon, but the segment does, then the start index is -1 (no overlap)
    final int inclusiveIndexPositiveDirectionStart = findInclusiveExonIndex(transcript.getGenomicStrand(), exons, segmentStart, true);

    // If the end of the segment does not overlap the last exon, but the segment does, then the end index is -1 (no overlap)
    final int inclusiveIndexPositiveDirectionEnd = findInclusiveExonIndex(transcript.getGenomicStrand(), exons, segmentEnd, false);

    // Construct the final strings
    final String startResult = inclusiveIndexPositiveDirectionStart != NO_EXON_OVERLAP ?
            inclusiveIndexPositiveDirectionStart + determineSegmentOverlapDirection(transcript.getGenomicStrand(), true)
            : NOT_IN_TRANSCRIPT;
    final String endResult = inclusiveIndexPositiveDirectionEnd != NO_EXON_OVERLAP ?
            inclusiveIndexPositiveDirectionEnd + determineSegmentOverlapDirection(transcript.getGenomicStrand(), false)
            : NOT_IN_TRANSCRIPT;

    return new SegmentExonOverlaps(startResult, endResult);
}
 

/**
 * Create a cDNA string for an intronic variant.
 * The cDNA string contains information about the coding position of a variant and the alleles involved.
 * If the transcript in which the variant occurs contains at least 1 exon, the string will be non-empty and of
 * the form:
 *     c.e[EXON NUMBER][+|-][BASES FROM EXON][REF ALLELE]>[ALT ALLELE]
 * Concretely:
 *     c.e2-1A>G
 * Where:
 *      2 = the number of the exon to which the given variant start is closest
 *     -1 = number of bases away from the exon (1 before)
 *      A = Reference allele
 *      G = Alternate allele
 * @param variantStart The start position (1-based, inclusive) in genomic coordinates of the variant.
 * @param exonList The {@link List<Locatable>} representing the exons in the transcript in which this variant occurs.  Must not be {@code null}.
 * @param strandCorrectedRefAllele A {@link String} containing the bases of the reference allele, which are correct for strandedness (i.e. if on a transcript on the {@link Strand#NEGATIVE} strand, the string has already been reverse-complemented).  Must not be {@code null}.
 * @param strandCorrectedAltAllele A {@link String} containing the bases of the alternate allele, which are correct for strandedness (i.e. if on a transcript on the {@link Strand#NEGATIVE} strand, the string has already been reverse-complemented).  Must not be {@code null}.
 * @return A {@link String} representing the cDNA change for the given data.  Will be empty if the given {@code exonList} is empty.
 */
public static String createIntronicCDnaString(final int variantStart,
                                              final List<? extends Locatable> exonList,
                                              final String strandCorrectedRefAllele,
                                              final String strandCorrectedAltAllele) {

    Utils.nonNull(exonList);
    Utils.nonNull(strandCorrectedRefAllele);
    Utils.nonNull(strandCorrectedAltAllele);

    // Get the exon that is closest to our variant:
    final int exonIndex = getClosestExonIndex(variantStart, exonList);

    if ( exonIndex != -1 ) {
        final Locatable closestExon = exonList.get(exonIndex);

        final int startDiff = variantStart - closestExon.getStart();
        final int endDiff =  variantStart - closestExon.getEnd();

        // Get the offset from our start:
        final int exonOffset;
        if ( Math.abs(startDiff) <= Math.abs(endDiff) ) {
            exonOffset = startDiff;
        }
        else {
            exonOffset = endDiff;
        }

        // Get the cDNA string itself:
        return "c.e" + (exonIndex+1) + (exonOffset < 0 ? "-" : "+") + Math.abs(exonOffset) + strandCorrectedRefAllele + ">" + strandCorrectedAltAllele;
    }
    else {
        return "NA";
    }
}
 

/**
 * Gets the distance, in base pairs between two locatables.  Looks for the closest endpoints, so the order of loc1
 *  and loc2 do not matter.
 *
 * Overlapping (or abutting) locatables will always return 0.
 *
 * Locatables on different contigs will have a distance of {@link Long#MAX_VALUE}
 * @param loc1 Never {@code null}
 * @param loc2 Never {@code null}
 * @return 0 or a positive number.
 */
private static long getDistance(final Locatable loc1, final Locatable loc2) {
    Utils.nonNull(loc1);
    Utils.nonNull(loc2);

    if (!loc1.getContig().equals(loc2.getContig())) {
        return Long.MAX_VALUE;
    }

    if (IntervalUtils.overlaps(loc1, loc2)) {
        return 0;
    }

    return (loc1.getEnd() < loc2.getStart()) ? (loc2.getStart() - loc1.getEnd()) : (loc1.getStart() - loc2.getEnd());
}
 

/**
 * Determines whether this interval contains the entire region represented by other
 * (in other words, whether it covers it).
 *
 * @param other interval to check
 * @return true if this interval contains all of the bases spanned by other, otherwise false
 */
public boolean contains( final Locatable other ) {
    if ( other == null || other.getContig() == null ) {
        return false;
    }

    return this.contig.equals(other.getContig()) && this.start <= other.getStart() && this.end >= other.getEnd();
}
 

/**
 * Determines whether this interval comes within "margin" of overlapping the provided locatable.
 * This is the same as plain overlaps if margin=0.
 *
 * @param other interval to check
 * @param margin how many bases may be between the two interval for us to still consider them overlapping; must be non-negative
 * @return true if this interval overlaps other, otherwise false
 * @throws IllegalArgumentException if margin is negative
 */
public boolean overlapsWithMargin(final Locatable other, final int margin) {
    if ( margin < 0 ) {
        throw new IllegalArgumentException("given margin is negative: " + margin +
                "\tfor this: " + toString() + "\tand that: " + (other == null ? "other is null" : other.toString()));
    }
    if ( other == null || other.getContig() == null ) {
        return false;
    }

    return this.contig.equals(other.getContig()) && this.start <= other.getEnd() + margin && other.getStart() - margin <= this.end;
}
 

protected static boolean[] calculateKnownSites( final GATKRead read, final Iterable<? extends Locatable> knownSites ) {
    final int readLength = read.getLength();
    final boolean[] knownSitesArray = new boolean[readLength];//initializes to all false
    final Cigar cigar = read.getCigar();
    final int softStart = read.getSoftStart();
    final int softEnd = read.getSoftEnd();
    for ( final Locatable knownSite : knownSites ) {
        if (knownSite.getEnd() < softStart || knownSite.getStart() > softEnd) {
            // knownSite is outside clipping window for the read, ignore
            continue;
        }
        final Pair<Integer, CigarOperator> featureStartAndOperatorOnRead = ReadUtils.getReadIndexForReferenceCoordinate(read, knownSite.getStart());
        int featureStartOnRead = featureStartAndOperatorOnRead.getLeft() == ReadUtils.READ_INDEX_NOT_FOUND ? 0 : featureStartAndOperatorOnRead.getLeft();
        if (featureStartAndOperatorOnRead.getRight() == CigarOperator.DELETION) {
            featureStartOnRead--;
        }

        final Pair<Integer, CigarOperator> featureEndAndOperatorOnRead = ReadUtils.getReadIndexForReferenceCoordinate(read, knownSite.getEnd());
        int featureEndOnRead = featureEndAndOperatorOnRead.getLeft() == ReadUtils.READ_INDEX_NOT_FOUND ? readLength : featureEndAndOperatorOnRead.getLeft();

        if( featureStartOnRead > readLength ) {
            featureStartOnRead = featureEndOnRead = readLength;
        }

        Arrays.fill(knownSitesArray, Math.max(0, featureStartOnRead), Math.min(readLength, featureEndOnRead + 1), true);
    }
    return knownSitesArray;
}
 

/**
 * Tests whether the first Locatable starts after the end of the second Locatable
 *
 * @param first first Locatable
 * @param second second Locatable
 * @param dictionary sequence dictionary used to determine contig ordering
 * @return true if first starts after the end of second, otherwise false
 */
public static boolean isAfter(final Locatable first, final Locatable second, final SAMSequenceDictionary dictionary) {
    Utils.nonNull(first);
    Utils.nonNull(second);
    Utils.nonNull(dictionary);

    final int contigComparison = compareContigs(first, second, dictionary);
    return contigComparison == 1 || (contigComparison == 0 && first.getStart() > second.getEnd());
}
 

/**
 * Tests whether the first Locatable ends before the start of the second Locatable
 *
 * @param first first Locatable
 * @param second second Locatable
 * @param dictionary sequence dictionary used to determine contig ordering
 * @return true if first ends before the start of second, otherwise false
 */
public static boolean isBefore(final Locatable first, final Locatable second, final SAMSequenceDictionary dictionary) {
    Utils.nonNull(first);
    Utils.nonNull(second);
    Utils.nonNull(dictionary);

    final int contigComparison = compareContigs(first, second, dictionary);
    return contigComparison == -1 || (contigComparison == 0 && first.getEnd() < second.getStart());
}
 

/**
 * Create a new Haplotype derived from this one that exactly spans the provided location
 *
 * Note that this haplotype must have a contain a genome loc for this operation to be successful.  If no
 * GenomeLoc is contained than @throws an IllegalStateException
 *
 * Also loc must be fully contained within this Haplotype's genomeLoc.  If not an IllegalArgumentException is
 * thrown.
 *
 * @param loc a location completely contained within this Haplotype's location
 * @return a new Haplotype within only the bases spanning the provided location, or null for some reason the haplotype would be malformed if
 */
public Haplotype trim(final Locatable loc) {
    Utils.nonNull( loc, "Loc cannot be null");
    Utils.nonNull(genomeLocation, "Cannot trim a Haplotype without containing GenomeLoc");
    Utils.validateArg(new SimpleInterval(genomeLocation).contains(loc), () -> "Can only trim a Haplotype to a containing span.  My loc is " + genomeLocation + " but wanted trim to " + loc);
    Utils.nonNull( getCigar(), "Cannot trim haplotype without a cigar " + this);

    final int newStart = loc.getStart() - this.genomeLocation.getStart();
    final int newStop = newStart + loc.getEnd() - loc.getStart();

    // note: the following returns null if the bases covering the ref interval start or end in a deletion.
    final byte[] newBases = AlignmentUtils.getBasesCoveringRefInterval(newStart, newStop, getBases(), 0, getCigar());

    if ( newBases == null || newBases.length == 0 ) { // we cannot meaningfully chop down the haplotype, so return null
        return null;
    }

    // note: trimCigarByReference does not remove leading or trailing indels, while getBasesCoveringRefInterval does remove bases
    // of leading and trailing insertions.  We must remove leading and trailing insertions from the Cigar manually.
    // we keep leading and trailing deletions because these are necessary for haplotypes to maintain consistent reference coordinates
    final Cigar newCigar = AlignmentUtils.trimCigarByReference(getCigar(), newStart, newStop).getCigar();
    final boolean leadingInsertion = !newCigar.getFirstCigarElement().getOperator().consumesReferenceBases();
    final boolean trailingInsertion = !newCigar.getLastCigarElement().getOperator().consumesReferenceBases();
    final int firstIndexToKeepInclusive = leadingInsertion ? 1 : 0;
    final int lastIndexToKeepExclusive = newCigar.numCigarElements() - (trailingInsertion ? 1 : 0);

    if (lastIndexToKeepExclusive <= firstIndexToKeepInclusive) {    // edge case of entire cigar is insertion
        return null;
    }

    final Cigar leadingIndelTrimmedNewCigar = !(leadingInsertion || trailingInsertion)  ? newCigar :
            new CigarBuilder(false).addAll(newCigar.getCigarElements().subList(firstIndexToKeepInclusive, lastIndexToKeepExclusive)).make();

    final Haplotype ret = new Haplotype(newBases, isReference());
    ret.setCigar(leadingIndelTrimmedNewCigar);
    ret.setGenomeLocation(loc);
    ret.setScore(score);
    ret.setAlignmentStartHapwrtRef(newStart + getAlignmentStartHapwrtRef());
    return ret;
}
 

/**
 * getVariantShardsFromInterval calculates *all* shards overlapping location.
 * @param location the range of sites determines which shards are overlapping
 * @return All overlapping VariantShards
 */
public static List<VariantShard> getVariantShardsFromInterval(final Locatable location) {
    if (location.getContig()==null) {
        // don't feed me unmapped reads!
        throw new GATKException("getVariantShardsFromInterval requires locations to be mapped");
    }
    List<VariantShard> shardList = new ArrayList<>();
    // Get all of the shard numbers that span the start and end of the interval.
    int startShard = location.getStart()/ VARIANT_SHARDSIZE;
    int endShard = location.getEnd()/ VARIANT_SHARDSIZE;
    for (int i = startShard; i <= endShard; ++i) {
        shardList.add(new VariantShard(i, location.getContig()));
    }
    return shardList;
}
 

public static double calculateDistance(final Locatable from, final Locatable to) {
    if (!from.getContig().equals(to.getContig())) {
        return Double.POSITIVE_INFINITY;
    } else {
        final double toMidpoint = (to.getStart() + to.getEnd())/2;
        final double fromMidpoint = (from.getStart() + from.getEnd())/2;
        return Math.abs(toMidpoint - fromMidpoint);
    }
}
 

private boolean contiguous(final Locatable that) {
    Utils.nonNull(that);
    return this.getContig().equals(that.getContig()) && this.getStart() <= that.getEnd() + 1 && that.getStart() <= this.getEnd() + 1;
}
 

/**
 * Get the position (1-based, inclusive) of the given {@link VariantContext} start relative to the transcript it appears in.
 * The transcript is specified by {@code sortedTranscriptExonList}.
 * @param variant The {@link VariantContext} of which to find the start position in the given transcript (must not be {@code null}).
 * @param exons {@link List} of {@link Locatable}s representing the exons in the transcript in which the given {@code variant} occurs.
 * @param strand The {@link Strand} on which the {@code variant} occurs.
 * @return The start position (1-based, inclusive) of the given {@code variant} in the transcript in which it appears.
 */
public static int getTranscriptAlleleStartPosition(final VariantContext variant,
                                                   final List<? extends Locatable> exons,
                                                   final Strand strand) {
    Utils.nonNull(variant);
    Utils.nonNull(exons);
    assertValidStrand(strand);

    // Set up our position variable:
    int position;

    // NOTE: We don't need to worry about UTRs in here - all UTRs occur somewhere in an exon in GENCODE.

    // Filter the elements by whether they come before the variant in the transcript and
    // then sort them by their order in the transcript:
    final List<Locatable> sortedFilteredExons;
    if ( strand == Strand.POSITIVE) {
        sortedFilteredExons = exons.stream()
                .filter(e -> e.getStart() <= variant.getStart())
                .sorted(Comparator.comparingInt(Locatable::getStart))
                .collect(Collectors.toList());

        // We are guaranteed that the variant occurs in the last element of sortedTranscriptElements because of the sorting.
        // Add 1 to position to account for inclusive values:
        position = variant.getStart() - sortedFilteredExons.get(sortedFilteredExons.size()-1).getStart() + 1;
    }
    else {
        sortedFilteredExons = exons.stream()
                .filter(e -> e.getEnd() >= variant.getStart())
                .sorted(Comparator.comparingInt(Locatable::getStart).reversed())
                .collect(Collectors.toList());

        // We are guaranteed that the variant occurs in the last element of sortedTranscriptElements because of the sorting.
        // Add 1 to position to account for inclusive values:
        position = sortedFilteredExons.get(sortedFilteredExons.size()-1).getEnd() - variant.getStart() + 1;
    }

    // Add up the lengths of all exons before the last one:
    final int numExonsBeforeLast = sortedFilteredExons.size() - 1;
    for ( int i = 0; i < numExonsBeforeLast; ++i ) {
        final Locatable exon = sortedFilteredExons.get(i);

        // Add 1 to position to account for inclusive values:
        position += exon.getEnd() - exon.getStart() + 1;
    }

    return position;
}
 

/**
 * Get the overlapping exon start/stop as a {@link SimpleInterval} for the given altAllele / reference.
 * @param refAllele {@link Allele} for the given {@code altAllele}.  Must not be {@code null}.
 * @param altAllele {@link Allele} to locate on an exon.  Must not be {@code null}.
 * @param contig Contig on which the altAllele occurs.  Must not be {@code null}.
 * @param variantStart Start position (1-based, inclusive) of the given {@code altAllele}.  Must be > 0.
 * @param variantEnd End position (1-based, inclusive) of the given {@code altAllele}.  Must be > 0.
 * @param exonPositionList List of exon start / stop positions to cross-reference with the given {@code altAllele}.  Must not be {@code null}.
 * @param strand The {@link Strand} on which the {@code altAllele} is located.  Must not be {@code null}.  Must not be {@link Strand#NONE}.
 * @return A {@link SimpleInterval} containing the extents of the exon that overlaps the given altAllele.  {@code null} if no overlap occurs.
 */
public static SimpleInterval getOverlappingExonPositions(final Allele refAllele,
                                                         final Allele altAllele,
                                                         final String contig,
                                                         final int variantStart,
                                                         final int variantEnd,
                                                         final Strand strand,
                                                         final List<? extends htsjdk.samtools.util.Locatable> exonPositionList) {
    Utils.nonNull( refAllele );
    Utils.nonNull( altAllele );
    Utils.nonNull( contig );
    Utils.nonNull( exonPositionList );
    assertValidStrand( strand );

    ParamUtils.isPositive( variantStart, "Genome positions must be > 0." );
    ParamUtils.isPositive( variantEnd, "Genome positions must be > 0." );

    // Get the correct start / end positions:
    int alleleStart = variantStart;
    int alleleEnd   = variantEnd;

    if ( refAllele.length() > refAllele.length() ) {
        final int lengthAdjustment = (refAllele.length() - altAllele.length());
        if ( strand == Strand.NEGATIVE ) {
            alleleStart -= lengthAdjustment;
        }
        else {
            alleleEnd   += lengthAdjustment;
        }
    }

    // Create an interval so we can check for overlap:
    final SimpleInterval variantInterval = new SimpleInterval( contig, alleleStart, alleleEnd );

    int exonStart = -1;
    int exonStop  = -1;

    // Check for overlap:
    for ( final Locatable exon : exonPositionList ) {
        if ( variantInterval.overlaps(exon) ) {
            exonStart = exon.getStart();
            exonStop  = exon.getEnd();
        }
    }

    // Since we set the start/stop together we only need to check one of these:
    if ( exonStart == -1 ) {
        return null;
    }
    else {
        return new SimpleInterval( contig, exonStart, exonStop );
    }
}
 

/**
 * Check whether two locatables overlap.
 * <p>
 *    Two locatables overlap if the share the same contig and they have at least one
 *    base in common based on their start and end positions.
 * </p>
 * <p>
 *    This method returns {@code false} if either input {@link Locatable} has a {@code null}
 *    contig.
 * </p>
 *
 * @param left first locatable.
 * @param right second locatable.
 * @throws IllegalArgumentException if either {@code left} or {@code right} locatable
 *  is {@code null}.
 * @return {@code true} iff there is an overlap between both locatables.
 */
public static boolean overlaps(final Locatable left, final Locatable right) {
    Utils.nonNull(left,"the left locatable is null");
    Utils.nonNull(right,"the right locatable is null");
    if (left.getContig() == null || right.getContig() == null) {
        return false;
    } else if (!left.getContig().equals(right.getContig())) {
        return false;
    } else {
        return left.getStart() <= right.getEnd() && right.getStart() <= left.getEnd();
    }
}
 

/**
 * Create a new SimpleInterval from a {@link Locatable}
 * @param locatable any Locatable
 * @throws IllegalArgumentException if locatable violates any of the SimpleInterval constraints or is null
 */
public SimpleInterval(final Locatable locatable){
    this(Utils.nonNull(locatable).getContig(),
            locatable.getStart(), locatable.getEnd());
}