Java Code Examples for htsjdk.samtools.CigarOperator#M

/**
 * Determine whether the provided cigar is okay to merge into the reference path
 *
 * @param cigar                the cigar to analyze
 * @param requireFirstElementM if true, require that the first cigar element be an M operator in order for it to be okay
 * @param requireLastElementM  if true, require that the last cigar element be an M operator in order for it to be okay
 * @return true if it's okay to merge, false otherwise
 */
@VisibleForTesting
static boolean cigarIsOkayToMerge(final Cigar cigar, final boolean requireFirstElementM, final boolean requireLastElementM) {
    final List<CigarElement> elements = cigar.getCigarElements();
    final int numElements = elements.size();

    // don't allow more than a couple of different ops
    if (numElements == 0 || numElements > MAX_CIGAR_COMPLEXITY) {
        return false;
    }

    // the first element must be an M
    if (requireFirstElementM && elements.get(0).getOperator() != CigarOperator.M) {
        return false;
    }

    // the last element must be an M
    if (requireLastElementM && elements.get(numElements - 1).getOperator() != CigarOperator.M) {
        return false;
    }

    // note that there are checks for too many mismatches in the dangling branch later in the process

    return true;
}
 

/**
 * Determine whether the provided cigar is okay to merge into the reference path
 *
 * @param cigar    the cigar to analyze
 * @param requireFirstElementM if true, require that the first cigar element be an M operator in order for it to be okay
 * @param requireLastElementM  if true, require that the last cigar element be an M operator in order for it to be okay
 * @return true if it's okay to merge, false otherwise
 */
@VisibleForTesting
static boolean cigarIsOkayToMerge(final Cigar cigar, final boolean requireFirstElementM, final boolean requireLastElementM) {
    final List<CigarElement> elements = cigar.getCigarElements();
    final int numElements = elements.size();

    // don't allow more than a couple of different ops
    if ( numElements == 0 || numElements > MAX_CIGAR_COMPLEXITY ) {
        return false;
    }

    // the first element must be an M
    if ( requireFirstElementM && elements.get(0).getOperator() != CigarOperator.M ) {
        return false;
    }

    // the last element must be an M
    if ( requireLastElementM && elements.get(numElements - 1).getOperator() != CigarOperator.M ) {
        return false;
    }

    // note that there are checks for too many mismatches in the dangling branch later in the process

    return true;
}
 

/**
 * Rejects reads if the sum of the cigar string bases is less than M_BASES_IN_CIGAR, reject the read.
 * Don't process if M_BASES_IN_CIGAR is -1.
 * @return return false if the sum of the matching bases in the cigar is greater than the threshold.
 */
boolean rejectOnCigar(final SAMRecord r) {
	if (this.SUM_MATCHING_BASES==null) return (false);
	Cigar c = r.getCigar();
	int count=0;
	for (CigarElement ce: c.getCigarElements())
		if (ce.getOperator()==CigarOperator.M)
			count+=ce.getLength();
	if (count>=this.SUM_MATCHING_BASES)
		return false;
	return true;
}
 

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

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

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

private SplitPos getSplitPosition(GATKRead read) {
    if (read.getCigar().getFirstCigarElement().getOperator() == CigarOperator.M) {
        final int matchLength = read.getCigar().getCigarElements().stream().filter(e -> e.getOperator().consumesReferenceBases()).mapToInt(CigarElement::getLength).sum();
        return new SplitPos(read.getStart() + matchLength, POSITION.RIGHT);
    } else if (read.getCigar().getLastCigarElement().getOperator() == CigarOperator.M) {
        return new SplitPos(read.getStart(), POSITION.LEFT);
    }

    return new SplitPos(-1, POSITION.MIDDLE);
}
 

private static CigarElement makeElement(final State state, final int length) {
    CigarOperator op = null;
    switch (state) {
        case MATCH: op = CigarOperator.M; break;
        case INSERTION: op = CigarOperator.I; break;
        case DELETION: op = CigarOperator.D; break;
        case CLIP: op = CigarOperator.S; break;
    }
    return new CigarElement(length, op);
}
 

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

public boolean isVisible() {
	return cigarElement.getOperator() == CigarOperator.M || cigarElement.getOperator() == CigarOperator.X || cigarElement.getOperator() == CigarOperator.EQ;
}
 

/**
 * For a read on this pileup, get the base and quality of the base that is at the same
 * position as the SNP.  If the read does not overlap the interval, then return an empty array.
 * @param r The read
 * @return A length 2 byte array containing the base and quality.  Empty if the read does not overlap the interval.
 */
public byte [] getBaseAndQualityOverlappingInterval (final SAMRecord r) {
	byte [] result = new byte [2];

	List<CigarElement> elements = r.getCigar().getCigarElements();
	Iterator<AlignmentBlock> blocks = r.getAlignmentBlocks().iterator();

	int finalSNPLocalPosition=-1;
	int lengthTraversed=0;

	for (CigarElement ce: elements) {
		CigarOperator co = ce.getOperator();
		// you're in an alignment block
		if (co==CigarOperator.M) {
			// get the next alignment block
			AlignmentBlock b = blocks.next();
			int refStart = b.getReferenceStart();
			int snpLocalPos=this.getPosition() - refStart +1;
			int blockLength=b.getLength();

			// is the local position inside this alignment block?
			// if not, move onto the next block.
			if (snpLocalPos >0 && snpLocalPos<=blockLength) {
				// found it!  Done.
				finalSNPLocalPosition=snpLocalPos+lengthTraversed;
				break;
			}
		}
		// consume the bases if necessary and move on to the next element.
		if (co.consumesReadBases())
			lengthTraversed+=ce.getLength();
	}

	// if the position is assigned, then add to the pileup.
	if (finalSNPLocalPosition!=-1) {
		// arrays 0 based.
		byte [] quals = r.getBaseQualities();
		byte [] bases = r.getReadBases();
		byte base = bases[finalSNPLocalPosition-1];
		// char baseChar = StringUtil.byteToChar(base);
		byte qual = quals[finalSNPLocalPosition-1];
		result[0]=base;
		result[1]=qual;
		return (result);
	}
	return (ArrayUtils.EMPTY_BYTE_ARRAY);
}
 

@Test
public void  testForIdenticalAlignmentsWithDifferingFlankLengths() {
    //This test is designed to ensure that the indels are correctly placed
    //if the region flanking these indels is extended by a varying amount.
    //It checks for problems caused by floating point rounding leading to different
    //paths being selected.

    //Create two versions of the same sequence with different flanking regions.
    final byte[] paddedRef = "GCGTCGCAGTCTTAAGGCCCCGCCTTTTCAGACAGCTTCCGCTGGGCCTGGGCCGCTGCGGGGCGGTCACGGCCCCTTTAAGCCTGAGCCCCGCCCCCTGGCTCCCCGCCCCCTCTTCTCCCCTCCCCCAAGCCAGCACCTGGTGCCCCGGCGGGTCGTGCGGCGCGGCGCTCCGCGGTGAGCGCCTGACCCCGAGGGGGCCCGGGGCCGCGTCCCTGGGCCCTCCCCACCCTTGCGGTGGCCTCGCGGGTCCCAGGGGCGGGGCTGGAGCGGCAGCAGGGCCGGGGAGATGGGCGGTGGGGAGCGCGGGAGGGACCGGGCCGAGCCGGGGGAAGGGCTCCGGTGACT"
            .getBytes();
    final byte[] paddedHap = "GCGTCGCAGTCTTAAGGCCCCGCCTTTTCAGACAGCTTCCGCTGGGCCTGGGCCGCTGCGGGGCGGTCACGGCCCCTTTAAGCCTGAGCCCCGCCCCCTGGCTCCCCGCCCCCTCTTCTCCCCTCCCCCAAGCCAGCACCTGGTGCCCCGGCGGGTCGTGCGGCGCGGCGCTCCGCGGTGAGCGCCTGACCCCGA--GGGCC---------------GGGCCCTCCCCACCCTTGCGGTGGCCTCGCGGGTCCCAGGGGCGGGGCTGGAGCGGCAGCAGGGCCGGGGAGATGGGCGGTGGGGAGCGCGGGAGGGACCGGGCCGAGCCGGGGGAAGGGCTCCGGTGACT"
            .replace("-", "")
            .getBytes();
    final byte[] notPaddedRef = "CTTTAAGCCTGAGCCCCGCCCCCTGGCTCCCCGCCCCCTCTTCTCCCCTCCCCCAAGCCAGCACCTGGTGCCCCGGCGGGTCGTGCGGCGCGGCGCTCCGCGGTGAGCGCCTGACCCCGAGGGGGCCCGGGGCCGCGTCCCTGGGCCCTCCCCACCCTTGCGGTGGCCTCGCGGGTCCCAGGGGCGGGGCTGGAGCGGCAGCAGGGCCGGGGAGATGGGCGGTGGGGAGCGCGGGAGGGA"
            .getBytes();
    final byte[] notPaddedHap = "CTTTAAGCCTGAGCCCCGCCCCCTGGCTCCCCGCCCCCTCTTCTCCCCTCCCCCAAGCCAGCACCTGGTGCCCCGGCGGGTCGTGCGGCGCGGCGCTCCGCGGTGAGCGCCTGACCCCGA---------GGGCC--------GGGCCCTCCCCACCCTTGCGGTGGCCTCGCGGGTCCCAGGGGCGGGGCTGGAGCGGCAGCAGGGCCGGGGAGATGGGCGGTGGGGAGCGCGGGAGGGA"
            .replace("-", "")
            .getBytes();
    //a simplified version of the getCigar routine in the haplotype caller to align these
    final String SW_PAD = "NNNNNNNNNN";
    final String paddedsRef = SW_PAD + new String(paddedRef) + SW_PAD;
    final String paddedsHap = SW_PAD + new String(paddedHap) + SW_PAD;
    final String notPaddedsRef = SW_PAD + new String(notPaddedRef) + SW_PAD;
    final String notpaddedsHap = SW_PAD + new String(notPaddedHap) + SW_PAD;
    final SmithWatermanAlignment paddedAlignment;
    final SmithWatermanAlignment notPaddedAlignment;
    try (final SmithWatermanAligner aligner = getAligner()) {
        paddedAlignment = aligner.align(paddedsRef.getBytes(), paddedsHap.getBytes(), CigarUtils.NEW_SW_PARAMETERS, SWOverhangStrategy.SOFTCLIP);
        notPaddedAlignment = aligner.align(notPaddedsRef.getBytes(), notpaddedsHap.getBytes(), CigarUtils.NEW_SW_PARAMETERS, SWOverhangStrategy.SOFTCLIP);
    }
    //Now verify that the two sequences have the same alignment and not match positions.
    final Cigar rawPadded = paddedAlignment.getCigar();
    final Cigar notPadded = notPaddedAlignment.getCigar();
    final List<CigarElement> paddedC = rawPadded.getCigarElements();
    final List<CigarElement> notPaddedC = notPadded.getCigarElements();
    Assert.assertEquals(paddedC.size(), notPaddedC.size());
    for (int i = 0; i < notPaddedC.size(); i++) {
        final CigarElement pc = paddedC.get(i);
        final CigarElement npc = notPaddedC.get(i);
        if (pc.getOperator() == CigarOperator.M && npc.getOperator() == CigarOperator.M) {
            continue;
        }
        final int l1 = pc.getLength();
        final int l2 = npc.getLength();
        Assert.assertEquals(l1, l2);
        Assert.assertEquals(pc.getOperator(), npc.getOperator());
    }
}
 

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

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 List<Feature> getExtraJunctions(ContigAlignerResult result, List<Feature> junctions, List<Feature> junctions2) {
	// Look for junctions with adjacent deletions.
	// Treat these as putative novel junctions.
	Set<Feature> junctionSet = new HashSet<Feature>(junctions);
	junctionSet.addAll(junctions2);
	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;
	
	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;
			} else {
				gapLength += elem.getLength();
				numElems += 1;
			}
		} else {
			if (isInGap) {
				
				if (numElems > 1) {
					long start = result.getGenomicPos() + gapStart;
					long end = start + gapLength;
					Feature junc = new Feature(result.getChromosome(), start, end-1);
					if (!junctionSet.contains(junc)) {
						Logger.info("Extra junction idenfified: %s", junc);
						extraJunctions.add(junc);
					}
				}
				
				isInGap = false;
				gapStart = -1;
				gapLength = -1;
				numElems = -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();
		}
	}
	
	return extraJunctions;
}
 

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

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

public static final List<int[]> generateMappedCoords(final Cigar cigar, int posStart)
{
    final List<int[]> mappedCoords = Lists.newArrayList();

    // first establish whether the read is split across 2 distant regions, and if so which it maps to
    int posOffset = 0;
    boolean continueRegion = false;

    for(CigarElement element : cigar.getCigarElements())
    {
        if(element.getOperator() == CigarOperator.S)
        {
            // nothing to skip
        }
        else if(element.getOperator() == D)
        {
            posOffset += element.getLength();
            continueRegion = true;
        }
        else if(element.getOperator() == CigarOperator.I)
        {
            // nothing to skip
            continueRegion = true;
        }
        else if(element.getOperator() == CigarOperator.N)
        {
            posOffset += element.getLength();
            continueRegion = false;
        }
        else if(element.getOperator() == CigarOperator.M)
        {
            int readStartPos = posStart + posOffset;
            int readEndPos = readStartPos + element.getLength() - 1;

            if(continueRegion && !mappedCoords.isEmpty())
            {
                int[] lastRegion = mappedCoords.get(mappedCoords.size() - 1);
                lastRegion[SE_END] = readEndPos;
            }
            else
            {
                mappedCoords.add(new int[] { readStartPos, readEndPos });
            }

            posOffset += element.getLength();
            continueRegion = false;
        }
    }

    return mappedCoords;
}