static private void elementStraddlesClippedRead(List<CigarElement> newCigar, CigarElement c, int relativeClippedPosition, int clippedBases){ final CigarOperator op = c.getOperator(); int clipAmount = clippedBases; if (op.consumesReadBases()){ if (op.consumesReferenceBases() & relativeClippedPosition > 0){ newCigar.add(new CigarElement(relativeClippedPosition, op)); } if (!op.consumesReferenceBases()){ clipAmount = clippedBases + relativeClippedPosition; } } else if (relativeClippedPosition != 0){ throw new SAMException("Unexpected non-0 relativeClippedPosition " + relativeClippedPosition); } newCigar.add(new CigarElement(clipAmount, CigarOperator.S)); // S is always last element }
static private void elementStraddlesClippedRead(List<CigarElement> newCigar, CigarElement c, int relativeClippedPosition, int clippedBases){ final CigarOperator op = c.getOperator(); int clipAmount = clippedBases; if (op.consumesReadBases()){ if (op.consumesReferenceBases() & relativeClippedPosition > 0){ newCigar.add(new CigarElement(relativeClippedPosition, op)); } if (!op.consumesReferenceBases()){ clipAmount = clippedBases + relativeClippedPosition; } } else if (relativeClippedPosition != 0){ throw new SAMException("Unexpected non-0 relativeClippedPosition " + relativeClippedPosition); } newCigar.add(new CigarElement(clipAmount, CigarOperator.S)); // S is always last element }
static private void elementStraddlesClippedRead(List<CigarElement> newCigar, CigarElement c, int relativeClippedPosition, int clippedBases){ final CigarOperator op = c.getOperator(); int clipAmount = clippedBases; if (op.consumesReadBases()){ if (op.consumesReferenceBases() & relativeClippedPosition > 0){ newCigar.add(new CigarElement(relativeClippedPosition, op)); } if (!op.consumesReferenceBases()){ clipAmount = clippedBases + relativeClippedPosition; } } else if (relativeClippedPosition != 0){ throw new SAMException("Unexpected non-0 relativeClippedPosition " + relativeClippedPosition); } newCigar.add(new CigarElement(clipAmount, CigarOperator.S)); // S is always last element }
/** * Returns the number of reference bases 'consumed' by the given cigar * @param list * @return base count */ public static int referenceLength(List<CigarElement> list) { int length = 0; for (CigarElement e : list) { if (e.getOperator().consumesReferenceBases()) { length += e.getLength(); } } return length; } /**
/** * Check if cigar start with an insertion, ignoring other operators that do not consume references bases * * @param cigar the cigar * @return <code>true</code> if the first operator to consume reference bases or be an insertion, is an insertion; <code>false</code> otherwise */ protected static boolean startWithInsertion(final Cigar cigar) { for (final CigarElement element : cigar.getCigarElements()) { if (element.getOperator() == CigarOperator.I) return true; if (!element.getOperator().consumesReferenceBases()) continue; break; } return false; }
/** * Check if cigar start with an insertion, ignoring other operators that do not consume references bases * * @param cigar the cigar * @return <code>true</code> if the first operator to consume reference bases or be an insertion, is an insertion; <code>false</code> otherwise */ protected static boolean startWithInsertion(final Cigar cigar) { for (final CigarElement element : cigar.getCigarElements()) { if (element.getOperator() == CigarOperator.I) return true; if (!element.getOperator().consumesReferenceBases()) continue; break; } return false; }
/** * Determines if a cigar has any element that both consumes read bases and consumes reference bases * (e.g. is not all soft-clipped) */ public static boolean cigarMapsNoBasesToRef(final Cigar cigar) { for (final CigarElement el : cigar.getCigarElements()) { if (el.getOperator().consumesReadBases() && el.getOperator().consumesReferenceBases()) { return false; } } return true; }
/** * Determines if a cigar has any element that both consumes read bases and consumes reference bases * (e.g. is not all soft-clipped) */ public static boolean cigarMapsNoBasesToRef(final Cigar cigar) { for (final CigarElement el : cigar.getCigarElements()) { if (el.getOperator().consumesReadBases() && el.getOperator().consumesReferenceBases()) { return false; } } return true; }
public static int countMappedBases(List<CigarElement> cigar) { int mapped = 0; for (CigarElement e : cigar) { if (e.getOperator().consumesReferenceBases() && e.getOperator().consumesReadBases()) { mapped += e.getLength(); } } return mapped; } /**
/** * Determines if a cigar has any element that both consumes read bases and consumes reference bases * (e.g. is not all soft-clipped) */ public static boolean cigarMapsNoBasesToRef(final Cigar cigar) { for (final CigarElement el : cigar.getCigarElements()) { if (el.getOperator().consumesReadBases() && el.getOperator().consumesReferenceBases()) { return false; } } return true; }
/** * Determines if a cigar has any element that both consumes read bases and consumes reference bases * (e.g. is not all soft-clipped) */ public static boolean cigarMapsNoBasesToRef(final Cigar cigar) { for (final CigarElement el : cigar.getCigarElements()) { if (el.getOperator().consumesReadBases() && el.getOperator().consumesReferenceBases()) { return false; } } return true; }
/** * Calculates the reference coordinate for a read coordinate * * @param read the read * @param offset the base in the read (coordinate in the read) * @return the reference coordinate correspondent to this base */ public static long getReferenceCoordinateForReadCoordinate(GATKSAMRecord read, int offset) { if (offset > read.getReadLength()) throw new ReviewedGATKException(String.format(OFFSET_OUT_OF_BOUNDS_EXCEPTION, offset, read.getReadLength())); long location = read.getAlignmentStart(); Iterator<CigarElement> cigarElementIterator = read.getCigar().getCigarElements().iterator(); while (offset > 0 && cigarElementIterator.hasNext()) { CigarElement cigarElement = cigarElementIterator.next(); long move = 0; if (cigarElement.getOperator().consumesReferenceBases()) move = (long) Math.min(cigarElement.getLength(), offset); location += move; offset -= move; } if (offset > 0 && !cigarElementIterator.hasNext()) throw new ReviewedGATKException(OFFSET_NOT_ZERO_EXCEPTION); return location; }
/** * Returns the alignment offset of the given base relative to the starting alignment * @param cigar alignment CIGAR * @param readBaseOffset base offset * @return offset relative to first alignment */ public static int offsetOf(Cigar cigar, int readBaseOffset) { List<CigarElement> cl = cigar.getCigarElements(); int basesLeft = readBaseOffset; int currentAlignmentOffset = 0; for (int i = 0; i < cl.size(); i++) { CigarElement ce = cl.get(i); if (ce.getOperator().consumesReadBases() && ce.getOperator().consumesReferenceBases()) { if (basesLeft < ce.getLength()) { return currentAlignmentOffset + basesLeft; } } if (ce.getOperator().consumesReferenceBases()) { currentAlignmentOffset += ce.getLength(); } if (ce.getOperator().consumesReadBases()) { basesLeft = Math.max(0, basesLeft - ce.getLength()); } } throw new IllegalArgumentException(String.format("Offset of %d base not defined for %s", readBaseOffset, cigar)); } /**
public static List<SimpleBEDFeature> getSpanningDeletion(final SAMRecord r, final int minMapq, final int minIndelSize, final int maxMappedBases, final int minIndelComponentSize) { List<SimpleBEDFeature> list = new ArrayList<SimpleBEDFeature>(); if (r.getMappingQuality() < minMapq) return list; List<CigarElement> ce = r.getCigar().getCigarElements(); int offset = 0; for (int i = 0; i < ce.size(); i++) { CigarElement e = ce.get(i); if (e.getOperator() == CigarOperator.D) { if (maxMappedBases > 0) { throw new RuntimeException("Indel merging via MAX_INTERVENING_MAPPED_BASES not yet implemented."); } if (e.getLength() > minIndelSize) { list.add(new SimpleBEDFeature(r.getAlignmentStart() + offset, r.getAlignmentStart() + offset + e.getLength(), r.getReferenceName())); } } if (e.getOperator().consumesReferenceBases()) { offset += e.getLength(); } } return list; } public static SimpleBEDFeature getSplitReadDeletion(SAMRecord r, ChimericAlignment ca, double minLengthPortion, int minMapq, int minIndelSize) {
/** * Counts the number of bases mapped to any reference position in both alignment * @param cigar1 * @param cigar2 * @return number of bases mapped to reference in both alignments */ public static int commonReferenceBases(Cigar cigar1, Cigar cigar2) { int count = 0; int readLength = cigar1.getReadLength(); assert(readLength == cigar2.getReadLength()); CigarOperatorIterator it1 = new CigarOperatorIterator(cigar1); CigarOperatorIterator it2 = new CigarOperatorIterator(cigar2); for (int readBaseOffset = 0; readBaseOffset < readLength; readBaseOffset++) { CigarOperator op1; CigarOperator op2; do { op1 = it1.next(); } while (!op1.consumesReadBases()); do { op2 = it2.next(); } while (!op2.consumesReadBases()); if (op1.consumesReferenceBases() && op2.consumesReferenceBases()) { count++; } } return count; } public static int getEndClipLength(List<CigarElement> elements) {
private static ImmutableRangeSet<Integer> getMappedBases(int start, Cigar cigar) { Builder<Integer> builder = ImmutableRangeSet.builder(); int position = start; for (CigarElement op : cigar.getCigarElements()) { if (op.getOperator().consumesReferenceBases()) { if (op.getOperator().consumesReadBases()) { builder.add(Range.closedOpen(position, position + op.getLength())); } position += op.getLength(); } } return builder.build(); } /**
private int countDeletions(SAMRecord r, int minDeletionSize, int startPosition, int endPosition) { int position = r.getAlignmentStart(); int count = 0; for (CigarElement ce : r.getCigar().getCigarElements()) { if (position >= endPosition) return count; if (position >= startPosition && position <= endPosition) { if (ce.getOperator() == CigarOperator.D && ce.getLength() > minDeletionSize && position + ce.getLength() < endPosition) { // deletion contained within the interval in question count += ce.getLength(); } } if (ce.getOperator().consumesReferenceBases()) { position += ce.getLength(); } } return count; } }
private static void updateTracks(final List<SAMRecord> samRecords, final ReferenceTracks tracks) { for (final SAMRecord samRecord : samRecords) { if (samRecord.getAlignmentStart() != SAMRecord.NO_ALIGNMENT_START) { int refPos = samRecord.getAlignmentStart(); int readPos = 0; for (final CigarElement cigarElement : samRecord.getCigar().getCigarElements()) { if (cigarElement.getOperator().consumesReferenceBases()) { for (int elementIndex = 0; elementIndex < cigarElement.getLength(); elementIndex++) tracks.addCoverage(refPos + elementIndex, 1); } switch (cigarElement.getOperator()) { case M: case X: case EQ: for (int pos = readPos; pos < cigarElement.getLength(); pos++) { final byte readBase = samRecord.getReadBases()[readPos + pos]; final byte refBase = tracks.baseAt(refPos + pos); if (readBase != refBase) tracks.addMismatches(refPos + pos, 1); } break; default: break; } readPos += cigarElement.getOperator().consumesReadBases() ? cigarElement.getLength() : 0; refPos += cigarElement.getOperator().consumesReferenceBases() ? cigarElement.getLength() : 0; } } } }
/**List of positions on screen where the skipped bases (cigar op: N) start * @throws IOException * @throws InvalidGenomicCoordsException * */ private void setTextPositionsOfSkippedBases() throws InvalidGenomicCoordsException, IOException{ int genomicPosition= this.getSamRecord().getAlignmentStart(); List<CigarElement> cigar = this.getSamRecord().getCigar().getCigarElements(); for(CigarElement el : cigar){ if(el.getOperator().equals(CigarOperator.SKIPPED_REGION)){ int[] textPositions= new int[2]; // +1 because textPosition is 1-based textPositions[0]= Utils.getIndexOfclosestValue(genomicPosition, this.gc.getMapping()) + 1; textPositions[1]= Utils.getIndexOfclosestValue(genomicPosition + el.getLength(), this.gc.getMapping()) + 1; this.textPositionsOfSkippedBases.add(textPositions); }; if(el.getOperator().consumesReferenceBases()){ genomicPosition += el.getLength(); } } }
private static void processAnchor(RangeSet<Integer> anchoredBases, SAMRecord record) { if (record.getReadUnmappedFlag()) return; if (CigarUtil.widthOfImprecision(record.getCigar()) > 0) return; // unanchored int pos = record.getAlignmentStart(); for (CigarElement ce : record.getCigar().getCigarElements()) { switch (ce.getOperator()) { case M: case EQ: case X: anchoredBases.add(Range.closedOpen(pos, pos + ce.getLength())); break; default: break; } if (ce.getOperator().consumesReferenceBases()) { pos += ce.getLength(); } } } /**