Java Code Examples for htsjdk.samtools.CigarOperator#I

The following examples show how to use htsjdk.samtools.CigarOperator#I . You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.

Example 1

Source File: IndelShifter.java From abra2 with MIT License

5 votes

private int firstIndelOffset(Cigar cigar) {
	int pos = 0;
	
	for (CigarElement elem : cigar.getCigarElements()) {
		if ((elem.getOperator() == CigarOperator.D) || (elem.getOperator() == CigarOperator.I)) {
			return pos;
		} else if (elem.getOperator() != CigarOperator.S) {
			pos += elem.getLength();
		}
	}
	
	throw new IllegalArgumentException("No indel for record: [" + cigar + "]");
}

Example 2

Source File: CigarModifierTest.java From VarDictJava with MIT License

5 votes

@Test(dataProvider = "cigar")
public void getCigarOperatorTest(Object cigarObject) {
    Cigar cigar = (Cigar) cigarObject;
    CigarOperator[] operators = new CigarOperator[] {
            CigarOperator.M,
            CigarOperator.S,
            CigarOperator.I,
            CigarOperator.D,
            CigarOperator.N,
            CigarOperator.H,
    };
    for (int i = 0; i < cigar.numCigarElements(); i++) {
        Assert.assertEquals(CigarParser.getCigarOperator(cigar, i), operators[i]);
    }
}

Example 3

Source File: SAMRecordUtils.java From abra2 with MIT License

5 votes

public static int getNumIndels(SAMRecord read) {
	int numGaps = 0;
	
	for (CigarElement element : read.getCigar().getCigarElements()) {
		if ((element.getOperator() == CigarOperator.D) || (element.getOperator() == CigarOperator.I)) {
			numGaps += 1;
		}
	}
	
	return numGaps;
}

Example 4

Source File: SAMRecordUtils.java From abra2 with MIT License

5 votes

/**
 *  Returns total length of deletions and insertions for the input read. 
 */
public static int getNumIndelBases(SAMRecord read) {
	int numIndelBases = 0;
	
	for (CigarElement element : read.getCigar().getCigarElements()) {
		if ((element.getOperator() == CigarOperator.D) || (element.getOperator() == CigarOperator.I)) {
			numIndelBases += element.getLength();
		}
	}
	
	return numIndelBases;
}

Example 5

Source File: LocusIteratorByStateBaseTest.java From gatk with BSD 3-Clause "New" or "Revised" License

5 votes

private LIBSTest makePermutationTest(final List<CigarElement> elements) {
    CigarElement last = null;
    boolean hasMatch = false;

    // starts with D => bad
    if ( startsWithDeletion(elements) )
        return null;

    // ends with D => bad
    if ( elements.get(elements.size()-1).getOperator() == CigarOperator.D )
        return null;

    // make sure it's valid
    String cigar = "";
    for ( final CigarElement ce : elements ) {
        if ( ce.getOperator() == CigarOperator.N )
            return null; // TODO -- don't support N

        // abort on a bad cigar
        if ( last != null ) {
            if ( ce.getOperator() == last.getOperator() )
                return null;
            if ( isIndel(ce) && isIndel(last) )
                return null;
        }

        cigar += ce.getLength() + ce.getOperator().toString();
        last = ce;
        hasMatch = hasMatch || ce.getOperator() == CigarOperator.M;
    }

    if ( ! hasMatch && elements.size() == 1 &&
            ! (last.getOperator() == CigarOperator.I || last.getOperator() == CigarOperator.S))
        return null;

    return new LIBSTest(cigar);
}

Example 6

Source File: NativeAssembler.java From abra2 with MIT License

5 votes

private int firstIndelLength(Cigar cigar) {
	int len = 0;
	
	for (CigarElement elem : cigar.getCigarElements()) {
		if (elem.getOperator() == CigarOperator.D || elem.getOperator() == CigarOperator.I) {
			len = elem.getLength();
			break;
		}
	}
	
	return len;
}

Example 7

Source File: NativeAssembler.java From abra2 with MIT License

5 votes

private int countGaps(Cigar cigar) {
	int count = 0;
	for (CigarElement elem : cigar.getCigarElements()) {
		if (elem.getOperator() == CigarOperator.D || elem.getOperator() == CigarOperator.I || elem.getOperator() == CigarOperator.N) {
			count += 1;
		}
	}
	
	return count;
}

Example 8

Source File: IndelShifter.java From abra2 with MIT License

5 votes

private boolean containsIndel(Cigar cigar) {
	for (CigarElement elem : cigar.getCigarElements()) {
		if ((elem.getOperator() == CigarOperator.D) || (elem.getOperator() == CigarOperator.I)) {
			return true;
		}
	}
	
	return false;
}

Example 9

Source File: ForwardShiftInsertIterator.java From abra2 with MIT License

5 votes

public int getAlignmentStart() {
	int start = read.getAlignmentStart();
	
	if (read.getCigarLength() > 0 && read.getCigar().getCigarElement(0).getOperator() == CigarOperator.I) {
		start = start - 1;
	}
	
	return start;
}

Example 10

Source File: SimpleAlleleCounter.java From abra2 with MIT License

5 votes

private IndelInfo checkForIndelAtLocus(int alignmentStart, Cigar cigar, int refPos) {
	
	IndelInfo ret = null;
	
	int readIdx = 0;
	int currRefPos = alignmentStart;
	for (CigarElement element : cigar.getCigarElements()) {
		if (element.getOperator() == CigarOperator.M) {
			readIdx += element.getLength();
			currRefPos += element.getLength();
		} else if (element.getOperator() == CigarOperator.I) {
			if (currRefPos == refPos+1) {
				ret = new IndelInfo(element, readIdx);
				break;
			}
			readIdx += element.getLength();
		} else if (element.getOperator() == CigarOperator.D) {
			if (currRefPos == refPos+1) {
				ret = new IndelInfo(element, readIdx);
				break;
			}				
			currRefPos += element.getLength();
		} else if (element.getOperator() == CigarOperator.S) {
			readIdx += element.getLength();
		} else if (element.getOperator() == CigarOperator.N) {
			currRefPos += element.getLength();
		}
		
		if (currRefPos > refPos+1) {
			break;
		}
	}
	
	return ret;
}

Example 11

Source File: CadabraProcessor.java From abra2 with MIT License

5 votes

private IndelInfo checkForIndelAtLocus(int alignmentStart, Cigar cigar, int refPos) {
	
	IndelInfo ret = null;
	
	int readIdx = 0;
	int currRefPos = alignmentStart;
	for (CigarElement element : cigar.getCigarElements()) {
		if (element.getOperator() == CigarOperator.M) {
			readIdx += element.getLength();
			currRefPos += element.getLength();
		} else if (element.getOperator() == CigarOperator.I) {
			if (currRefPos == refPos+1) {
				ret = new IndelInfo(element, readIdx);
				break;
			}
			readIdx += element.getLength();
		} else if (element.getOperator() == CigarOperator.D) {
			if (currRefPos == refPos+1) {
				ret = new IndelInfo(element, readIdx);
				break;
			}				
			currRefPos += element.getLength();
		} else if (element.getOperator() == CigarOperator.S) {
			readIdx += element.getLength();
		} else if (element.getOperator() == CigarOperator.N) {
			currRefPos += element.getLength();
		}
		
		if (currRefPos > refPos+1) {
			break;
		}
	}
	
	return ret;
}

Example 12

Source File: PileupElement.java From gatk with BSD 3-Clause "New" or "Revised" License

5 votes

/**
 * Get the bases for an insertion that immediately follows this alignment state, or null if none exists
 *
 * @see #getLengthOfImmediatelyFollowingIndel() for details on the meaning of immediately.
 *
 * If the immediately following state isn't an insertion, returns null
 *
 * @return actual sequence of inserted bases, or a null if the event is a deletion or if there is no event in the associated read.
 */
public String getBasesOfImmediatelyFollowingInsertion() {
    final CigarElement element = getNextIndelCigarElement();
    if ( element != null && element.getOperator() == CigarOperator.I ) {
        final int getFrom = offset + 1;
        final byte[] bases = Arrays.copyOfRange(read.getBases(), getFrom, getFrom + element.getLength());
        return new String(bases);
    } else {
        return null;
    }
}

Example 13

Source File: ReadLocusReader.java From abra2 with MIT License

5 votes

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;
}

Example 14

Source File: SimpleCaller.java From abra2 with MIT License

4 votes

private BaseInfo getBaseAtReferencePosition(SAMRecord read, int refPos) {
	boolean isNearEdge = false;
	boolean isIndel = false;
	int alignmentStart = read.getAlignmentStart();
	Cigar cigar = read.getCigar();
	
	char base = 'N';
	
	int readIdx = 0;
	int currRefPos = alignmentStart;
	
	for (CigarElement element : cigar.getCigarElements()) {
					
		if (element.getOperator() == CigarOperator.M) {
			readIdx += element.getLength();
			currRefPos += element.getLength();
			
			if (currRefPos > refPos) {  // TODO: double check end of read base here...
				
				int offset = currRefPos - refPos;
				
				if ((offset < 3) || (offset+3 >= element.getLength())) {
					// This position is within 3 bases of start/end of alignment or clipping or indel.
					// Tag this base for evaluation downstream.
					isNearEdge = true;
				}
				
				readIdx -= offset;
				if ((readIdx < 0) || (readIdx >= read.getReadBases().length)) {
					System.err.println("Read index out of bounds for read: " + read.getSAMString());
					break;
				}
				
				if (read.getBaseQualities()[readIdx] >= MIN_BASE_QUALITY) {
					base = (char) read.getReadBases()[readIdx];
				}
				break;
			}
		} else if (element.getOperator() == CigarOperator.I) {
			if (currRefPos == refPos) {
				//TODO: Handle insertions
				isIndel = true;
				break;
			}
			readIdx += element.getLength();
		} else if (element.getOperator() == CigarOperator.D) {
			if (refPos >= currRefPos && refPos <= currRefPos+element.getLength()) {
				//TODO: Handle deletions
				isIndel = true;
				break;
			}				
			currRefPos += element.getLength();
		} else if (element.getOperator() == CigarOperator.S) {
			readIdx += element.getLength();
		}			
	}
	
	return new BaseInfo(Character.toUpperCase(base), isNearEdge, isIndel);
}

Example 15

Source File: PileupElementUnitTest.java From gatk with BSD 3-Clause "New" or "Revised" License

4 votes

@Test(dataProvider = "PileupElementTest")
public void testPileupElementTest(final LIBSTest params) {
    final GATKRead read = params.makeRead();
    final AlignmentStateMachine state = new AlignmentStateMachine(read);
    final LIBS_position tester = new LIBS_position(read);

    while (state.stepForwardOnGenome() != null) {
        tester.stepForwardOnGenome();
        final PileupElement pe = state.makePileupElement();

        Assert.assertEquals(pe.getRead(), read);
        Assert.assertEquals(pe.getMappingQual(), read.getMappingQuality());
        Assert.assertEquals(pe.getOffset(), state.getReadOffset());

        Assert.assertEquals(pe.isDeletion(), state.getCigarOperator() == D);
        Assert.assertEquals(pe.isAfterInsertion(), tester.isAfterInsertion);
        Assert.assertEquals(pe.isBeforeInsertion(), tester.isBeforeInsertion);
        Assert.assertEquals(pe.isNextToSoftClip(), tester.isNextToSoftClip);

        if (!hasNeighboringPaddedOps(params.getElements(), pe.getCurrentCigarOffset())) {
            Assert.assertEquals(pe.isAfterDeletionEnd(), tester.isAfterDeletionEnd);
            Assert.assertEquals(pe.isBeforeDeletionStart(), tester.isBeforeDeletionStart);
        }

        Assert.assertEquals(pe.getOffsetInCurrentCigar(), tester.getCurrentPositionOnOperatorBase0(), "CigarElement index failure");
        Assert.assertTrue(pe.getOffsetInCurrentCigar() >= 0, "Offset into current cigar too small");
        Assert.assertTrue(pe.getOffsetInCurrentCigar() < pe.getCurrentCigarElement().getLength(), "Offset into current cigar too big");

        Assert.assertNotNull(pe.toString());

        Assert.assertEquals(pe.atEndOfCurrentCigar(), state.getOffsetIntoCurrentCigarElement() == state.getCurrentCigarElement().getLength() - 1, "atEndOfCurrentCigar failed");
        Assert.assertEquals(pe.atStartOfCurrentCigar(), state.getOffsetIntoCurrentCigarElement() == 0, "atStartOfCurrentCigar failed");

        Assert.assertEquals(pe.getBase(), pe.isDeletion() ? PileupElement.DELETION_BASE : read.getBases()[state.getReadOffset()]);
        Assert.assertEquals(pe.getQual(), pe.isDeletion() ? PileupElement.DELETION_QUAL : read.getBaseQualities()[state.getReadOffset()]);

        Assert.assertEquals(pe.getCurrentCigarElement(), state.getCurrentCigarElement());
        Assert.assertEquals(pe.getCurrentCigarOffset(), state.getCurrentCigarElementOffset());

        final int lengthOfImmediatelyFollowingIndel = pe.getLengthOfImmediatelyFollowingIndel();
        final String basesOfImmediatelyFollowingInsertion = pe.getBasesOfImmediatelyFollowingInsertion();
        Assert.assertTrue(lengthOfImmediatelyFollowingIndel != 0 || basesOfImmediatelyFollowingInsertion == null);
        Assert.assertTrue(basesOfImmediatelyFollowingInsertion == null || basesOfImmediatelyFollowingInsertion.length() == lengthOfImmediatelyFollowingIndel);

        // Don't test -- pe.getBaseIndex();
        if ( pe.atEndOfCurrentCigar() && state.getCurrentCigarElementOffset() < read.numCigarElements() - 1 ) {
            final CigarElement nextElement = read.getCigar().getCigarElement(state.getCurrentCigarElementOffset() + 1);
            if (nextElement.getOperator() == CigarOperator.I) {
                Assert.assertTrue(pe.getBetweenNextPosition().size() >= 1);
                Assert.assertEquals(pe.getBetweenNextPosition().get(0), nextElement);
            }
            if (nextElement.getOperator() == M) {
                Assert.assertTrue(pe.getBetweenNextPosition().isEmpty());
            }
        } else {
            Assert.assertTrue(pe.getBetweenNextPosition().isEmpty());
        }

        if (pe.atStartOfCurrentCigar() && state.getCurrentCigarElementOffset() > 0) {
            final CigarElement prevElement = read.getCigar().getCigarElement(state.getCurrentCigarElementOffset() - 1);
            if (prevElement.getOperator() == CigarOperator.I) {
                Assert.assertTrue(pe.getBetweenPrevPosition().size() >= 1);
                Assert.assertEquals(pe.getBetweenPrevPosition().get(pe.getBetweenPrevPosition().size() - 1), prevElement);
            }
            if (prevElement.getOperator() == M) {
                Assert.assertTrue(pe.getBetweenPrevPosition().isEmpty());
            }
        } else {
            Assert.assertTrue(pe.getBetweenPrevPosition().isEmpty());
        }

        // TODO -- add meaningful tests
        pe.getBaseInsertionQual();
        pe.getBaseDeletionQual();
    }
}

Example 16

Source File: SimpleAlleleCounter.java From abra2 with MIT License

4 votes

private IndelInfo checkForIndelAtLocus(SAMRecord read, int refPos) {
		IndelInfo elem = null;
		
//		if (refPos == 105243047 && read.getReadName().equals("D7T4KXP1:400:C5F94ACXX:5:2302:20513:30410")) {
//			System.out.println("bar");
//		}
		
		String contigInfo = read.getStringAttribute("YA");
		if (contigInfo != null) {
			// Get assembled contig info.
			String[] fields = contigInfo.split(":");
			int contigPos = Integer.parseInt(fields[1]);
			
			Cigar contigCigar = TextCigarCodec.decode(fields[2]);
			
			// Check to see if contig contains indel at current locus
			elem = checkForIndelAtLocus(contigPos, contigCigar, refPos);
			
			if (elem != null) {
				// Now check to see if this read supports the indel
				IndelInfo readElem = checkForIndelAtLocus(read.getAlignmentStart(),
						read.getCigar(), refPos);
				
				// Allow partially overlapping indels to support contig
				// (Should only matter for inserts)
				if (readElem == null || readElem.getCigarElement().getOperator() != elem.getCigarElement().getOperator()) {
					// Read element doesn't match contig indel
					elem = null;
				} else {
					elem.setReadIndex(readElem.getReadIndex());
					
					// If this read overlaps the entire insert, capture the bases.
					if (elem.getCigarElement().getOperator() == CigarOperator.I) {

						if (elem.getCigarElement().getLength() == readElem.getCigarElement().getLength()) {
					
							String insertBases = read.getReadString().substring(readElem.getReadIndex(), readElem.getReadIndex()+readElem.getCigarElement().getLength());
							elem.setInsertBases(insertBases);
						} else if (readElem.getCigarElement().getLength() < elem.getCigarElement().getLength()) {
							
							int lengthDiff = elem.getCigarElement().getLength() - readElem.getCigarElement().getLength();
							
							if (readElem.getReadIndex() == 0) {
								elem.setReadIndex(readElem.getReadIndex() - lengthDiff);
							} else if (readElem.getReadIndex() == read.getReadLength()-1) {
								elem.setReadIndex(readElem.getReadIndex() + lengthDiff);
							}
						}
					}
				}
			}
		}
		
		return elem;
	}

Example 17

Source File: LocusIteratorByStateBaseTest.java From gatk with BSD 3-Clause "New" or "Revised" License

4 votes

private boolean isIndel(final CigarElement ce) {
    return ce.getOperator() == CigarOperator.D || ce.getOperator() == CigarOperator.I;
}

Example 18

Source File: ContigAligner.java From abra2 with MIT License

4 votes

public ContigAlignerResult align(String seq) {
		
		ContigAlignerResult result = null;
		
		SemiGlobalAligner.Result sgResult = aligner.align(seq, ref);
		Logger.trace("SG Alignment [%s]:\t%s, possible: %d to: %s", seq, sgResult, seq.length()*MATCH, ref);
		if (sgResult.score > MIN_ALIGNMENT_SCORE && sgResult.score > sgResult.secondBest && sgResult.endPosition > 0) {
			Cigar cigar = TextCigarCodec.decode(sgResult.cigar);
			
			CigarElement first = cigar.getFirstCigarElement();
			CigarElement last = cigar.getLastCigarElement();
		
			if (first.getOperator() == CigarOperator.I) {
				Logger.trace("Contig with leading insert: [%s] - [%s]", sgResult, seq);
			}
			
			// Do not allow indels at the edges of contigs.
			if (minAnchorLength > 0 && 
				(first.getOperator() != CigarOperator.M || first.getLength() < minAnchorLength || 
				last.getOperator() != CigarOperator.M || last.getLength() < minAnchorLength)) {

//				if ((first.getOperator() != CigarOperator.M || last.getOperator() != CigarOperator.M) &&
//						cigar.toString().contains("I")) {
//					Logger.trace("INDEL_NEAR_END: %s", cigar.toString());
//					return ContigAlignerResult.INDEL_NEAR_END;
//				}
//				
//				if ((first.getLength() < 5 || last.getLength() < 5) && 
//						cigar.toString().contains("I") &&
//						minAnchorLength >= 5) {
//					Logger.trace("INDEL_NEAR_END: %s", cigar.toString());
//					return ContigAlignerResult.INDEL_NEAR_END;						
//				}

				return null;
			}
				
			int endPos = sgResult.position + cigar.getReferenceLength();
			
			// Require first and last minAnchorLength bases of contig to be similar to ref
			int mismatches = 0;
			for (int i=0; i<minAnchorLength; i++) {
				if (seq.charAt(i) != ref.charAt(sgResult.position+i)) {
					mismatches += 1;
				}
			}
			
			if (mismatches > maxAnchorMismatches) {
				Logger.trace("Mismatches at beginning of: %s", seq);
			} else {
			
				mismatches = 0;
				for (int i=minAnchorLength; i>0; i--) {
					
					int seqIdx = seq.length()-i;
					int refIdx = endPos-i;
					
					if (seq.charAt(seqIdx) != ref.charAt(refIdx)) {
						mismatches += 1;
					}
				}
				
				if (mismatches > maxAnchorMismatches) {
					Logger.trace("Mismatches at end of: %s", seq);
				} else {
					result = finishAlignment(sgResult.position, endPos, sgResult.cigar, sgResult.score, seq);
				}
			}
		}
		
		return result;
	}

Example 19

Source File: LIBS_position.java From gatk with BSD 3-Clause "New" or "Revised" License

4 votes

/**
     * Steps forward on the genome.  Returns false when done reading the read, true otherwise.
     */
    @SuppressWarnings("fallthrough")
    public boolean stepForwardOnGenome() {
        if ( currentOperatorIndex == numOperators )
            return false;

        CigarElement curElement = read.getCigar().getCigarElement(currentOperatorIndex);
        if ( currentPositionOnOperator >= curElement.getLength() ) {
            if ( ++currentOperatorIndex == numOperators )
                return false;

            curElement = read.getCigar().getCigarElement(currentOperatorIndex);
            currentPositionOnOperator = 0;
        }

        switch ( curElement.getOperator() ) {
            case I: // insertion w.r.t. the reference
//                if ( !sawMop )
//                    break;
            case S: // soft clip
                currentReadOffset += curElement.getLength();
            case H: // hard clip
            case P: // padding
                currentOperatorIndex++;
                return stepForwardOnGenome();

            case D: // deletion w.r.t. the reference
            case N: // reference skip (looks and gets processed just like a "deletion", just different logical meaning)
                currentPositionOnOperator++;
                currentGenomeOffset++;
                break;

            case M:
            case EQ:
            case X:
                sawMop = true;
                currentReadOffset++;
                currentPositionOnOperator++;
                currentGenomeOffset++;
                break;
            default:
                throw new IllegalStateException("No support for cigar op: " + curElement.getOperator());
        }

        final boolean isFirstOp = currentOperatorIndex == 0;
        final boolean isLastOp = currentOperatorIndex == numOperators - 1;
        final boolean isFirstBaseOfOp = currentPositionOnOperator == 1;
        final boolean isLastBaseOfOp = currentPositionOnOperator == curElement.getLength();

        isBeforeDeletionStart = isBeforeOp(read.getCigar(), currentOperatorIndex, CigarOperator.D, isLastOp, isLastBaseOfOp);
        isBeforeDeletedBase = isBeforeDeletionStart || (!isLastBaseOfOp && curElement.getOperator() == CigarOperator.D);
        isAfterDeletionEnd = isAfterOp(read.getCigar(), currentOperatorIndex, CigarOperator.D, isFirstOp, isFirstBaseOfOp);
        isAfterDeletedBase  = isAfterDeletionEnd || (!isFirstBaseOfOp && curElement.getOperator() == CigarOperator.D);
        isBeforeInsertion = isBeforeOp(read.getCigar(), currentOperatorIndex, CigarOperator.I, isLastOp, isLastBaseOfOp)
                || (!sawMop && curElement.getOperator() == CigarOperator.I);
        isAfterInsertion = isAfterOp(read.getCigar(), currentOperatorIndex, CigarOperator.I, isFirstOp, isFirstBaseOfOp);
        isNextToSoftClip = isBeforeOp(read.getCigar(), currentOperatorIndex, CigarOperator.S, isLastOp, isLastBaseOfOp)
                || isAfterOp(read.getCigar(), currentOperatorIndex, CigarOperator.S, isFirstOp, isFirstBaseOfOp);

        return true;
    }

Example 20

Source File: RealignmentEngine.java From gatk with BSD 3-Clause "New" or "Revised" License

4 votes

private final static boolean mightSupportInsertion(final CigarElement cigarElement) {
    final CigarOperator cigarOperator = cigarElement.getOperator();
    return cigarOperator == CigarOperator.I || cigarOperator == CigarOperator.S;
}