Java Code Examples for htsjdk.samtools.CigarOperator#S

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

@Override
public boolean test(final GATKRead read) {
    int alignedLength = 0;
    int softClipBlocks = 0;
    final int minSoftClipBlocks = doNotRequireSoftClipsOnBothEnds ? 1 : 2;
    CigarOperator prevOperator = null;

    for ( final CigarElement element : read.getCigarElements() ) {
        if ( element.getOperator() == CigarOperator.S ) { // count total S blocks
            //Treat consecutive S blocks as a single one
            if(prevOperator != CigarOperator.S){
                softClipBlocks += 1;
            }

        } else if ( element.getOperator().consumesReadBases() ) {   // M, I, X, and EQ (S was already accounted for above)
            alignedLength += element.getLength();
        }
        prevOperator = element.getOperator();
    }

    return (alignedLength >= minAlignedBases || softClipBlocks < minSoftClipBlocks);
}
 

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

/**
 * Make sure that the SW didn't fail in some terrible way, and throw exception if it did
 */
private static boolean isSWFailure(final SmithWatermanAlignment alignment) {
    // check that the alignment starts at the first base, which it should given the padding
    if ( alignment.getAlignmentOffset() > 0 ) {
        return true;
    }

    // check that we aren't getting any S operators (which would be very bad downstream)
    for ( final CigarElement ce : alignment.getCigar().getCigarElements() ) {
        if ( ce.getOperator() == CigarOperator.S )
            return true;
        // soft clips at the end of the alignment are really insertions
    }

    return false;
}
 

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

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

private boolean testMinLeadingTrailingSoftClippedRatio(final GATKRead read) {

        if ( read.getCigarElements().size() < 1 ) {
            return false;
        }

        // Get the index of the last cigar element:
        final int lastCigarOpIndex = read.getCigarElements().size() - 1;

        // Calculate the number of bases here:
        // Note: we don't want to double-count if the read has only 1 cigar operator.
        final int numLeadingTrailingSoftClippedBases =
                (read.getCigarElement(0).getOperator() == CigarOperator.S ? read.getCigarElement(0).getLength() : 0)
                        +
                ((lastCigarOpIndex != 0) && (read.getCigarElement(lastCigarOpIndex).getOperator() == CigarOperator.S)
                        ? read.getCigarElement(lastCigarOpIndex).getLength()
                        : 0);

        // Determine length:
        final int totalLength = read.getCigarElements().stream()
                .mapToInt( CigarElement::getLength )
                .sum();

        // Calculate the ratio:
        final double softClipRatio = ((double)numLeadingTrailingSoftClippedBases / (double)totalLength);

        return softClipRatio > minimumLeadingTrailingSoftClippedRatio;
    }
 

/**
 * Calculate the index of the current pileup position against the reference-aligned read
 * This offset should be representative of the "IGV view" for the read where insertions are collapsed and deletions
 * are padded so that we can easily count the mismatches against the reference
 * @param p the PileupElement containing the offset as an index into the read base sequence
 * @return the new reference-aligned index/offset
 */
@VisibleForTesting
protected int getCigarModifiedOffset (final PileupElement p){
    final GATKRead read = p.getRead();
    int offset = (p.getCurrentCigarElement().getOperator().consumesReferenceBases() || p.getCurrentCigarElement().getOperator() == CigarOperator.S)? p.getOffsetInCurrentCigar() : 0;
    for (int i = 0; i < p.getCurrentCigarOffset(); i++) {
        final CigarElement elem = read.getCigarElement(i);
        if (elem.getOperator().consumesReferenceBases() || elem.getOperator() == CigarOperator.S) {
            offset += elem.getLength();
        }
    }
    return offset;
}
 

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 isLastNonSoftClippedElem(int idx, Cigar cigar) {
	int numElems = cigar.getCigarElements().size();
	
	// Last element, not soft clipped.
	if ((idx == numElems-1) && (cigar.getCigarElement(idx).getOperator() != CigarOperator.S)) {
		return true;
	}
	
	// Second to last element, with last element soft clipped.
	if ((idx == numElems-2) && (cigar.getCigarElement(numElems-1).getOperator() == CigarOperator.S)) {
		return true;
	}

	return false;
}
 

private boolean isFirstNonSoftClippedElem(int idx, Cigar cigar) {
	// First element, not soft clipped
	if ((idx == 0) && (cigar.getCigarElement(0).getOperator() != CigarOperator.S)) {
		return true;
	}
	
	// Second element, with first element soft clipped
	if ((idx == 1) && (cigar.getCigarElement(0).getOperator() == CigarOperator.S)) {
		return true;
	}
	
	return false;
}
 

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

/**
 * Determine the appropriate start and stop offsets in the reads for the bases given the cigar string
 * @param read
 * @return
 */
private final Pair<Integer,Integer> calculateQueryRange(final GATKRead read) {
    int queryStart = -1, queryStop = -1;
    int readI = 0;

    // iterate over the cigar elements to determine the start and stop of the read bases for the BAQ calculation
    for ( CigarElement elt : read.getCigarElements() ) {
        switch (elt.getOperator()) {
            case N:  return null; // cannot handle these
            case H : case P : case D: break; // ignore pads, hard clips, and deletions
            case I : case S: case M: case EQ: case X:
                int prev = readI;
                readI += elt.getLength();
                if ( elt.getOperator() != CigarOperator.S) {
                    if ( queryStart == -1 ) {
                        queryStart = prev;
                    }
                    queryStop = readI;
                }
                // in the else case we aren't including soft clipped bases, so we don't update
                // queryStart or queryStop
                break;
            default: throw new GATKException("BUG: Unexpected CIGAR element " + elt + " in read " + read.getName());
        }
    }

    if ( queryStop == queryStart ) {
        // this read is completely clipped away, and yet is present in the file for some reason
        // usually they are flagged as non-PF, but it's possible to push them through the BAM
        //System.err.printf("WARNING -- read is completely clipped away: " + read.format());
        return null;
    }

    return new MutablePair<>(queryStart, queryStop);
}
 

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

private static boolean isPadding(final CigarElement elt) {
    return elt.getOperator() == CigarOperator.P || elt.getOperator() == CigarOperator.H || elt.getOperator() == CigarOperator.S;
}
 

private final static boolean mightSupportInsertion(final CigarElement cigarElement) {
    final CigarOperator cigarOperator = cigarElement.getOperator();
    return cigarOperator == CigarOperator.I || cigarOperator == CigarOperator.S;
}
 

private static boolean isClipOrPaddingOperator(final CigarOperator op) {
    return op == CigarOperator.S || op == CigarOperator.H || op == CigarOperator.P;
}
 

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