Java Code Examples for htsjdk.samtools.SAMRecord#getAlignmentEnd()

@NotNull
static Set<VariantHotspot> findInframeIndelsWithIncorrectRefs(@NotNull final SAMRecord record) {
    Set<VariantHotspot> result = Sets.newHashSet();
    if (containsInframeIndel(record)) {
        for (int refPosition = record.getAlignmentStart(); refPosition <= record.getAlignmentEnd(); refPosition++) {
            int readPosition = record.getReadPositionAtReferencePosition(refPosition);

            if (readPosition != 0) {
                int basesInserted = SAMRecords.basesInsertedAfterPosition(record, refPosition);
                if (basesInserted > 0 && basesInserted % 3 == 0) {
                    result.add(createInsert(record, readPosition, refPosition, basesInserted + 1));
                    continue;
                }

                int basesDeleted = SAMRecords.basesDeletedAfterPosition(record, refPosition);
                if (basesDeleted > 0 && basesDeleted % 3 == 0) {
                    result.add(createDelete(record, readPosition, refPosition, basesDeleted + 1));
                }
            }
        }
    }

    return result;
}
 

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

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 void addRead(final GcObject gcObj, final SAMRecord rec, final String group, final byte[] gc, final byte[] refBases) {
    if (!rec.getReadPairedFlag() || rec.getFirstOfPairFlag()) ++gcObj.totalClusters;
    final int pos = rec.getReadNegativeStrandFlag() ? rec.getAlignmentEnd() - scanWindowSize : rec.getAlignmentStart();
    ++gcObj.totalAlignedReads;
    if (pos > 0) {
        final int windowGc = gc[pos];
        if (windowGc >= 0) {
            ++gcObj.readsByGc[windowGc];
            gcObj.basesByGc[windowGc] += rec.getReadLength();
            gcObj.errorsByGc[windowGc] +=
                    SequenceUtil.countMismatches(rec, refBases, bisulfite) +
                            SequenceUtil.countInsertedBases(rec) + SequenceUtil.countDeletedBases(rec);
        }
    }
    if (gcObj.group == null) {
        gcObj.group = group;
    }
}
 

public Map<Gene, LocusFunction> getLocusFunctionForReadByGene (final SAMRecord rec, final OverlapDetector<Gene> geneOverlapDetector) {
	Map<Gene, LocusFunction> result = new HashMap<>();
	final Interval readInterval = new Interval(rec.getReferenceName(), rec.getAlignmentStart(), rec.getAlignmentEnd(), rec.getReadNegativeStrandFlag(), rec.getReadName());
	final Collection<Gene> overlappingGenes = geneOverlapDetector.getOverlaps(readInterval);

       for (Gene g: overlappingGenes) {
       	LocusFunction f = getLocusFunctionForRead(rec, g);
           result.put(g, f);
       }
       return result;
}
 

public static int basesInsertedAfterPosition(@NotNull final SAMRecord record, int position) {
    int startReadPosition = record.getReadPositionAtReferencePosition(position);
    assert (startReadPosition != 0);
    int nextReadPosition = record.getReadPositionAtReferencePosition(position + 1);

    return nextReadPosition == 0 && record.getAlignmentEnd() == position
            ? record.getReadString().length() - startReadPosition
            : Math.max(0, nextReadPosition - startReadPosition - 1);
}
 

public static int basesDeletedAfterPosition(@NotNull final SAMRecord record, int position) {
    int startReadPosition = record.getReadPositionAtReferencePosition(position);
    assert (startReadPosition != 0);

    int nextReferencePosition = record.getReferencePositionAtReadPosition(startReadPosition + 1);
    return nextReferencePosition == 0 && startReadPosition == record.getReadLength()
            ? record.getAlignmentEnd() - position
            : Math.max(0, nextReferencePosition - position - 1);
}
 

static boolean indel(int bafPosition, int readPosition, @NotNull final SAMRecord samRecord) {
    if (samRecord.getAlignmentEnd() > bafPosition) {

        // Delete?
        if (samRecord.getReadPositionAtReferencePosition(bafPosition + 1) == 0) {
            return true;
        }

        // Insert?
        return samRecord.getReferencePositionAtReadPosition(readPosition + 1) != bafPosition + 1;
    }

    return false;
}
 

static int getBaseQuality(@NotNull final GenomePosition position, @NotNull final SAMRecord samRecord) {
    // Get quality of base after del if necessary
    for (int pos = (int) position.position(); pos <= samRecord.getAlignmentEnd(); pos++) {
        int readPosition = samRecord.getReadPositionAtReferencePosition(pos);
        if (readPosition != 0) {
            return SAMRecords.getBaseQuality(samRecord, readPosition);
        }
    }

    return 0;
}
 

public static Map<Integer, AlignmentSpan> getSpans(List<SAMRecord> samRecords) {
	Map<Integer, AlignmentSpan> spans = new HashMap<Integer, AlignmentSpan>();
	int unmapped = 0;
	for (final SAMRecord r : samRecords) {

		int refId = r.getReferenceIndex();
		if (refId == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX) {
			unmapped++;
			continue;
		}

		int start = r.getAlignmentStart();
		int end = r.getAlignmentEnd();

		if (spans.containsKey(refId)) {
			spans.get(refId).add(start, end - start, 1);
		} else {
			spans.put(refId, new AlignmentSpan(start, end - start));
		}
	}
	if (unmapped > 0) {
		AlignmentSpan span = new AlignmentSpan(AlignmentSpan.UNMAPPED_SPAN.getStart(),
				AlignmentSpan.UNMAPPED_SPAN.getSpan(), unmapped);
		spans.put(SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX, span);
	}
	return spans;
}
 

private boolean reachedDepthLimit(@NotNull final SAMRecord record) {
    int alignmentStart = record.getAlignmentStart();
    int alignmentEnd = record.getAlignmentEnd();

    RefContext startRefContext = candidates.refContext(bounds.chromosome(), alignmentStart);
    RefContext endRefContext = candidates.refContext(bounds.chromosome(), alignmentEnd);

    return startRefContext.reachedLimit() && endRefContext.reachedLimit();
}
 

private void setMateInfo(SAMRecord read, Map<MateKey, SAMRecord> mates) {
	
	if (read.getReadPairedFlag()) {
		SAMRecord mate = mates.get(getMateKey(read));
		if (mate != null) {
			// Only update mate info if a read has been modified
			if (read.getAttribute("YO") != null || mate.getAttribute("YO") != null) {
				read.setMateAlignmentStart(mate.getAlignmentStart());
				read.setMateUnmappedFlag(mate.getReadUnmappedFlag());
				read.setMateNegativeStrandFlag(mate.getReadNegativeStrandFlag());
				 
				int start = read.getAlignmentStart() < mate.getAlignmentStart() ? read.getAlignmentStart() : mate.getAlignmentStart();
				int stop  = read.getAlignmentEnd() > mate.getAlignmentEnd() ? read.getAlignmentEnd() : mate.getAlignmentEnd();
				
				int insert = stop-start+1;
				
				if (read.getAlignmentStart() > mate.getAlignmentStart()) {
					insert *= -1;
				} else if (read.getAlignmentStart() == mate.getAlignmentStart() && mate.getFirstOfPairFlag()) {
					insert *= -1;
				}
				
				read.setInferredInsertSize(insert);
				
				if (read.getStringAttribute("MC") != null) {
					read.setAttribute("MC", mate.getCigarString());
				}
				
				if (!mate.getReadUnmappedFlag()) {
					read.setMateReferenceName(mate.getReferenceName());
				}
			}
		}
	}
}
 

private boolean getCachedReadsAtCurrentLocus(List<SAMRecord> reads) {
	
	reads.clear();
	
	Iterator<SAMRecord> cacheIter = readCache.iterator();
	
	String nextChr = null;
	int nextPos = -1;
	
	while (cacheIter.hasNext()) {
		SAMRecord read = cacheIter.next();
		
		if (read.getAlignmentEnd() < currentPos && read.getReferenceName().equals(currentChr)) {
			// We've gone past the end of this read, so remove from cache.
			cacheIter.remove();
		} else if (getAlignmentStart(read) <= currentPos && read.getAlignmentEnd() >= currentPos) {
			// This read spans the current locus of interest.
			reads.add(read);
		} else {
			// This read is beyond the current locus.
			if (nextChr == null) {
				nextChr = read.getReferenceName();
				nextPos = getAlignmentStart(read);
			}
		}
	}
	
	if (reads.isEmpty() && nextChr != null) {
		currentChr = nextChr;
		currentPos = nextPos;
		
		return true;
	}
	
	return false;
}
 

/**
 * Find reads in a given range.
 * 
 * <p>
 * TODO add cigar to the list of returned values
 * <p>
 * TODO add pair information to the list of returned values
 * 
 * @param request
 * @return
 * @throws InterruptedException 
 */
@Override
protected void processDataRequest(DataRequest request) throws InterruptedException {
	
	// Read the given region
	CloseableIterator<SAMRecord> iterator = dataSource.query(request.start.chr, request.start.bp.intValue(), request.end.bp.intValue());
	
	// Produce results
	while (iterator.hasNext()) {

		List<Feature> responseList = new LinkedList<Feature>();

		// Split results into chunks
		for (int c = 0; c < RESULT_CHUNK_SIZE && iterator.hasNext(); c++) {
			SAMRecord record = iterator.next();

			// Region for this read
			Region recordRegion = new Region((long) record.getAlignmentStart(), (long) record.getAlignmentEnd(), request.start.chr);

			// Values for this read
			LinkedHashMap<DataType, Object> values = new LinkedHashMap<DataType, Object>();

			Feature read = new Feature(recordRegion, values);

			if (request.getRequestedContents().contains(DataType.ID)) {
				values.put(DataType.ID, record.getReadName());
			}

			if (request.getRequestedContents().contains(DataType.STRAND)) {
				values.put(DataType.STRAND, BamUtils.getStrand(record, coverageType));					
			}

			if (request.getRequestedContents().contains(DataType.QUALITY)) {
				values.put(DataType.QUALITY, record.getBaseQualityString());
			}

			if (request.getRequestedContents().contains(DataType.CIGAR)) {
				Cigar cigar = new Cigar(read, record.getCigar());
				values.put(DataType.CIGAR, cigar);
			}

			// TODO Deal with "=" and "N" in read string
			if (request.getRequestedContents().contains(DataType.SEQUENCE)) {
				String seq = record.getReadString();
				values.put(DataType.SEQUENCE, seq);
			}

			if (request.getRequestedContents().contains(DataType.MATE_POSITION)) {
				
				BpCoord mate = new BpCoord((Long)(long)record.getMateAlignmentStart(),
						new Chromosome(record.getMateReferenceName()));
				
				values.put(DataType.MATE_POSITION, mate);
			}
			
			if (request.getRequestedContents().contains(DataType.BAM_TAG_NH)) {
				Object ng = record.getAttribute("NH");
				if (ng != null) {
					values.put(DataType.BAM_TAG_NH, (Integer)record.getAttribute("NH"));
				}
			}
			
			/*
			 * NOTE! RegionContents created from the same read area has to be equal in methods equals, hash and compareTo. Primary types
			 * should be ok, but objects (including tables) has to be handled in those methods separately. Otherwise tracks keep adding
			 * the same reads to their read sets again and again.
			 */
			responseList.add(read);
		}

		// Send result			
		super.createDataResult(new DataResult(request.getStatus(), responseList));			
	}

	// We are done
	iterator.close();
}
 

SAMRecord advance() {
	while (true) {
		// Pull next record from stream
		if (!wrappedIterator.hasNext())
			return null;

		final SAMRecord record = wrappedIterator.next();
		// If beyond the end of this reference sequence, end iteration
		final int referenceIndex = record.getReferenceIndex();
		if (referenceIndex != mReferenceIndex) {
			if (referenceIndex < 0 || referenceIndex > mReferenceIndex) {
				return null;
			}
			// If before this reference sequence, continue
			continue;
		}
		if (mRegionStart == 0 && mRegionEnd == Integer.MAX_VALUE) {
			// Quick exit to avoid expensive alignment end calculation
			return record;
		}
		final int alignmentStart = record.getAlignmentStart();
		// If read is unmapped but has a coordinate, return it if the
		// coordinate is within
		// the query region, regardless of whether the mapped mate will
		// be returned.
		final int alignmentEnd;
		if (mQueryType == QueryType.STARTING_AT) {
			alignmentEnd = -1;
		} else {
			alignmentEnd = (record.getAlignmentEnd() != SAMRecord.NO_ALIGNMENT_START ? record.getAlignmentEnd()
					: alignmentStart);
		}

		if (alignmentStart > mRegionEnd) {
			// If scanned beyond target region, end iteration
			return null;
		}
		// Filter for overlap with region
		if (mQueryType == QueryType.CONTAINED) {
			if (alignmentStart >= mRegionStart && alignmentEnd <= mRegionEnd) {
				return record;
			}
		} else if (mQueryType == QueryType.OVERLAPPING) {
			if (alignmentEnd >= mRegionStart && alignmentStart <= mRegionEnd) {
				return record;
			}
		} else {
			if (alignmentStart == mRegionStart) {
				return record;
			}
		}
	}
}
 

private void populateNext(boolean force) throws IOException {
  final SAMRecord previousRecord = mNextRecord;
  mNextRecord = null;
  if (force) {
    advanceSubIterator();
  }
  while (mCurrentIt != null) {
    if (!mCurrentIt.hasNext()) {  // Only happens when stream is exhausted, so effectively just closes things out.
      advanceSubIterator();
    } else {
      final SAMRecord rec = nextOrBuffered();
      final int refId = rec.getReferenceIndex();

      if (refId > mCurrentTemplate) { // current offset has exceeded region and block overlapped next template
        recordPreviousAndAdvance(previousRecord, rec);
      } else {

        if (refId < mCurrentTemplate) { // Current block may occasionally return records from the previous template if the block overlaps
          continue;
        }

        final int alignmentStart = rec.getAlignmentStart() - 1; // to 0-based
        int alignmentEnd = rec.getAlignmentEnd(); // to 0-based exclusive = noop
        if (alignmentEnd == 0) {
          // Use the read length to get an end point for unmapped reads
          alignmentEnd = rec.getAlignmentStart() + rec.getReadLength();
        }

        if (alignmentEnd <= mCurrentRegion.getStart()) {  // before region
          continue;
        }

        if (alignmentStart <= mPreviousAlignmentStart) { // this record would have been already returned by an earlier region
          //Diagnostic.developerLog("Ignoring record from earlier block at " + rec.getReferenceName() + ":" + rec.getAlignmentStart());
          ++mDoubleFetched;
          if (mDoubleFetched % 100000 == 0) {
            Diagnostic.developerLog("Many double-fetched records for source " + mLabel + " noticed at " + rec.getReferenceName() + ":" + rec.getAlignmentStart() + " in region " + mCurrentRegion + " (skipping through to " + mPreviousAlignmentStart + ")");
          }
          continue;
        }

        if (alignmentStart >= mCurrentRegion.getEnd()) {  // past current region, advance the iterator and record the furtherest we got
          recordPreviousAndAdvance(previousRecord, rec);
          continue;
        }

        mNextRecord = rec;
        break;
      }
    }
  }
}
 

private void calculateCoverage(DataRequest request, BpCoord from, BpCoord to) throws InterruptedException {	
			
	//query data for full average bins, because merging them later would be difficult
	long start = CoverageTool.getBin(from.bp);		
	long end = CoverageTool.getBin(to.bp) + CoverageTool.BIN_SIZE - 1;

	//if end coordinate is smaller than 1 Picard returns a full chromosome and we'll run out of memory
	if (end < 1) {
		end = 1;
	}
			
	CloseableIterator<SAMRecord> iterator = dataSource.query(from.chr, (int)start, (int)end);		
	
	BaseStorage forwardBaseStorage = new BaseStorage();
	BaseStorage reverseBaseStorage = new BaseStorage();
	
	while (iterator.hasNext()) {
		
		SAMRecord record = iterator.next();
		
		LinkedHashMap<DataType, Object> values = new LinkedHashMap<DataType, Object>();

		Region recordRegion = new Region((long) record.getAlignmentStart(), (long) record.getAlignmentEnd(), request.start.chr);
		
		Feature read = new Feature(recordRegion, values);

		values.put(DataType.ID, record.getReadName());
		
		values.put(DataType.STRAND, BamUtils.getStrand(record, coverageType));
		
		Cigar cigar = new Cigar(read, record.getCigar());
		values.put(DataType.CIGAR, cigar);
		
		String seq = record.getReadString();
		values.put(DataType.SEQUENCE, seq);
		
		
		// Split read into continuous blocks (elements) by using the cigar
		List<ReadPart> parts = Cigar.splitElements(read);
		
		for (ReadPart part : parts) {				 
			
			if (read.values.get(DataType.STRAND) == Strand.FORWARD) {
				forwardBaseStorage.addNucleotideCounts(part);
			} else if (read.values.get(DataType.STRAND) == Strand.REVERSE) {
				reverseBaseStorage.addNucleotideCounts(part);			
			}
		}						
	}
			
	// We are done
	iterator.close();
			
	/* Reads that overlap query regions create nucleotide counts outside the query region.
	 * Remove those extra nucleotide counts, because they don't contain all reads of those regions and would show
	 * wrong information. 
	 */
	Region filterRegion = new Region(start, end, from.chr);
	forwardBaseStorage.filter(filterRegion);
	reverseBaseStorage.filter(filterRegion);

	// Send result		
	LinkedList<Feature> resultList = new LinkedList<Feature>();					
	
	createResultList(from, forwardBaseStorage, resultList, Strand.FORWARD);
	createResultList(from, reverseBaseStorage, resultList, Strand.REVERSE);		
	
	LinkedList<Feature> averageCoverage = CoverageTool.average(resultList, from.chr);
	
	if (request.getRequestedContents().contains(DataType.COVERAGE)) {		
		super.createDataResult(new DataResult(request, resultList));
	}
	
	if (request.getRequestedContents().contains(DataType.COVERAGE_AVERAGE)) {
		super.createDataResult(new DataResult(request, averageCoverage));
	}
}
 

private boolean containsAllMNVPositions(@NotNull final SAMRecord record) {
    final VariantContext lastVariant = region().variants().get(region().variants().size() - 1);
    return record.getAlignmentStart() <= region().start()
            && record.getAlignmentEnd() >= lastVariant.getStart() + lastVariant.getReference().length() - 1;
}
 

public void select(final SAMRecord record, final Consumer<P> handler) {
    long startWithSoftClip = record.getAlignmentStart() - SAMRecords.leftSoftClip(record);
    long endWithSoftClip = record.getAlignmentEnd() + SAMRecords.rightSoftClip(record);

    super.select(startWithSoftClip, endWithSoftClip, handler);
}
 

private boolean samRecordOverlapsVariant(int start, int end, @NotNull final SAMRecord record) {
    return record.getAlignmentStart() <= start && record.getAlignmentEnd() >= end;
}