private boolean isAfter(BreakendSummary breakendSummary, SAMRecord position) { return position.getReferenceIndex() > breakendSummary.referenceIndex || (position.getReferenceIndex() == breakendSummary.referenceIndex && position.getUnclippedStart() > breakendSummary.end + windowSize); } @Override
public SgNode(LinearGenomicCoordinate lgc, Read read, int offset) { this.inferredPosition = lgc.getLinearCoordinate(read.getRead().getReferenceIndex(), read.getRead().getUnclippedStart()) + offset; this.read = read; } /**
/** * Determine whether the two reads are the same. This is useful for * filtering out read pairs which, due to sequencing chemistry error, both * sequence the same end of the fragment. * * @param r1 * read to compare * @param r2 * read to compare * @return true if the reads appear to be the same, false otherwise */ public static boolean areSameRead(SAMRecord r1, SAMRecord r2) { // TODO: compare read sequences return !r1.getReadUnmappedFlag() && !r2.getReadUnmappedFlag() && r1.getReferenceIndex().equals(r2.getReferenceIndex()) && (r1.getUnclippedStart() == r2.getUnclippedStart() || r1.getUnclippedEnd() == r2.getUnclippedEnd()) && r1.getReadNegativeStrandFlag() == r2.getReadNegativeStrandFlag(); }
public int getSoftUnclippedAlignmentStart(SAMRecord rec){ List<CigarElement> cigar= rec.getCigar().getCigarElements(); if(cigar.size() == 0){ return rec.getUnclippedStart(); } if(cigar.size() == 1 && cigar.get(0).getOperator().equals(CigarOperator.SOFT_CLIP)){ return rec.getUnclippedStart(); } int offset= 0; for(int i= 0; i < cigar.size(); i++){ CigarElement op = cigar.get(i); if(op.getOperator().equals(CigarOperator.HARD_CLIP)){ continue; } else if(op.getOperator().equals(CigarOperator.SOFT_CLIP)){ offset += op.getLength(); } else { break; } } int start= rec.getAlignmentStart() - offset; // if(start < 1){ // start= 1; // } return start; }
public Alignment buildAlignment(SAMRecord record, boolean compareReference) { if(compareReference) { String seq; try { seq = adaptor.getSequence(new Region(record.getReferenceName(), record.getUnclippedStart(), record.getUnclippedEnd())); } catch (IOException e) { e.printStackTrace(); return null; } return buildAlignment(record, seq); } else { return buildAlignment(record); } }
@Requires("r1 != null && r2 != null") @Ensures("result == 0 || result == 1 || result == -1") public int compare(SAMRecord r1, SAMRecord r2) { int result; if (r1 == r2) result = 0; else if (r1.getReadUnmappedFlag()) result = 1; else if (r2.getReadUnmappedFlag()) result = -1; else { final int cmpContig = compareContigs(r1.getReferenceIndex(), r2.getReferenceIndex()); if (cmpContig != 0) result = cmpContig; else { if (r1.getUnclippedStart() < r2.getUnclippedStart()) result = -1; else result = 1; } } return result; } }
/** * Estimates the size of sequenced fragment * */ public static int calculateFragmentSize(SAMRecord record1, SAMRecord record2, PairOrientation expectedOrientation) { if (expectedOrientation != PairOrientation.FR) throw new RuntimeException("NYI"); if (record1.getReadUnmappedFlag() || record2.getReadUnmappedFlag() || !record1.getReferenceIndex().equals(record2.getReferenceIndex()) || // FR assumption record1.getReadNegativeStrandFlag() == record2.getReadNegativeStrandFlag()) { return 0; } // Assuming FR orientation and adapter sequences have been removed if (record1.getReadNegativeStrandFlag()) { // <--record int r1end = record1.getUnclippedEnd(); int r2start = record2.getUnclippedStart(); return r1end - r2start + 1; } else { int r1start = record1.getUnclippedStart(); int r2end = record2.getUnclippedEnd(); return r2end - r1start + 1; } }
/** * Create a genome loc, given a read using its unclipped alignment. If the read is unmapped, *and* yet the read has a contig and start position, * then a GenomeLoc is returned for contig:start-start, otherwise an UNMAPPED GenomeLoc is returned. * * @param read the read from which to create a genome loc * * @return the GenomeLoc that was created */ @Requires("read != null") @Ensures("result != null") public GenomeLoc createGenomeLocUnclipped(final SAMRecord read) { if ( read.getReadUnmappedFlag() && read.getReferenceIndex() == -1 ) // read is unmapped and not placed anywhere on the genome return GenomeLoc.UNMAPPED; else { // Use Math.max to ensure that end >= start (Picard assigns the end to reads that are entirely within an insertion as start-1) final int end = read.getReadUnmappedFlag() ? read.getUnclippedEnd() : Math.max(read.getUnclippedEnd(), read.getUnclippedStart()); return createGenomeLoc(read.getReferenceName(), read.getReferenceIndex(), read.getUnclippedStart(), end, false); } }
private long getWindowPosition(Read r) { // TODO: update to inferred position interval return lgc.getLinearCoordinate(r.getRead().getReferenceIndex(), (r.getRead().getUnclippedEnd() + r.getRead().getUnclippedStart()) / 2); } private long currentTransformPosition() { return currentPosition - windowSize; }
/** * Populates the set of transient attributes on SAMRecords if they are not already there. */ private void populateTransientAttributes(final SAMRecord... recs) { for (final SAMRecord rec : recs) { if (rec.getTransientAttribute(Attr.LibraryId) != null) continue; rec.setTransientAttribute(Attr.LibraryId, getLibraryId(rec)); rec.setTransientAttribute(Attr.ReadCoordinate, rec.getReadNegativeStrandFlag() ? rec.getUnclippedEnd() : rec.getUnclippedStart()); rec.setTransientAttribute(Attr.MateCoordinate, getMateCoordinate(rec)); } }
/** * Populates the set of transient attributes on SAMRecords if they are not already there. */ private void populateTransientAttributes(final SAMRecord... recs) { for (final SAMRecord rec : recs) { if (rec.getTransientAttribute(Attr.LibraryId) != null) continue; rec.setTransientAttribute(Attr.LibraryId, getLibraryId(rec)); rec.setTransientAttribute(Attr.ReadCoordinate, rec.getReadNegativeStrandFlag() ? rec.getUnclippedEnd() : rec.getUnclippedStart()); rec.setTransientAttribute(Attr.MateCoordinate, getMateCoordinate(rec)); } }
positionClosestToBreakpoint = local.getAlignmentStart(); intervalDirection = -1; intervalExtendedReadDueToLocalClipping = local.getAlignmentStart() - local.getUnclippedStart(); intervalReducedDueToRemoteMapping -= remote.getUnclippedEnd() - remote.getAlignmentEnd(); } else { intervalReducedDueToRemoteMapping -= remote.getAlignmentStart() - remote.getUnclippedStart();
public BAQCalculationResult calcBAQFromHMM(final SAMRecord read, final ReferenceSequenceFile refReader) { // start is alignment start - band width / 2 - size of first I element, if there is one. Stop is similar int offset = getBandWidth() / 2; long readStart = includeClippedBases ? read.getUnclippedStart() : read.getAlignmentStart(); long start = Math.max(readStart - offset - ReadUtils.getFirstInsertionOffset(read), 1); long stop = (includeClippedBases ? read.getUnclippedEnd() : read.getAlignmentEnd()) + offset + ReadUtils.getLastInsertionOffset(read); if ( stop > refReader.getSequenceDictionary().getSequence(read.getReferenceName()).getSequenceLength() ) { return null; } else { // now that we have the start and stop, get the reference sequence covering it ReferenceSequence refSeq = refReader.getSubsequenceAt(read.getReferenceName(), start, stop); return calcBAQFromHMM(read, refSeq.getBases(), (int)(start - readStart)); } }
startPosition = local.getUnclippedStart() + minFragSize - remote.getReadLength(); endPosition = local.getUnclippedStart() + maxFragSize - remote.getReadLength(); } else {
final int read5PrimeStart = (rec.getReadNegativeStrandFlag()) ? rec.getUnclippedEnd() : rec.getUnclippedStart(); final int mate5PrimeStart = (rec.getMateNegativeStrandFlag()) ? SAMUtils.getMateUnclippedEnd(rec) : SAMUtils.getMateUnclippedStart(rec);
final int read5PrimeStart = (rec.getReadNegativeStrandFlag()) ? rec.getUnclippedEnd() : rec.getUnclippedStart(); final int mate5PrimeStart = (rec.getMateNegativeStrandFlag()) ? SAMUtils.getMateUnclippedEnd(rec) : SAMUtils.getMateUnclippedStart(rec);
return r1end - r2start + 1; } else { int r1start = record.getUnclippedStart(); int r2end = mc == null ?
public static Alignment buildAlignment(SAMRecord record, Map<String, Object> attributes, String referenceSequence) { List<Alignment.AlignmentDifference> differences; differences = AlignmentUtils.getDifferencesFromCigar(record, referenceSequence, Integer.MAX_VALUE); Alignment alignment = new Alignment(record.getReadName(), record.getReferenceName(), record.getAlignmentStart(), record.getAlignmentEnd(), record.getUnclippedStart(), record.getUnclippedEnd(), record.getReadLength(), record.getMappingQuality(), record.getBaseQualityString(),//.replace("\\", "\\\\").replace("\"", "\\\""), record.getMateReferenceName(), record.getMateAlignmentStart(), record.getInferredInsertSize(), record.getFlags(), differences, attributes); return alignment; } public static Alignment buildAlignment(SAMRecord record, String referenceSequence) {
public VariantEvidence(int k, NonReferenceReadPair pair, LinearGenomicCoordinate lgc) { this.evidenceID = pair.getEvidenceID(); boolean shouldReverseComplement = !pair.onExpectedStrand(); byte[] bases = pair.getNonReferenceRead().getReadBases(); byte[] quals = pair.getNonReferenceRead().getBaseQualities(); int chrPos; if (pair.getBreakendSummary().direction == BreakendDirection.Forward) { // Assumes FR orientation chrPos = pair.getLocalledMappedRead().getUnclippedStart() + pair.getEvidenceSource().getExpectedFragmentSize() - bases.length; } else { chrPos = pair.getLocalledMappedRead().getUnclippedEnd() - pair.getEvidenceSource().getExpectedFragmentSize() + 1; } this.expectedAnchorPos = lgc.getLinearCoordinate(pair.getBreakendSummary().referenceIndex, chrPos); this.referenceKmerStartOffset = -1; this.referenceKmerEndOffset = -1; this.startSkipKmerCount = 0; this.endSkipKmerCount = 0; this.isExact = pair.isBreakendExact(); this.be = pair.getBreakendSummary(); this.category = pair.getEvidenceSource().getSourceCategory(); this.kmers = new PackedKmerList(k, bases, quals, shouldReverseComplement, shouldReverseComplement); this.ambiguous = markAmbiguous(k, bases); } public VariantEvidence(int k, SingleReadEvidence softClipEvidence, LinearGenomicCoordinate lgc) {
this.read1Coordinate = record.getReadNegativeStrandFlag() ? record.getUnclippedEnd() : record.getUnclippedStart(); if (record.getReadUnmappedFlag()) { throw new PicardException("Found an unexpected unmapped read");