Java Code Examples for htsjdk.samtools.CigarElement#getLength()

private int getReadOffset(CigarElement elem) {
	int offset = -1;
	switch(elem.getOperator()) {
		case M:
		case I:
			offset = elem.getLength();
			break;
		case D:
			offset = 0;
			break;
		default:
			throw new IllegalArgumentException("Unexpected Cigar operator: " + elem.getOperator());
	}
	
	return offset;
}
 

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

@Override
public void handleDelete(@NotNull final SAMRecord record, @NotNull final CigarElement e, final int readIndex, final int refPosition) {
    if (result != null) {
        return;
    }

    int refPositionEnd = refPosition + e.getLength();
    if (refPosition == variant.position()) {
        boolean altSupport = isDelete && e.getLength() == variant.ref().length() - 1 && matchesFirstBase(record,
                readIndex,
                variant.ref());
        int baseQuality = altSupport ? baseQuality(readIndex, record, 2) : 0;
        result = RawContext.indel(readIndex, altSupport, baseQuality);
    } else if (refPositionEnd >= variant.position()) {
        result = RawContext.inDelete(readIndex);
    }
}
 

public int getAdjustedAlignmentEnd() {
	int end = -1;
	
	if (adjustedAlignmentEnd > -1) {
		end = adjustedAlignmentEnd;
	} else {
	
		if (samRecord.getReadUnmappedFlag()) {
			end = samRecord.getAlignmentStart() + samRecord.getReadLength();
		} else {
			// Use standard alignment end and pad for soft clipping if necessary
			end = samRecord.getAlignmentEnd();
			
			if (samRecord.getCigar().numCigarElements() > 0) {
				CigarElement elem = samRecord.getCigar().getCigarElement(samRecord.getCigar().numCigarElements()-1);
				if (elem.getOperator() == CigarOperator.S) {
					end += elem.getLength();
				}
			}
		}
	}

	return end;
}
 

@Override
public void handleAlignment(@NotNull final SAMRecord record, @NotNull final CigarElement e, final int readIndex, final int refPosition) {
    if (result != null) {
        return;
    }

    long refPositionEnd = refPosition + e.getLength() - 1;
    if (refPosition <= variant.position() && variant.position() <= refPositionEnd) {
        int readIndexOffset = (int) (variant.position() - refPosition);
        int variantReadIndex = readIndex + readIndexOffset;

        int baseQuality = record.getBaseQualities()[variantReadIndex];
        boolean altSupport = isSNV && refPositionEnd >= variant.end() && matchesString(record, variantReadIndex, variant.alt());
        boolean refSupport = !altSupport && matchesFirstBase(record, variantReadIndex, variant.ref());
        result = RawContext.snv(variantReadIndex, altSupport, refSupport, baseQuality);
    }
}
 

public static int getUnclippedLength(SAMRecord read) {
	int softClipLen = 0;
	for (CigarElement elem : read.getCigar().getCigarElements()) {
		if (elem.getOperator() == CigarOperator.S) {
			softClipLen += elem.getLength();
		}
	}
	
	return read.getReadLength() - softClipLen;
}
 

private boolean hasIndel(SAMRecord read, LocusInfo locus) {
	int readPos = 0;
	int refPosInRead = read.getAlignmentStart();
	int cigarElementIdx = 0;
	
	while (refPosInRead <= locus.posStop && cigarElementIdx < read.getCigar().numCigarElements() && readPos < read.getReadLength()) {
		CigarElement elem = read.getCigar().getCigarElement(cigarElementIdx++);
		
		switch(elem.getOperator()) {
			case H: //NOOP
				break;
			case I:
				if (isWithin(refPosInRead, locus.posStart, locus.posStop)) {
					return true;
				}
				// Fall through to next case
			case S:
				readPos += elem.getLength();
				break;
			case D:
				if (isWithin(refPosInRead, locus.posStart, locus.posStop)) {
					return true;
				}
				// Fall through to next case
			case N:
				refPosInRead += elem.getLength();
				break;
			case M:
				readPos += elem.getLength();
				refPosInRead += elem.getLength();
				break;
			default:
				throw new IllegalArgumentException("Invalid Cigar Operator: " + elem.getOperator() + " for read: " + read.getSAMString());					
		}
	}
	
	return false;
}
 

private List<SpliceJunction> getJunctions(SAMRecord read) {
	List<SpliceJunction> junctions = new ArrayList<SpliceJunction>();
	
	int pos = read.getAlignmentStart();
	
	for (CigarElement elem : read.getCigar()) {
		switch (elem.getOperator()) {
			case D:
			case M:
				pos += elem.getLength();
				break;
			case N:
				junctions.add(new SpliceJunction(read.getReferenceName(), pos, pos+elem.getLength()-1));
				pos += elem.getLength();
				break;
			case S:
			case I:
			case H:
				// NOOP
				break;
			default:
				throw new UnsupportedOperationException("Unsupported Cigar Operator: " + elem.getOperator());
		}
	}
	
	return junctions;
}
 

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

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

@Override
public void handleInsert(@NotNull final SAMRecord record, @NotNull final CigarElement e, final int readIndex, final int refPosition) {
    if (result != null) {
        return;
    }

    if (refPosition == variant.position()) {
        boolean altSupport = isInsert && e.getLength() == variant.alt().length() - 1 && matchesString(record, readIndex, variant.alt());
        int baseQuality = altSupport ? baseQuality(readIndex, record, variant.alt().length()) : 0;
        result = RawContext.indel(readIndex, altSupport, baseQuality);
    }

}
 

public List<Span> getSpanningRegions() {
	
	List<Span> spans = new ArrayList<Span>();
	
	int start = getAdjustedAlignmentStart();
	
	if (samRecord.getReadUnmappedFlag()) {
		spans.add(new Span(start, getAdjustedAlignmentEnd()));
	} else {
		int end = start;
		for (CigarElement elem : samRecord.getCigar().getCigarElements()) {
			switch (elem.getOperator()) {
			case M:
			case S:
			case D:
				end += elem.getLength();
				break;
			case I:
				break;
			case H:
				break;
			case N:
				spans.add(new Span(start, end));
				start = end + elem.getLength();
				end = start;
				break;
			default:
				throw new UnsupportedOperationException("Unhandled cigar operator: " + elem.getOperator() + " in: " + 
						samRecord.getReadName() + " : " + samRecord.getCigarString());
			}
		}
		
		spans.add(new Span(start, end));
	}
	
	return spans;
}
 

@Override
public void handleAlignment(@NotNull final SAMRecord r, @NotNull final CigarElement e, final int startReadIndex, final int refPos) {
    for (int i = 0; i < e.getLength(); i++) {
        int readIndex = startReadIndex + i;
        int position = refPos + i;

        if (position > bounds.end()) {
            return;
        }

        if (position < bounds.start()) {
            continue;
        }

        byte ref = refGenome.base(position);
        byte alt = r.getReadBases()[readIndex];
        byte quality = r.getBaseQualities()[readIndex];
        byte[] trinucleotideContext = refGenome.trinucleotideContext(position);

        if (alt != N && isValid(trinucleotideContext)) {
            final QualityCounterKey key = ImmutableQualityCounterKey.builder()
                    .ref(ref)
                    .alt(alt)
                    .qual(quality)
                    .position(position)
                    .trinucleotideContext(trinucleotideContext)
                    .build();
            qualityMap.computeIfAbsent(key, QualityCounter::new).increment();
        }
    }
}
 

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

public static int getMappedLength(SAMRecord read) {
	int length = read.getReadLength();
	for (CigarElement elem : read.getCigar().getCigarElements()) {
		if (elem.getOperator() == CigarOperator.S) {
			length -= elem.getLength();
		}
	}
	
	return length;
}
 

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

public String getAlternateReference(int refStart, int refEnd, String chromosome, String seq, Cigar cigar) {
	String alt = null;
	
	if (refEnd < getRefLength(chromosome)) {
		
		StringBuffer altBuf = new StringBuffer(seq.length());
	
		int readIdx = 0;
		int refIdx = refStart-1;
		for (CigarElement element : cigar.getCigarElements()) {
			if (element.getOperator() == CigarOperator.M) {
				
				for (int i=0; i<element.getLength(); i++) {

					if (refIdx >= getRefLength(chromosome)) {
						// You're off the edge of the map matey.  Monsters be here!
						// This read has aligned across chromosomes.  Do not proceed.
						return null;
					}
					
					char refBase = getRefBase(refIdx, chromosome);
					
					altBuf.append(refBase);
										
					readIdx++;
					refIdx++;
				}
			} else if (element.getOperator() == CigarOperator.I) {
				altBuf.append(seq.substring(readIdx, readIdx + element.getLength()));
				readIdx += element.getLength();
			} else if (element.getOperator() == CigarOperator.D) {
				refIdx += element.getLength();
			} else if (element.getOperator() == CigarOperator.S) {
				readIdx += element.getLength();
			}
		}
		
		alt = altBuf.toString();
	}
	
	return alt;
}
 

public static void traverseCigar(@NotNull final SAMRecord record, @NotNull final CigarHandler handler) {
    final Cigar cigar = record.getCigar();

    int readIndex = 0;
    int refBase = record.getAlignmentStart();

    for (int i = 0; i < cigar.numCigarElements(); i++) {
        final CigarElement e = cigar.getCigarElement(i);
        switch (e.getOperator()) {
            case H:
                break; // ignore hard clips
            case P:
                break; // ignore pads
            case S:
                if (i == 0) {
                    handler.handleLeftSoftClip(record, e);
                } else if (i == cigar.numCigarElements() - 1) {
                    handler.handleRightSoftClip(record, e, readIndex, refBase);
                }
                readIndex += e.getLength();
                break; // soft clip read bases
            case N:
                handler.handleSkippedReference(record, e, readIndex - 1, refBase - 1);
                refBase += e.getLength();
                break;  // reference skip
            case D:
                handler.handleDelete(record, e, readIndex - 1, refBase - 1);
                refBase += e.getLength();
                break;
            case I:
                // TODO: Handle 1I150M
                int refIndex = refBase - 1 - record.getAlignmentStart();
                if (refIndex >= 0) {
                    handler.handleInsert(record, e, readIndex - 1, refBase - 1);
                }
                readIndex += e.getLength();
                break;
            case M:
            case EQ:
            case X:
                boolean isFollowedByIndel = i < cigar.numCigarElements() - 1 && cigar.getCigarElement(i + 1).getOperator().isIndel();
                final CigarElement element = isFollowedByIndel ? new CigarElement(e.getLength() - 1, e.getOperator()) : e;
                handler.handleAlignment(record, element, readIndex, refBase);
                readIndex += e.getLength();
                refBase += e.getLength();
                break;
            default:
                throw new IllegalStateException("Case statement didn't deal with op: " + e.getOperator() + "in CIGAR: " + cigar);
        }
    }
}
 

private List<Feature> getExonSkippingJunctions(ContigAlignerResult result, List<Feature> junctions) {

		// Handles special case where an exon skipping junction causes large gaps with a tiny (~1)
		// number of bases mapped somewhere in the gap
		Set<Feature> junctionSet = new HashSet<Feature>(junctions);
		List<Feature> extraJunctions = new ArrayList<Feature>();
		
		Cigar cigar = TextCigarCodec.decode(result.getCigar());
		
		boolean isInGap = false;
		int gapStart = -1;
		int gapLength = -1;
		int numElems = -1;
		int refOffset = 0;
		
		int maxBasesInMiddle = 5;
		int middleBases = -1;
		
		CigarOperator prev = CigarOperator.M;
		
		for (CigarElement elem : cigar.getCigarElements()) {
			if (elem.getOperator() == CigarOperator.D || elem.getOperator() == CigarOperator.N) {
				if (!isInGap) {
					isInGap = true;
					gapStart = refOffset;
					gapLength = elem.getLength();
					numElems = 1;
					middleBases = 0;
				} else {
					gapLength += elem.getLength();
					numElems += 1;
				}
			} else {
				if (isInGap) {
					
					if (elem.getLength() + middleBases < maxBasesInMiddle) {
						middleBases += elem.getLength();
					} else {
						
						if (numElems > 1 && middleBases > 0 && (prev == CigarOperator.D || prev == CigarOperator.N)) {
							long start = result.getGenomicPos() + gapStart;
							long end = start + gapLength;
							
							// Find junction start / end points closest to gap start / end point
							long closestStart = junctions.get(0).getStart();
							long closestEnd = junctions.get(0).getEnd();
							for (Feature junction : junctions) {
								if (Math.abs(junction.getStart()-start) < Math.abs(closestStart-start)) {
									closestStart = junction.getStart();
								}

								if (Math.abs(junction.getEnd()-end) < Math.abs(closestEnd-end)) {
									closestEnd = junction.getEnd();
								}
							}
							
							if (closestEnd > closestStart+1) {
								Feature junc = new Feature(result.getChromosome(), closestStart, closestEnd);
								if (!junctionSet.contains(junc)) {
									Logger.info("Potential exon skipping junction idenfified: %s", junc);
									extraJunctions.add(junc);
								}
							}
						}
						
						isInGap = false;
						gapStart = -1;
						gapLength = -1;
						numElems = -1;
						middleBases = -1;
					}
				}
			}
			
			if (elem.getOperator() == CigarOperator.M || 
				elem.getOperator() == CigarOperator.D || 
				elem.getOperator() == CigarOperator.N ||
				elem.getOperator() == CigarOperator.X ||
				elem.getOperator() == CigarOperator.EQ) {
				
				refOffset += elem.getLength();
			}
			
			prev = elem.getOperator();
		}
		
		return extraJunctions;
	}
 

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