htsjdk.samtools.CigarElement Java Examples

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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