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

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

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

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

/**
 * Replace hard clips with soft clips.
 */
public static void replaceHardClips(SAMRecord read) {
	Cigar cigar = read.getCigar();
	
	if (cigar.getCigarElements().size() > 0) {
		CigarElement firstElement = cigar.getCigarElement(0);
		CigarElement lastElement  = cigar.getCigarElement(cigar.numCigarElements()-1);
		
		if ((firstElement.getOperator() == CigarOperator.H) ||
			(lastElement.getOperator() == CigarOperator.H)) {
			
			Cigar newCigar = new Cigar();
			
			boolean isFirst = true;
			
			for (CigarElement element : cigar.getCigarElements()) {
				if (element.getOperator() != CigarOperator.H) {
					newCigar.add(element);
				} else {
					CigarElement softClip = new CigarElement(element.getLength(), CigarOperator.S);
					newCigar.add(softClip);
					
					if (isFirst) {
						read.setReadString(padBases(element.getLength()) + read.getReadString());
					} else {
						read.setReadString(read.getReadString() + padBases(element.getLength()));							
					}
				}
				
				isFirst = false;
			}
			
			read.setCigar(newCigar);
		}
	}
}
 

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 Character getBaseAtPosition(SAMRecord read, int refPos) {
	int readPos = 0;
	int refPosInRead = read.getAlignmentStart();
	int cigarElementIdx = 0;
	
	while (refPosInRead <= refPos && cigarElementIdx < read.getCigar().numCigarElements() && readPos < read.getReadLength()) {
		CigarElement elem = read.getCigar().getCigarElement(cigarElementIdx++);
		
		switch(elem.getOperator()) {
			case H: //NOOP
				break;
			case S:
			case I:
				readPos += elem.getLength();
				break;
			case D:
			case N:
				refPosInRead += elem.getLength();
				break;
			case M:
				if (refPos < (refPosInRead + elem.getLength())) {
					readPos += refPos - refPosInRead;
					if (readPos < read.getReadLength()) {
						// Found the base.  Return it
						return read.getReadString().charAt(readPos);
					}
				} else {
					readPos += elem.getLength();
					refPosInRead += elem.getLength();
				}
				break;
			default:
				throw new IllegalArgumentException("Invalid Cigar Operator: " + elem.getOperator() + " for read: " + read.getSAMString());					
		}
	}
	
	return null;
}
 

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

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

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

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

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

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

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

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 boolean qcCheck(SAMRecord sr){
Cigar cigar = sr.getCigar();
int rLen = sr.getReadLength();
int effectiveLen = 0;
if(cigar==null) 
    return false;
else{
    for(final CigarElement ce : cigar.getCigarElements()){
	CigarOperator op = ce.getOperator();
	int cigarLen = ce.getLength();
	switch(op)
	    {
	    case M:
		{
		    effectiveLen += cigarLen;
		    break;
		}
	    case I:
		{
		    effectiveLen += cigarLen;
		    break;
		}
	    default: 
		break;
	    }
    }
}
boolean readdebug = false;
if(readdebug){
    HLA.log.appendln(sr.getSAMString());
    HLA.log.appendln("EffectiveLen:\t" + effectiveLen);
    HLA.log.appendln("ReadLen:\t" + rLen);
}
Integer i = sr.getIntegerAttribute("NM");
int nm = 0;
if(i!=null)
    nm = i.intValue();
if(readdebug)
    HLA.log.appendln("NM=\t" + nm);
if(nm < 16){
    if(readdebug)
	HLA.log.appendln("PASSWED QC");
    return true;
}
if(readdebug){
    HLA.log.appendln("FAILED QC");
    HLA.log.appendln(sr.getSAMString());
}
return false;
   }
 

private static boolean isClip(CigarElement elem) {
	return elem.getOperator() == CigarOperator.S || elem.getOperator() == CigarOperator.H;
}
 

public static boolean hasPossibleAdapterReadThrough(SAMRecord read, Map<String, SAMRecordWrapper> firstReads, Map<String, SAMRecordWrapper> secondReads) {
	
	boolean hasReadThrough = false;
	
	// Check for fragment read through in paired end
	if (read.getReadPairedFlag() && !read.getReadUnmappedFlag() && !read.getMateUnmappedFlag() &&
			read.getAlignmentStart() == read.getMateAlignmentStart() &&
			read.getReadNegativeStrandFlag() != read.getMateNegativeStrandFlag()) {
		
		SAMRecordWrapper pair = null;
		
		if (read.getFirstOfPairFlag()) {
			pair = secondReads.get(read.getReadName() + "_" + read.getMateAlignmentStart());
		} else {
			pair = firstReads.get(read.getReadName() + "_" + read.getMateAlignmentStart());
		}
		
		if (pair != null && read.getCigar().getCigarElements().size() > 0 && pair.getSamRecord().getCigar().getCigarElements().size() > 0) {
			
			// Looking for something like:
			//     -------->
			//  <--------
			SAMRecord first = null;
			SAMRecord second = null;
			if (read.getReadNegativeStrandFlag()) {
				first = read;
				second = pair.getSamRecord();
			} else {
				first = pair.getSamRecord();
				second = read;
			}
			
			CigarElement firstElement = first.getCigar().getFirstCigarElement();
			CigarElement lastElement = second.getCigar().getLastCigarElement();
			
			if (firstElement.getOperator() == CigarOperator.S && lastElement.getOperator() == CigarOperator.S) {
				// We likely have read through into adapter here.
				hasReadThrough = true;
			}
		}
	}

	return hasReadThrough;
}