public static int calcFragmentLength(final ReadRecord read1, final ReadRecord read2)
{
    int insertSize = abs(read1.fragmentInsertSize());

    if(!read1.containsSplit() && !read2.containsSplit())
        return insertSize;

    // find unique split lengths and remove them

    List<Integer> splitLengths = read1.Cigar.getCigarElements().stream()
            .filter(x -> x.getOperator() == CigarOperator.N).map(x -> x.getLength()).collect(Collectors.toList());

    for(final CigarElement element : read2.Cigar.getCigarElements())
    {
        if(element.getOperator() == CigarOperator.N && !splitLengths.contains(element.getLength()))
            splitLengths.add(element.getLength());
    }

    int totalSplitLength = splitLengths.stream().mapToInt(x -> x).sum();

    return insertSize - totalSplitLength;
}
 

public static Cigar createCigar(int preSc, int match, int postSc)
{
    if (preSc == 0 && match == 0 && postSc == 0)
        return null;

    Cigar cigar = new Cigar();

    if (preSc > 0)
        cigar.add(new CigarElement(preSc, CigarOperator.SOFT_CLIP));

    if (match > 0)
        cigar.add(new CigarElement(match, CigarOperator.MATCH_OR_MISMATCH));

    if (postSc > 0)
        cigar.add(new CigarElement(postSc, CigarOperator.SOFT_CLIP));

    return cigar;
}
 

public static Cigar createCigar(int preSc, int preMatch, int nonRef, int postMatch, int postSc)
{
    if (preSc == 0 && preMatch == 0 && nonRef == 0 && postMatch == 0 && postSc == 0)
        return null;

    Cigar cigar = new Cigar();

    if (preSc > 0)
        cigar.add(new CigarElement(preSc, CigarOperator.SOFT_CLIP));

    if (preMatch > 0)
        cigar.add(new CigarElement(preMatch, CigarOperator.MATCH_OR_MISMATCH));

    if (nonRef > 0)
        cigar.add(new CigarElement(nonRef, CigarOperator.SKIPPED_REGION));

    if (postMatch > 0)
        cigar.add(new CigarElement(postMatch, CigarOperator.MATCH_OR_MISMATCH));

    if (postSc > 0)
        cigar.add(new CigarElement(postSc, CigarOperator.SOFT_CLIP));

    return cigar;
}
 

public static int numberOfEvents(@NotNull final SAMRecord record) {
    Object nm = record.getAttribute("NM");
    if (!(nm instanceof Integer)) {
        return 0;
    }

    int additionalIndels = 0;
    int softClips = 0;
    for (CigarElement cigarElement : record.getCigar()) {
        switch (cigarElement.getOperator()) {
            case D:
            case I:
                additionalIndels += (cigarElement.getLength() - 1);
                break;
            case S:
                softClips++;
                break;
        }
    }

    return (Integer) nm - additionalIndels + softClips;
}
 

@Nullable
private AltRead processInsert(@NotNull final CigarElement e, @NotNull final SAMRecord record, int readIndex, int refPosition,
        final IndexedBases refBases, int numberOfEvents) {
    int refIndex = refBases.index(refPosition);

    if (refPosition <= bounds.end() && refPosition >= bounds.start()) {
        final String ref = new String(refBases.bases(), refIndex, 1);
        final String alt = new String(record.getReadBases(), readIndex, e.getLength() + 1);
        boolean findReadContext = findReadContext(readIndex, record);

        final RefContext refContext = candidates.refContext(record.getContig(), refPosition);
        if (!refContext.reachedLimit()) {
            final int baseQuality = baseQuality(readIndex, record, alt.length());
            final ReadContext readContext =
                    findReadContext ? readContextFactory.createInsertContext(alt, refPosition, readIndex, record, refBases) : null;
            return new AltRead(refContext, ref, alt, baseQuality, numberOfEvents, readContext);
        }
    }

    return null;
}
 

@Nullable
private AltRead processDel(@NotNull final CigarElement e, @NotNull final SAMRecord record, int readIndex, int refPosition,
        final IndexedBases refBases, int numberOfEvents) {
    int refIndex = refBases.index(refPosition);

    if (refPosition <= bounds.end() && refPosition >= bounds.start()) {
        final String ref = new String(refBases.bases(), refIndex, e.getLength() + 1);
        final String alt = new String(record.getReadBases(), readIndex, 1);
        boolean findReadContext = findReadContext(readIndex, record);

        final RefContext refContext = candidates.refContext(record.getContig(), refPosition);
        if (!refContext.reachedLimit()) {
            final int baseQuality = baseQuality(readIndex, record, 2);
            final ReadContext readContext =
                    findReadContext ? readContextFactory.createDelContext(ref, refPosition, readIndex, record, refBases) : null;
            return new AltRead(refContext, ref, alt, baseQuality, numberOfEvents, readContext);
        }
    }

    return null;
}
 

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

    long refPositionEnd = refPosition + element.getLength() - 1;
    if (refPositionEnd < variant.position()) {
        throw new IllegalStateException("Variant is after record");
    }

    if (variant.position() >= refPosition && variant.position() <= refPositionEnd) {
        int alignmentEnd = record.getAlignmentEnd();
        int actualIndex = record.getReadPositionAtReferencePosition(alignmentEnd) - 1 - alignmentEnd + (int) variant.position();
        result = RawContext.inSoftClip(actualIndex);
    }
}
 

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

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

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

    if (e.getLength() > maxSkippedReferenceRegions) {
        int refPositionEnd = refPosition + e.getLength();
        if (refPositionEnd >= variant.position()) {
            result = RawContext.inSkipped(readIndex);
        }
    }

    handleDelete(record, e, readIndex, refPosition);
}
 

private void mergeAbuttingElements(List<CigarElement> elems) {
	int origSize = elems.size();
	List<CigarElement> orig = Logger.LEVEL == Level.TRACE ? new ArrayList<CigarElement>(elems) : null;
	
	for (int i=1; i<elems.size(); i++) {
		if (elems.get(i).getOperator() == elems.get(i-1).getOperator()) {
			int length = elems.get(i).getLength() + elems.get(i-1).getLength();
			elems.set(i-1, new CigarElement(length, elems.get(i-1).getOperator()));
			elems.remove(i);
			i = i-1;
		}
	}
	
	if (Logger.LEVEL == Level.TRACE && origSize != elems.size()) {
		Cigar old = new Cigar(orig);
		Cigar curr = new Cigar(elems);
		Logger.trace("Cigar abutting elements merged: [%s] -> [%s]", old, curr);
	}
}
 

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

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

public int getAdjustedAlignmentStart() {
	
	int start = 0;
	
	if (adjustedAlignmentStart > -1) {
		start = adjustedAlignmentStart;
	} else {
	
		start = samRecord.getAlignmentStart();
		
		if (samRecord.getCigar().numCigarElements() > 0) {
			CigarElement elem = samRecord.getCigar().getCigarElement(0);
			if (elem.getOperator() == CigarOperator.S) {
				start -= elem.getLength();
				if (start < 1) {
					start = 1;
				}
			}
		}
	}
	
	return start;
}
 

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

public static String getMappedReadPortion(SAMRecord read) {
	int start = 0; 
	
	List<CigarElement> elems = read.getCigar().getCigarElements();
	int i = 0;
	while (i < elems.size() && isClip(elems.get(i))) {
		if (elems.get(i).getOperator() == CigarOperator.S) {
			start += elems.get(i).getLength();
		}
		
		i += 1;
	}
	
	int stop = read.getReadLength();
	i = elems.size()-1;
	while (i >= 0 && isClip(elems.get(i))) {
		if (elems.get(i).getOperator() == CigarOperator.S) {
			stop -= elems.get(i).getLength();
		}
		
		i -= 1;
	}
	
	return read.getReadString().substring(start, stop);
}
 

public static String getLeadingClips(SAMRecord read) {
	List<CigarElement> elems = read.getCigar().getCigarElements();
	
	List<CigarElement> leading = new ArrayList<CigarElement>();
	
	for (CigarElement elem : elems) {
		if (isClip(elem)) {
			leading.add(elem);
		} else {
			break;
		}
	}
	
	String ret = "";
	if (leading.size() > 0) {
		Cigar cigar = new Cigar(leading);
		ret = TextCigarCodec.encode(cigar);
	}
	
	return ret;
}
 

public static String getTrailingClips(SAMRecord read) {
	List<CigarElement> elems = read.getCigar().getCigarElements();
	
	List<CigarElement> trailing = new ArrayList<CigarElement>();
	boolean isNonClippedReached = false;
	
	for (CigarElement elem : elems) {
		if (isClip(elem)) {
			if (isNonClippedReached) {
				trailing.add(elem);
			}
		} else {
			isNonClippedReached = true;
		}
	}

	String ret = "";
	if (trailing.size() > 0) {
		Cigar cigar = new Cigar(trailing);
		ret =  TextCigarCodec.encode(cigar);
	}
	
	return ret;
}
 

@Test (groups = "unit" )
public void testShiftCigarLeft_insertAtTail() {
	Cigar cigar = new Cigar();
	
	cigar.add(new CigarElement(40, CigarOperator.M));
	cigar.add(new CigarElement(10, CigarOperator.I));
	
	Cigar newCigar;

	newCigar = indelShifter.shiftCigarLeft(cigar, 40);
	Assert.assertEquals(newCigar.toString(), "10I40M");
	
	newCigar = indelShifter.shiftCigarLeft(cigar, 39);
	Assert.assertEquals(newCigar.toString(), "1M10I39M");

	newCigar = indelShifter.shiftCigarLeft(cigar, 30);
	Assert.assertEquals(newCigar.toString(), "10M10I30M");
	
	newCigar = indelShifter.shiftCigarLeft(cigar, 1);
	Assert.assertEquals(newCigar.toString(), "39M10I1M");
}
 

@Test (groups = "unit" )
public void testShiftCigarLeft_multipleIndels() {
	Cigar cigar = new Cigar();
	
	cigar.add(new CigarElement(20, CigarOperator.M));
	cigar.add(new CigarElement(1, CigarOperator.I));
	cigar.add(new CigarElement(5, CigarOperator.M));
	cigar.add(new CigarElement(3, CigarOperator.D));
	cigar.add(new CigarElement(24, CigarOperator.M));
	
	Cigar newCigar;
	
	newCigar = indelShifter.shiftCigarLeft(cigar, 20);
	Assert.assertEquals(newCigar.toString(), "1I5M3D44M");
	
	newCigar = indelShifter.shiftCigarLeft(cigar, 10);
	Assert.assertEquals(newCigar.toString(), "10M1I5M3D34M");

	newCigar = indelShifter.shiftCigarLeft(cigar, 1);
	Assert.assertEquals(newCigar.toString(), "19M1I5M3D25M");
}
 

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

@Test
public void testGetConsistentExons () {
	// need to build a new read that falls into 2 genes to check if it's inconsistent.
	AnnotationUtils u = AnnotationUtils.getInstance();
	OverlapDetector<Gene> geneOverlapDetector = getGeneOverlapDetector(BAM, GTF);
	Map<String, Gene> map = getGeneMap(geneOverlapDetector);
	// expected UTR and coding.
	SAMRecord rec = getReadByName("000000000-AMY9M:1:2104:20987:12124", BAM);
	// change the cigar so the read maps to two genes exons.
	// rec.setCigar(cigar);
	rec.setAlignmentStart(6269430);
	Cigar c = new Cigar ();
	c.add(new CigarElement(20, CigarOperator.M));
	c.add(new CigarElement(9741377, CigarOperator.N));
	c.add(new CigarElement(30, CigarOperator.M));
	rec.setCigar(c);
	Set<Gene> genes = new HashSet<>();
	genes.add(map.get("RPL22"));
	genes.add(map.get("PLEKHM2"));

	// This read spans 2 genes, so it's not consistent under the strict test.
	Set<Gene> actual = u.getConsistentExons(rec, genes, false);
	Set<Gene> expected = new HashSet<>();
	Assert.assertEquals(actual, expected);

	// under the less strict option, we allow the read to span 2 different genes.
	actual = u.getConsistentExons(rec, genes, true);
	expected = genes;
	Assert.assertEquals(actual, expected);

}
 

private static boolean containsInframeIndel(@NotNull final SAMRecord record) {
    for (final CigarElement cigarElement : record.getCigar().getCigarElements()) {
        if (cigarElement.getOperator().isIndel() && cigarElement.getLength() % 3 == 0) {
            return true;
        }
    }

    return false;
}
 

public static int leftSoftClip(@NotNull final SAMRecord record) {
    Cigar cigar = record.getCigar();
    if (cigar.numCigarElements() > 0) {
        CigarElement firstElement = cigar.getCigarElement(0);
        if (firstElement.getOperator() == CigarOperator.S) {
            return firstElement.getLength();
        }
    }

    return 0;
}
 

public static int rightSoftClip(@NotNull final SAMRecord record) {
    Cigar cigar = record.getCigar();
    if (cigar.numCigarElements() > 0) {
        CigarElement lastLement = cigar.getCigarElement(cigar.numCigarElements() - 1);
        if (lastLement.getOperator() == CigarOperator.S) {
            return lastLement.getLength();
        }
    }

    return 0;
}
 

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

@Override
public void handleLeftSoftClip(@NotNull final SAMRecord record, @NotNull final CigarElement element) {
    if (variant.position() < record.getAlignmentStart()) {
        int readIndex = record.getReadPositionAtReferencePosition(record.getAlignmentStart()) - 1 - record.getAlignmentStart()
                + (int) variant.position() - variant.alt().length() + variant.ref().length();
        result = RawContext.inSoftClip(readIndex);
    }
}
 

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

}
 

private int indelLength(@NotNull final SAMRecord record) {
    int result = 0;
    for (CigarElement cigarElement : record.getCigar()) {
        switch (cigarElement.getOperator()) {
            case I:
            case D:
                result += cigarElement.getLength();
        }

    }

    return result;
}