Java Code Examples for htsjdk.samtools.util.Interval#getStart()

/**
 * Converts an interval in SimpleInterval format into an htsjdk QueryInterval.
 *
 * In doing so, a header lookup is performed to convert from contig name to index
 *
 * @param interval interval to convert
 * @param sequenceDictionary sequence dictionary used to perform the conversion
 * @return an equivalent interval in QueryInterval format
 */
private static QueryInterval convertSimpleIntervalToQueryInterval( final Interval interval,	final SAMSequenceDictionary sequenceDictionary ) {
	if (interval == null) {
		throw new IllegalArgumentException("interval may not be null");
	}
	if (sequenceDictionary == null) {
		throw new IllegalArgumentException("sequence dictionary may not be null");
	}

	final int contigIndex = sequenceDictionary.getSequenceIndex(interval.getContig());
	if ( contigIndex == -1 ) {
		throw new IllegalArgumentException("Contig " + interval.getContig() + " not present in reads sequence " +
				"dictionary");
	}

	return new QueryInterval(contigIndex, interval.getStart(), interval.getEnd());
}
 

public static int compare (final Interval i1, final Interval i2, final SAMSequenceDictionary dict) {
 	int result = 0;
 	// if there's a sequence dictionary, compare on the index of the sequence instead of the contig name.
 	if (dict!=null) {
 		int seqIdx1 = dict.getSequenceIndex(i1.getContig());
 		int seqIdx2 = dict.getSequenceIndex(i2.getContig());
 		result = seqIdx1 - seqIdx2;
 	} else
result = i1.getContig().compareTo(i2.getContig());
 	// same as Interval.compareTo
 	if (result == 0) {
         if (i1.getStart() == i2.getStart())
	result = i1.getEnd() - i2.getEnd();
else
	result = i1.getStart() - i2.getStart();
         // added bonus, sort on interval names to tie break, if both intervals have names.
         if (result==0) {
         	String n1 = i1.getName();
         	String n2 = i2.getName();
         	if (n1!=null && n2!=null)
		result = n1.compareTo(n2);
         }
     }
 	return result;
 }
 

@Override
List<Bait> design(final BaitDesigner designer, final Interval target, final ReferenceSequence reference) {
    final List<Bait> baits = new LinkedList<Bait>();
    final int baitSize = designer.BAIT_SIZE;
    final int baitOffset = designer.BAIT_OFFSET;
    final int lastPossibleBaitStart = Math.min(target.getEnd(), reference.length() - baitSize);
    final int baitCount = 1 + (int) Math.floor((lastPossibleBaitStart - target.getStart()) / (double) baitOffset);

    int i = 0;
    for (int start = target.getStart(); start < lastPossibleBaitStart; start += baitOffset) {
        final Bait bait = new Bait(target.getContig(),
                start,
                CoordMath.getEnd(start, baitSize),
                target.isNegativeStrand(),
                designer.makeBaitName(target.getName(), ++i, baitCount));
        bait.addBases(reference, designer.DESIGN_ON_TARGET_STRAND);
        baits.add(bait);
    }
    return baits;
}
 

/** Emits a file with per base coverage if an output file has been set. */
private void emitPerBaseCoverageIfRequested() {
    if (this.perBaseOutput == null) return;

    final PrintWriter out = new PrintWriter(IOUtil.openFileForBufferedWriting(this.perBaseOutput));
    out.println("chrom\tpos\ttarget\tcoverage");
    for (final Map.Entry<Interval,Coverage> entry : this.highQualityCoverageByTarget.entrySet()) {
        final Interval interval = entry.getKey();
        final String chrom = interval.getContig();
        final int firstBase = interval.getStart();

        final int[] cov = entry.getValue().getDepths();
        for (int i = 0; i < cov.length; ++i) {
            out.print(chrom);
            out.print('\t');
            out.print(firstBase + i);
            out.print('\t');
            out.print(interval.getName());
            out.print('\t');
            out.print(cov[i]);
            out.println();
        }
    }

    out.close();
}
 

/**
 * Tests whether the variant is within one of the pseudo-autosomal regions
 */
private boolean isInPseudoAutosomalRegion(final String chr, final int pos) {
    for (final Interval par : parIntervals) {
        if (par.getContig().equals(chr) && pos >= par.getStart() && pos <= par.getEnd()) return true;
    }
    return false;
}
 

private static IntervalList generateBins(final IntervalList preparedIntervalList, final int binLength, final SAMSequenceDictionary sequenceDictionary) {
    if (binLength == 0) {
        return IntervalList.copyOf(preparedIntervalList);
    }
    final IntervalList bins = new IntervalList(sequenceDictionary);
    for (final Interval interval : preparedIntervalList) {
        for (int binStart = interval.getStart(); binStart <= interval.getEnd(); binStart += binLength) {
            final int binEnd = FastMath.min(binStart + binLength - 1, interval.getEnd());
            bins.add(new Interval(interval.getContig(), binStart, binEnd));
        }
    }
    return bins;
}
 

private static Interval lookupIntervalContainingLocus(final IntervalList source, final String chromosome, final int startIndex) {
    for (final Interval interval : source) {
        if (interval.getContig().equals(chromosome) && startIndex >= interval.getStart() && startIndex <= interval.getEnd()) {
            return interval;
        }
    }
    throw new NoSuchElementException();
}
 

protected static VariantContextBuilder reverseComplementVariantContext(final VariantContext source, final Interval target, final ReferenceSequence refSeq) {
    if (target.isPositiveStrand()) {
        throw new IllegalArgumentException("This should only be called for negative strand liftovers");
    }
    final int start = target.getStart();
    final int stop = target.getEnd();

    final List<Allele> origAlleles = new ArrayList<>(source.getAlleles());
    final VariantContextBuilder vcb = new VariantContextBuilder(source);

    vcb.rmAttribute(VCFConstants.END_KEY)
            .chr(target.getContig())
            .start(start)
            .stop(stop)
            .alleles(reverseComplementAlleles(origAlleles));

    if (isIndelForLiftover(source)) {
        // check that the reverse complemented bases match the new reference
        if (referenceAlleleDiffersFromReferenceForIndel(vcb.getAlleles(), refSeq, start, stop)) {
            return null;
        }
        leftAlignVariant(vcb, start, stop, vcb.getAlleles(), refSeq);
    }

    vcb.genotypes(fixGenotypes(source.getGenotypes(), origAlleles, vcb.getAlleles()));

    return vcb;
}
 

@Override
public List<Interval> takeSome(final Interval interval, final long idealSplitWeight, final long currentSize, final double projectSizeOfRemaining) {
    final long amount = idealSplitWeight - currentSize;

    if (amount >= interval.length()) {
        return CollectionUtil.makeList(interval, null);
    }

    if (amount == 0) {
        return CollectionUtil.makeList(null, interval);
    }

    final Interval left = new Interval(
            interval.getContig(),
            interval.getStart(),
            interval.getStart() + (int) amount - 1,
            interval.isNegativeStrand(),
            interval.getName()
    );
    final Interval right = new Interval(
            interval.getContig(),
            interval.getStart() + (int) amount,
            interval.getEnd(),
            interval.isNegativeStrand(),
            interval.getName()
    );
    return CollectionUtil.makeList(left, right);
}
 

private List<GTFRecord> getGTFRecordsFromIntronIntervals (List<Interval> introns, GeneFromGTF.TranscriptFromGTF t) {
	List<GTFRecord> result = new ArrayList<>();
	GeneFromGTF g = t.getGene();
	for (Interval i: introns) {
		@SuppressWarnings("deprecation") GTFRecord intronRecord = new GTFRecord(g.getContig(), i.getStart(), i.getEnd(), g.isNegativeStrand(),
                   g.getGeneID(), g.getName(), t.getTranscriptName(), t.getTranscriptID(), t.getTranscriptType(),
                   INTRON_FEATURE_TYPE, g.getGeneVersion());
		result.add(intronRecord);
	}		
	return (result);
}
 

private void writeSequence (final ReferenceSequence rs, final List<Interval> intervalsToMask, final FastaSequenceFileWriter writer) {
	String sequence = rs.getBaseString();
	if (intervalsToMask!=null && intervalsToMask.size()>0) {
		char [] seqArray = sequence.toCharArray();
		for (Interval i: intervalsToMask)
			for (int pos=i.getStart()-1; pos<i.getEnd(); pos++)
				seqArray[pos]='N';
		sequence=new String (seqArray);
	}
	writer.writeSequence(rs.getName(), sequence);
}
 

public static byte [] setSequenceToN (final byte [] fastaRefBases, final Interval interval) {
	byte [] result = fastaRefBases;
	int startBase=interval.getStart();
	int endBase=interval.getEnd();
	// the byte [] is base 0, the coordinates are base 1.
	Arrays.fill(result, startBase-1, endBase, StringUtil.charToByte('N'));
	return (result);
}
 

public static String getSequence (final byte [] fastaRefBases, final Interval interval) {
	int startBase=interval.getStart();
	int endBase=interval.getEnd();
	byte [] bases=getSubArray (fastaRefBases, startBase-1, endBase-1);
	StringBuilder b= new StringBuilder();
	String baseString = new String(bases);
	b.append(baseString);
	return (b.toString());
}
 

/** Returns true if the position is within the start-end range inclusive of the given interval. */
private boolean overlapsInterval(final int pos, final Interval interval) {
    return pos >= interval.getStart() && pos <= interval.getEnd();
}
 

private static SortedMap<String, List<GenomeRange>> performLiftOver(List<GenomeRange> inRange) {

    SortedMap<String, List<GenomeRange>> toReturn = new TreeMap<>();

    LiftOver liftOverTool = new LiftOver(new File(ARGS.liftOverChainPath));

    int all = inRange.size();
    int current = 0;
    int currentMarker = 1;

    for (GenomeRange currRange : inRange) {
      int oneBasedClosedStart = (int) currRange.getStart() + 1;
      int oneBasedClosedEnd = (int) currRange.getEnd();
      if ((long) oneBasedClosedStart != currRange.getStart() + 1
          || (long) oneBasedClosedEnd != currRange.getEnd()) {
        throw new IllegalArgumentException("performLiftOver: cannot specify ranges which are "
            + "greater than Integer.MAX_VALUE");
      }

      // These intervals are 1-based, close ended
      final Interval currInterval = new Interval(currRange.getChromosome(), oneBasedClosedStart,
          oneBasedClosedEnd, !currRange.isPositiveStrand(), currRange.getName());

      Interval liftedInterval = liftOverTool.liftOver(currInterval, ARGS.minMatch);
      if (liftedInterval == null) {
        logger.log(Level.FINE, "failed to liftover an interval");
      } else {
        // Change from one based to 0 based
        GenomeRange toAdd = new GenomeRange(liftedInterval.getSequence(),
            liftedInterval.getStart() - 1, liftedInterval.getEnd(), liftedInterval.getName(),
            liftedInterval.isPositiveStrand());
        if (!toReturn.containsKey(liftedInterval.getSequence())) {
          toReturn.put(liftedInterval.getSequence(), new ArrayList<GenomeRange>());
        }
        toReturn.get(liftedInterval.getSequence()).add(toAdd);
      }

      if ((100 * current++) / all > currentMarker) {
        logger.log(Level.INFO, String.format("performing liftover: at %d%%", currentMarker++));
      }
    }

    return toReturn;
  }
 

@Test(enabled=true)
public void testGeneSNPSplitting () {
	List<String> cellBarcodeList = ParseBarcodeFile.readCellBarcodeFile(cellBCFile);
	IntervalList loci = IntervalList.fromFile(snpIntervals);
	// make a new SNP that is 10 bases before the first SNP.  Reads in the test data will hit both.
	Interval i1 = loci.getIntervals().iterator().next();
	Interval i2 = new Interval (i1.getContig(), i1.getStart()-10, i1.getStart()-10, true, "testSNP1");
	loci.add(i2);



	//A read that hits only the default SNP: start read at 76227020
	SAMRecord r1 = getReadByName("NS500217:67:H14GMBGXX:1:11308:22039:11268", smallBAMFile);
	Collection<SAMRecord> tempList = processOneRead(r1, loci, cellBarcodeList);
	Assert.assertTrue(tempList.size()==1);
	Assert.assertEquals(1, getNumGenes(tempList));
	Assert.assertEquals(1, getNumSNPs(tempList));

	//A read that hits both the default SNP and the new SNP: start read at 76227000
	r1 = getReadByName("NS500217:67:H14GMBGXX:1:13202:10555:15929", smallBAMFile);
	tempList = processOneRead(r1, loci, cellBarcodeList);
	Assert.assertTrue(tempList.size()==2);
	Assert.assertEquals(1, getNumGenes(tempList));
	Assert.assertEquals(2, getNumSNPs(tempList));


	//  A read that hits only the default SNP: start read at 76227020, has 2 genes.
	r1 = getReadByName("NS500217:67:H14GMBGXX:1:21206:20467:19854", smallBAMFile);
	tempList = processOneRead(r1, loci, cellBarcodeList);
	Assert.assertTrue(tempList.size()==2);
	Assert.assertEquals(2, getNumGenes(tempList));
	Assert.assertEquals(1, getNumSNPs(tempList));


	//A read that hits both the default SNP and the new SNP: start read at 76227000, has 2 genes.
	r1 = getReadByName("NS500217:67:H14GMBGXX:1:22302:3826:3320", smallBAMFile);
	tempList = processOneRead(r1, loci, cellBarcodeList);
	Assert.assertTrue(tempList.size()==4);
	Assert.assertEquals(2, getNumGenes(tempList));
	Assert.assertEquals(2, getNumSNPs(tempList));


}
 

/**
 * Converts an Interval object into a string representation that can be parsed by fromString.
 * Not using the Interval.toString method, because it has tabs, which might really screw up something else...
 * Instead, use ENCODE_DELIMITER as a delimiter for all fields except the range fields start-end.
 * @param i The interval to parse
 * @return A string representation of the interval.
 */
public static String toString (final Interval i) {
	// chr1:1-10	+	foo
	String result = i.getContig() + ENCODE_DELIMITER + i.getStart() + "-" + i.getEnd() + ENCODE_DELIMITER + (i.isNegativeStrand() ? '-' : '+') + ENCODE_DELIMITER + ((null == i.getName()) ? '.' : i.getName());
	return result;
}