private boolean isLastNonSoftClippedElem(int idx, Cigar cigar) { int numElems = cigar.getCigarElements().size(); // Last element, not soft clipped. if ((idx == numElems-1) && (cigar.getCigarElement(idx).getOperator() != CigarOperator.S)) { return true; } // Second to last element, with last element soft clipped. if ((idx == numElems-2) && (cigar.getCigarElement(numElems-1).getOperator() == CigarOperator.S)) { return true; } return false; }
private boolean isFirstNonSoftClippedElem(int idx, Cigar cigar) { // First element, not soft clipped if ((idx == 0) && (cigar.getCigarElement(0).getOperator() != CigarOperator.S)) { return true; } // Second element, with first element soft clipped if ((idx == 1) && (cigar.getCigarElement(0).getOperator() == CigarOperator.S)) { return true; } return false; }
/** * Checks to see if the provided Cigar could be considered the "sentinel cigar" that indicates * that the actual cigar is too long for the BAM spec and should be taken from the CG tag. This * was introduced in SAM v1.6. */ static boolean isSentinelCigar(final Cigar cigar, final int readLength) { // There's an implicit assumption here there readLength == length of read in cigar, unless readLength==0 return cigar.numCigarElements() == 2 && cigar.getCigarElement(1).getOperator() == CigarOperator.N && cigar.getCigarElement(0).getOperator() == CigarOperator.S && (cigar.getCigarElement(0).getLength() == readLength || readLength == 0) ; }
protected static Cigar reclipCigar(Cigar cigar, SAMRecord read) { ArrayList<CigarElement> elements = new ArrayList<CigarElement>(); int i = 0; int n = read.getCigar().numCigarElements(); while ( i < n && isClipOperator(read.getCigar().getCigarElement(i).getOperator()) ) elements.add(read.getCigar().getCigarElement(i++)); elements.addAll(cigar.getCigarElements()); i++; while ( i < n && !isClipOperator(read.getCigar().getCigarElement(i).getOperator()) ) i++; while ( i < n && isClipOperator(read.getCigar().getCigarElement(i).getOperator()) ) elements.add(read.getCigar().getCigarElement(i++)); return new Cigar(elements); }
/** * Checks to see if the provided Cigar could be considered the "sentinel cigar" that indicates * that the actual cigar is too long for the BAM spec and should be taken from the CG tag. This * was introduced in SAM v1.6. */ static boolean isSentinelCigar(final Cigar cigar, final int readLength) { // There's an implicit assumption here there readLength == length of read in cigar, unless readLength==0 return cigar.numCigarElements() == 2 && cigar.getCigarElement(1).getOperator() == CigarOperator.N && cigar.getCigarElement(0).getOperator() == CigarOperator.S && (cigar.getCigarElement(0).getLength() == readLength || readLength == 0) ; }
/** * Skip the insertions and deletions that are right after or before introns (they indicate of aligner problem) * @param cigar cigar of the read * @param ci current index of element in cigar * @return true if cigar operator must be skipped, false if not */ private static boolean skipIndelNextToIntron(Cigar cigar, int ci) { if ((cigar.numCigarElements() > ci && cigar.getCigarElement(ci + 1).getOperator() == CigarOperator.N) || (ci > 1 && cigar.getCigarElement(ci - 1).getOperator() == CigarOperator.N)) { return true; } return false; }
/** * Does cigar start or end with a deletion operation? * * @param cigar a non-null cigar to test * @return true if the first or last operator of cigar is a D */ public static boolean startsOrEndsWithInsertionOrDeletion(final Cigar cigar) { if ( cigar == null ) throw new IllegalArgumentException("Cigar cannot be null"); if ( cigar.isEmpty() ) return false; final CigarOperator first = cigar.getCigarElement(0).getOperator(); final CigarOperator last = cigar.getCigarElement(cigar.numCigarElements()-1).getOperator(); return first == CigarOperator.D || first == CigarOperator.I || last == CigarOperator.D || last == CigarOperator.I; }
static CigarOperator getCigarOperator(Cigar cigar, int ci) { CigarOperator operator = cigar.getCigarElement(ci).getOperator(); // Treat insertions at the edge as soft-clipping if ((ci == 0 || ci == cigar.numCigarElements() - 1) && operator == CigarOperator.I) { operator = CigarOperator.S; } return operator; }
private static boolean isBeforeOp(final Cigar cigar, final int currentOperatorIndex, final CigarOperator op, final boolean isLastOp, final boolean isLastBaseOfOp) { return !isLastOp && isLastBaseOfOp && cigar.getCigarElement(currentOperatorIndex+1).getOperator() == op; }
private static boolean isAfterOp(final Cigar cigar, final int currentOperatorIndex, final CigarOperator op, final boolean isFirstOp, final boolean isFirstBaseOfOp) { return !isFirstOp && isFirstBaseOfOp && cigar.getCigarElement(currentOperatorIndex-1).getOperator() == op; } }
/** * If a read starts in INSERTION, returns the first element length. * * Warning: If the read has Hard or Soft clips before the insertion this function will return 0. * * @param read * @return the length of the first insertion, or 0 if there is none (see warning). */ public final static int getFirstInsertionOffset(SAMRecord read) { CigarElement e = read.getCigar().getCigarElement(0); if ( e.getOperator() == CigarOperator.I ) return e.getLength(); else return 0; }
/** * Helper function to get the cigar element of the next or previous genomic position * @param direction the direction to look in * @return a CigarElement, or null if no such element exists */ @Ensures("result == null || ON_GENOME_OPERATORS.contains(result.getOperator())") private CigarElement getNeighboringOnGenomeCigarElement(final Direction direction) { final int increment = direction == Direction.NEXT ? 1 : -1; final int nCigarElements = read.getCigarLength(); for ( int i = currentCigarOffset + increment; i >= 0 && i < nCigarElements; i += increment) { final CigarElement elt = read.getCigar().getCigarElement(i); if ( ON_GENOME_OPERATORS.contains(elt.getOperator()) ) return elt; } // getting here means that you didn't find anything return null; }
public int getAlignmentStart() { int start = read.getAlignmentStart(); if (read.getCigarLength() > 0 && read.getCigar().getCigarElement(0).getOperator() == CigarOperator.I) { start = start - 1; } return start; }
private int getAlignmentStart(SAMRecord read) { int start = read.getAlignmentStart(); if (read.getCigarLength() > 0 && read.getCigar().getCigarElement(0).getOperator() == CigarOperator.I) { start = start - 1; } return start; }
/** * If a read ends in INSERTION, returns the last element length. * * Warning: If the read has Hard or Soft clips after the insertion this function will return 0. * * @param read * @return the length of the last insertion, or 0 if there is none (see warning). */ public final static int getLastInsertionOffset(SAMRecord read) { CigarElement e = read.getCigar().getCigarElement(read.getCigarLength() - 1); if ( e.getOperator() == CigarOperator.I ) return e.getLength(); else return 0; }
@Test(dataProvider = "AddCigarElementsData", enabled = !DEBUG) public void testAddCigarElements(final String cigarString, final int pos, final int start, final int end, final String expectedCigarString) { final Cigar cigar = TextCigarCodec.decode(cigarString); final CigarElement elt = cigar.getCigarElement(0); final Cigar expectedCigar = TextCigarCodec.decode(expectedCigarString); final List<CigarElement> elts = new LinkedList<CigarElement>(); final int actualEndPos = AlignmentUtils.addCigarElements(elts, pos, start, end, elt); Assert.assertEquals(actualEndPos, pos + elt.getLength()); Assert.assertEquals(AlignmentUtils.consolidateCigar(new Cigar(elts)), expectedCigar); }
private int maxSoftClipLength(SAMRecord read, Feature region) { int len = 0; if (read.getCigarLength() > 1) { CigarElement elem = read.getCigar().getCigarElement(0); if (elem.getOperator() == CigarOperator.S && read.getAlignmentStart() >= region.getStart()-readLength) { len = elem.getLength(); } elem = read.getCigar().getCigarElement(read.getCigarLength()-1); if (elem.getOperator() == CigarOperator.S && elem.getLength() > len && read.getAlignmentEnd() <= region.getEnd()+readLength) { len = elem.getLength(); } } return len; }
/** * Convert CIGAR from object representation to disk representation. * @return Array of unsigned ints, one for each element of CIGAR. */ int[] encode(final Cigar cigar) { if (cigar.numCigarElements() == 0) { return new int[0]; } // Binary rep can be no longer than 1/2 of text rep // Although this is documented as uint, I think lengths will never get that long, // and it's a pain in Java. final int[] binaryCigar = new int[cigar.numCigarElements()]; int binaryCigarLength = 0; for (int i = 0; i < cigar.numCigarElements(); ++i) { final CigarElement cigarElement = cigar.getCigarElement(i); final int op = CigarOperator.enumToBinary(cigarElement.getOperator()); binaryCigar[binaryCigarLength++] = cigarElement.getLength() << 4 | op; } return binaryCigar; }
@Test public void testRBPMappingQuals() { // create a read with high MQ final GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "read", 0, 1, 10); read.setReadBases(Utils.dupBytes((byte) 'A', 10)); read.setBaseQualities(Utils.dupBytes((byte) 30, 10)); read.setCigarString("10M"); read.setMappingQuality(200); // set a MQ higher than max signed byte // now create the RBP final List<PileupElement> elts = new LinkedList<>(); elts.add(new PileupElement(read, 0, read.getCigar().getCigarElement(0), 0, 0)); final Map<String, ReadBackedPileupImpl> pileupsBySample = new HashMap<>(); pileupsBySample.put("foo", new ReadBackedPileupImpl(loc, elts)); final ReadBackedPileup pileup = new ReadBackedPileupImpl(loc, pileupsBySample); Assert.assertEquals(pileup.getMappingQuals()[0], 200); } }
private List<PileupElement> makeReads( final int n, final int mapq, final String op ) { final int readLength = 3; final List<PileupElement> elts = new LinkedList<PileupElement>(); for ( int i = 0; i < n; i++ ) { GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "read", 0, 1, readLength); read.setReadBases(Utils.dupBytes((byte) 'A', readLength)); read.setBaseQualities(Utils.dupBytes((byte) 30, readLength)); read.setCigarString("1M1" + op + "1M"); read.setMappingQuality(mapq); final int baseOffset = op.equals("M") ? 1 : 0; final CigarElement cigarElement = read.getCigar().getCigarElement(1); elts.add(new PileupElement(read, baseOffset, cigarElement, 1, 0)); } return elts; }