Java Code Examples for htsjdk.samtools.SAMRecord#getReferenceIndex()
The following examples show how to use
htsjdk.samtools.SAMRecord#getReferenceIndex() .
You can vote up the ones you like or vote down the ones you don't like,
and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example 1
| Source File: UmiUtil.java From picard with MIT License | 6 votes |
|
/**
* Determines if the read represented by a SAM record belongs to the top or bottom strand
* or if it cannot determine strand position due to one of the reads being unmapped.
* Top strand is defined as having a read 1 unclipped 5' coordinate
* less than the read 2 unclipped 5' coordinate. If a read is unmapped
* we do not attempt to determine the strand to which the read or its mate belongs.
* If the mate belongs to a different contig from the read, then the reference
* index for the read and its mate is used in leu of the unclipped 5' coordinate.
* @param rec Record to determine top or bottom strand
* @return Top or bottom strand, unknown if it cannot be determined.
*/
static ReadStrand getStrand(final SAMRecord rec) {
if (rec.getReadUnmappedFlag() || rec.getMateUnmappedFlag()) {
return ReadStrand.UNKNOWN;
}
// If the read pair are aligned to different contigs we use
// the reference index to determine relative 5' coordinate ordering.
// Both the read and its mate should not have their unmapped flag set to true.
if (!rec.getReferenceIndex().equals(rec.getMateReferenceIndex())) {
if (rec.getFirstOfPairFlag() == (rec.getReferenceIndex() < rec.getMateReferenceIndex())) {
return ReadStrand.TOP;
} else {
return ReadStrand.BOTTOM;
}
}
final int read5PrimeStart = (rec.getReadNegativeStrandFlag()) ? rec.getUnclippedEnd() : rec.getUnclippedStart();
final int mate5PrimeStart = (rec.getMateNegativeStrandFlag()) ? SAMUtils.getMateUnclippedEnd(rec) : SAMUtils.getMateUnclippedStart(rec);
if (rec.getFirstOfPairFlag() == (read5PrimeStart <= mate5PrimeStart)) {
return ReadStrand.TOP;
} else {
return ReadStrand.BOTTOM;
}
}
Example 2
| Source File: MultiSamReader.java From abra2 with MIT License | 5 votes |
|
private SAMRecordWrapper getNext(int idx) {
SAMRecordWrapper record = null;
if (iterators[idx].hasNext()) {
SAMRecord read = iterators[idx].next();
// If no genomic location is assigned, we've reached the unmapped read pairs. Do not continue...
// TODO: Need to include these in final bam files
if (read.getReferenceIndex() >= 0) {
record = new SAMRecordWrapper(read, isFiltered(read), shouldAssemble(read), idx);
}
}
return record;
}
Example 3
| Source File: SmartSamWriter.java From rtg-tools with BSD 2-Clause "Simplified" License | 5 votes |
|
@Override
public boolean addRecord(SAMRecord r) throws IOException {
// Don't attempt to reorder all the fully-unmapped (i.e. no reference sequence) records
if (r.getReferenceIndex() == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX) {
if (!mRecordSet.isEmpty()) {
flush();
}
flushRecord(r);
return true;
} else {
return super.addRecord(r);
}
}
Example 4
| Source File: SamCompareUtils.java From rtg-tools with BSD 2-Clause "Simplified" License | 5 votes |
|
/**
* Compares SAM records on reference index, start position and if paired, mated vs unmated.
* @param a first record
* @param b second record
* @return -1 if first record comes before second record, 1 if the converse. 0 for equal on compared values
*/
public static int compareSamRecords(SAMRecord a, SAMRecord b) {
final int thisRef = a.getReferenceIndex() & 0x7FFFFFFF; // makes -1 a bignum
final int thatRef = b.getReferenceIndex() & 0x7FFFFFFF;
final int rc = Integer.compare(thisRef, thatRef);
if (rc != 0) {
return rc;
}
final int ac = Integer.compare(a.getAlignmentStart(), b.getAlignmentStart());
if (ac != 0) {
return ac;
}
// Do this ... (this one doesn't have same ordering as original)
//return Integer.compare(~a.getFlags(), ~b.getFlags());
// or this ...
// Following compares READ_PAIRED_FLAG, PORPER_PAIRED_FLAG, UNMAPPED_FLAG in that order.
// The ^3 toggles values to get the ordering we require
// The reverse makes sure comparison is in the order we require
return Integer.compare(reverse3bits(a.getFlags() ^ 3), reverse3bits(b.getFlags() ^ 3));
// or this ...
// final int rpc = Boolean.compare(b.getReadPairedFlag(), a.getReadPairedFlag());
// if (rpc != 0) {
// return rpc;
// }
//
// if (a.getReadPairedFlag()) {
// assert b.getReadPairedFlag();
// final int mc = Boolean.compare(b.getProperPairFlag(), a.getProperPairFlag());
// if (mc != 0) {
// return mc;
// }
// }
//
// return Boolean.compare(a.getReadUnmappedFlag(), b.getReadUnmappedFlag());
}
Example 5
| Source File: SamMultiRestrictingIterator.java From rtg-tools with BSD 2-Clause "Simplified" License | 5 votes |
|
private void recordPreviousAndAdvance(SAMRecord previousRecord, SAMRecord rec) throws IOException {
if (previousRecord != null && previousRecord.getReferenceIndex() == mCurrentTemplate) {
mPreviousAlignmentStart = previousRecord.getAlignmentStart() - 1; // to 0-based
} else {
mPreviousAlignmentStart = Integer.MIN_VALUE;
}
setBuffered(rec);
advanceSubIterator();
}
Example 6
| Source File: CollectJumpingLibraryMetrics.java From picard with MIT License | 5 votes |
|
/**
* Calculates the mode for outward-facing pairs, using the first SAMPLE_FOR_MODE
* outward-facing pairs found in INPUT
*/
private double getOutieMode() {
int samplePerFile = SAMPLE_FOR_MODE / INPUT.size();
Histogram<Integer> histo = new Histogram<Integer>();
for (File f : INPUT) {
SamReader reader = SamReaderFactory.makeDefault().open(f);
int sampled = 0;
for (Iterator<SAMRecord> it = reader.iterator(); it.hasNext() && sampled < samplePerFile; ) {
SAMRecord sam = it.next();
if (!sam.getFirstOfPairFlag()) {
continue;
}
// If we get here we've hit the end of the aligned reads
if (sam.getReadUnmappedFlag() && sam.getReferenceIndex() == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX) {
break;
} else if (sam.getReadUnmappedFlag() || sam.getMateUnmappedFlag()) {
continue;
} else if ((sam.getAttribute(SAMTag.MQ.name()) == null ||
sam.getIntegerAttribute(SAMTag.MQ.name()) >= MINIMUM_MAPPING_QUALITY) &&
sam.getMappingQuality() >= MINIMUM_MAPPING_QUALITY &&
sam.getMateNegativeStrandFlag() != sam.getReadNegativeStrandFlag() &&
sam.getMateReferenceIndex().equals(sam.getReferenceIndex()) &&
SamPairUtil.getPairOrientation(sam) == PairOrientation.RF) {
histo.increment(Math.abs(sam.getInferredInsertSize()));
sampled++;
}
}
CloserUtil.close(reader);
}
return histo.size() > 0 ? histo.getMode() : 0;
}
Example 7
| Source File: GcBiasMetricsCollector.java From picard with MIT License | 5 votes |
|
@Override
public void acceptRecord(final GcBiasCollectorArgs args) {
final SAMRecord rec = args.getRec();
if (logCounter < 100 && rec.getReadBases().length == 0) {
log.warn("Omitting read " + rec.getReadName() + " with '*' in SEQ field.");
if (++logCounter == 100) {
log.warn("There are more than 100 reads with '*' in SEQ field in file.");
}
return;
}
if (!rec.getReadUnmappedFlag()) {
if (referenceIndex != rec.getReferenceIndex() || gc == null) {
final ReferenceSequence ref = args.getRef();
refBases = ref.getBases();
StringUtil.toUpperCase(refBases);
final int refLength = refBases.length;
final int lastWindowStart = refLength - scanWindowSize;
gc = GcBiasUtils.calculateAllGcs(refBases, lastWindowStart, scanWindowSize);
referenceIndex = rec.getReferenceIndex();
}
addReadToGcData(rec, this.gcData);
if (ignoreDuplicates && !rec.getDuplicateReadFlag()) {
addReadToGcData(rec, this.gcDataNonDups);
}
} else {
updateTotalClusters(rec, this.gcData);
if (ignoreDuplicates && !rec.getDuplicateReadFlag()) {
updateTotalClusters(rec, this.gcDataNonDups);
}
}
}
Example 8
| Source File: CramSerilization.java From cramtools with Apache License 2.0 | 5 votes |
|
private static void updateTracks(final List<SAMRecord> samRecords, final ReferenceTracks tracks) {
for (final SAMRecord samRecord : samRecords) {
if (samRecord.getAlignmentStart() != SAMRecord.NO_ALIGNMENT_START
&& samRecord.getReferenceIndex() == tracks.getSequenceId()) {
int refPos = samRecord.getAlignmentStart();
int readPos = 0;
for (final CigarElement cigarElement : samRecord.getCigar().getCigarElements()) {
if (cigarElement.getOperator().consumesReferenceBases()) {
for (int elementIndex = 0; elementIndex < cigarElement.getLength(); elementIndex++)
tracks.addCoverage(refPos + elementIndex, 1);
}
switch (cigarElement.getOperator()) {
case M:
case X:
case EQ:
for (int pos = readPos; pos < cigarElement.getLength(); pos++) {
final byte readBase = samRecord.getReadBases()[readPos + pos];
final byte refBase = tracks.baseAt(refPos + pos);
if (readBase != refBase)
tracks.addMismatches(refPos + pos, 1);
}
break;
default:
break;
}
readPos += cigarElement.getOperator().consumesReadBases() ? cigarElement.getLength() : 0;
refPos += cigarElement.getOperator().consumesReferenceBases() ? cigarElement.getLength() : 0;
}
}
}
}
Example 9
| Source File: BAMRecordReader.java From Hadoop-BAM with MIT License | 5 votes |
|
/** Note: this is the only getKey function that handles unmapped reads
* specially!
*/
public static long getKey(final SAMRecord rec) {
final int refIdx = rec.getReferenceIndex();
final int start = rec.getAlignmentStart();
if (!(rec.getReadUnmappedFlag() || refIdx < 0 || start < 0))
return getKey(refIdx, start);
// Put unmapped reads at the end, but don't give them all the exact same
// key so that they can be distributed to different reducers.
//
// A random number would probably be best, but to ensure that the same
// record always gets the same key we use a fast hash instead.
//
// We avoid using hashCode(), because it's not guaranteed to have the
// same value across different processes.
int hash = 0;
byte[] var;
if ((var = rec.getVariableBinaryRepresentation()) != null) {
// Undecoded BAM record: just hash its raw data.
hash = (int)MurmurHash3.murmurhash3(var, hash);
} else {
// Decoded BAM record or any SAM record: hash a few representative
// fields together.
hash = (int)MurmurHash3.murmurhash3(rec.getReadName(), hash);
hash = (int)MurmurHash3.murmurhash3(rec.getReadBases(), hash);
hash = (int)MurmurHash3.murmurhash3(rec.getBaseQualities(), hash);
hash = (int)MurmurHash3.murmurhash3(rec.getCigarString(), hash);
}
return getKey0(Integer.MAX_VALUE, hash);
}
Example 10
| Source File: AlignerInstance.java From halvade with GNU General Public License v3.0 | 5 votes |
|
public int writePairedSAMRecordToContext(SAMRecord sam, boolean useCompact) throws IOException, InterruptedException {
int count = 0;
int read1Ref = sam.getReferenceIndex();
int read2Ref = sam.getMateReferenceIndex();
if (!sam.getReadUnmappedFlag() &&
(read1Ref == read2Ref || keepChrSplitPairs) &&
(read1Ref >= 0 || read2Ref >= 0)) {
context.getCounter(HalvadeCounters.OUT_BWA_READS).increment(1);
writableRecord.set(sam);
int beginpos = sam.getAlignmentStart();
HashSet<Integer> keys = null;
if (!mergeBam)
keys = splitter.getRegions(sam, read1Ref, read2Ref);
else {
keys = new HashSet<>();
keys.add(0);
}
for(Integer key : keys) {
if(useCompact) {
writeableCompactRegion.setRegion(key, beginpos);
context.write(writeableCompactRegion, stub);
} else {
writableRegion.setChromosomeRegion(read1Ref, beginpos, key);
context.write(writableRegion, writableRecord);
}
count++;
}
} else {
if(sam.getReadUnmappedFlag())
context.getCounter(HalvadeCounters.OUT_UNMAPPED_READS).increment(1);
else
context.getCounter(HalvadeCounters.OUT_DIFF_CHR_READS).increment(1);
}
return count;
}
Example 11
| Source File: AlignerInstance.java From halvade with GNU General Public License v3.0 | 5 votes |
|
public int writeSAMRecordToContext(SAMRecord sam, boolean useCompact) throws IOException, InterruptedException {
int count = 0;
if (!sam.getReadUnmappedFlag()){
int read1Ref = sam.getReferenceIndex();
context.getCounter(HalvadeCounters.OUT_BWA_READS).increment(1);
writableRecord.set(sam);
int beginpos = sam.getAlignmentStart();
HashSet<Integer> keys = null;
if (!mergeBam)
keys = splitter.getRegions(sam, read1Ref);
else {
keys = new HashSet<>();
keys.add(0);
}
for(Integer key : keys) {
if(useCompact) {
writeableCompactRegion.setRegion(key, beginpos);
context.write(writeableCompactRegion, stub);
} else {
writableRegion.setChromosomeRegion(read1Ref, beginpos, key);
context.write(writableRegion, writableRecord);
}
count++;
}
} else {
context.getCounter(HalvadeCounters.OUT_UNMAPPED_READS).increment(1);
}
return count;
}
Example 12
| Source File: CramSerilization.java From cramtools with Apache License 2.0 | 5 votes |
|
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;
}
Example 13
| Source File: SamRestrictingIterator.java From rtg-tools with BSD 2-Clause "Simplified" License | 4 votes |
|
private void populateNext() {
mNextRecord = null;
while (mIterator.hasNext()) {
final SAMRecord rec = mIterator.next();
final int refId = rec.getReferenceIndex();
if (refId != mTemplate) { // Get current sequence ranges
if (refId > mEndTemplate) { // Past the last template in requested list so we can bail out
//System.out.println("Aborting iterator for record at " + rec.getReferenceName() + ":" + rec.getAlignmentStart() + " as it is after last template");
break;
}
mTemplate = refId;
final RangeList<String> seqList = mRegions.get(mTemplate);
mSeqRegions = (seqList == null || seqList.getRangeList().size() == 0) ? null : seqList.getRangeList();
mCurrentRegionIdx = 0;
mCurrentRegion = mSeqRegions == null ? null : mSeqRegions.get(0);
}
if (mSeqRegions == null) { // No regions on this template
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();
}
// Advance active region if need be
while (mCurrentRegion.getEnd() <= alignmentStart && mCurrentRegionIdx < mSeqRegions.size() - 1) {
++mCurrentRegionIdx;
mCurrentRegion = mSeqRegions.get(mCurrentRegionIdx);
}
if (alignmentEnd <= mCurrentRegion.getStart()) { // before region
continue;
}
if (alignmentStart >= mCurrentRegion.getEnd()) { // past region
if (mTemplate == mEndTemplate) {
//System.out.println("Aborting iterator for record at " + rec.getReferenceName() + ":" + rec.getAlignmentStart() + " as it is after last region");
break;
} else {
continue;
}
}
mNextRecord = rec;
break;
}
}
Example 14
| Source File: Utils.java From cramtools with Apache License 2.0 | 4 votes |
|
/**
* Copied from net.sf.picard.sam.SamPairUtil. This is a more permissive
* version of the method, which does not reset alignment start and reference
* for unmapped reads.
*
* @param rec1
* @param rec2
* @param header
*/
public static void setLooseMateInfo(final SAMRecord rec1, final SAMRecord rec2, final SAMFileHeader header) {
if (rec1.getReferenceName() != SAMRecord.NO_ALIGNMENT_REFERENCE_NAME
&& rec1.getReferenceIndex() == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX)
rec1.setReferenceIndex(header.getSequenceIndex(rec1.getReferenceName()));
if (rec2.getReferenceName() != SAMRecord.NO_ALIGNMENT_REFERENCE_NAME
&& rec2.getReferenceIndex() == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX)
rec2.setReferenceIndex(header.getSequenceIndex(rec2.getReferenceName()));
// If neither read is unmapped just set their mate info
if (!rec1.getReadUnmappedFlag() && !rec2.getReadUnmappedFlag()) {
rec1.setMateReferenceIndex(rec2.getReferenceIndex());
rec1.setMateAlignmentStart(rec2.getAlignmentStart());
rec1.setMateNegativeStrandFlag(rec2.getReadNegativeStrandFlag());
rec1.setMateUnmappedFlag(false);
rec1.setAttribute(SAMTag.MQ.name(), rec2.getMappingQuality());
rec2.setMateReferenceIndex(rec1.getReferenceIndex());
rec2.setMateAlignmentStart(rec1.getAlignmentStart());
rec2.setMateNegativeStrandFlag(rec1.getReadNegativeStrandFlag());
rec2.setMateUnmappedFlag(false);
rec2.setAttribute(SAMTag.MQ.name(), rec1.getMappingQuality());
}
// Else if they're both unmapped set that straight
else if (rec1.getReadUnmappedFlag() && rec2.getReadUnmappedFlag()) {
rec1.setMateNegativeStrandFlag(rec2.getReadNegativeStrandFlag());
rec1.setMateUnmappedFlag(true);
rec1.setAttribute(SAMTag.MQ.name(), null);
rec1.setInferredInsertSize(0);
rec2.setMateNegativeStrandFlag(rec1.getReadNegativeStrandFlag());
rec2.setMateUnmappedFlag(true);
rec2.setAttribute(SAMTag.MQ.name(), null);
rec2.setInferredInsertSize(0);
}
// And if only one is mapped copy it's coordinate information to the
// mate
else {
final SAMRecord mapped = rec1.getReadUnmappedFlag() ? rec2 : rec1;
final SAMRecord unmapped = rec1.getReadUnmappedFlag() ? rec1 : rec2;
mapped.setMateReferenceIndex(unmapped.getReferenceIndex());
mapped.setMateAlignmentStart(unmapped.getAlignmentStart());
mapped.setMateNegativeStrandFlag(unmapped.getReadNegativeStrandFlag());
mapped.setMateUnmappedFlag(true);
mapped.setInferredInsertSize(0);
unmapped.setMateReferenceIndex(mapped.getReferenceIndex());
unmapped.setMateAlignmentStart(mapped.getAlignmentStart());
unmapped.setMateNegativeStrandFlag(mapped.getReadNegativeStrandFlag());
unmapped.setMateUnmappedFlag(false);
unmapped.setInferredInsertSize(0);
}
boolean firstIsFirst = rec1.getAlignmentStart() < rec2.getAlignmentStart();
int insertSize = firstIsFirst ? SamPairUtil.computeInsertSize(rec1, rec2) : SamPairUtil.computeInsertSize(rec2,
rec1);
rec1.setInferredInsertSize(firstIsFirst ? insertSize : -insertSize);
rec2.setInferredInsertSize(firstIsFirst ? -insertSize : insertSize);
}
Example 15
| Source File: BAMQueryFilteringIterator.java From cramtools with Apache License 2.0 | 4 votes |
|
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;
}
}
}
}
Example 16
| Source File: CollectIndependentReplicateMetrics.java From picard with MIT License | 4 votes |
|
private static QueryInterval queryIntervalFromSamRecord(final SAMRecord samRecord) {
return new QueryInterval(samRecord.getReferenceIndex(), samRecord.getStart(), samRecord.getEnd());
}
Example 17
| Source File: ReorderSam.java From picard with MIT License | 4 votes |
|
/**
* Helper function that writes reads from iterator it into writer out, updating each SAMRecord along the way
* according to the newOrder mapping from dictionary index -> index. Name is used for printing only.
*/
private void writeReads(final SAMFileWriter out,
final SAMRecordIterator it,
final Map<Integer, Integer> newOrder,
final String name) {
long counter = 0;
log.info(" Processing " + name);
while (it.hasNext()) {
counter++;
final SAMRecord read = it.next();
final int oldRefIndex = read.getReferenceIndex();
final int oldMateIndex = read.getMateReferenceIndex();
final int newRefIndex = newOrderIndex(read, oldRefIndex, newOrder);
read.setHeader(out.getFileHeader());
read.setReferenceIndex(newRefIndex);
// read becoming unmapped
if (oldRefIndex != NO_ALIGNMENT_REFERENCE_INDEX &&
newRefIndex == NO_ALIGNMENT_REFERENCE_INDEX) {
read.setAlignmentStart(NO_ALIGNMENT_START);
read.setReadUnmappedFlag(true);
read.setCigarString(SAMRecord.NO_ALIGNMENT_CIGAR);
read.setMappingQuality(SAMRecord.NO_MAPPING_QUALITY);
}
final int newMateIndex = newOrderIndex(read, oldMateIndex, newOrder);
if (oldMateIndex != NO_ALIGNMENT_REFERENCE_INDEX &&
newMateIndex == NO_ALIGNMENT_REFERENCE_INDEX) { // mate becoming unmapped
read.setMateAlignmentStart(NO_ALIGNMENT_START);
read.setMateUnmappedFlag(true);
read.setAttribute(SAMTag.MC.name(), null); // Set the Mate Cigar String to null
}
read.setMateReferenceIndex(newMateIndex);
out.addAlignment(read);
}
it.close();
log.info("Wrote " + counter + " reads");
}
Example 18
| Source File: SamMultiRestrictingIterator.java From rtg-tools with BSD 2-Clause "Simplified" License | 4 votes |
|
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;
}
}
}
}
Example 19
| Source File: AlignmentsTags.java From cramtools with Apache License 2.0 | 3 votes |
|
/**
* Convenience method. See
* {@link net.sf.cram.AlignmentsTags#calculateMdAndNmTags(SAMRecord, byte[], int, boolean, boolean)}
* for details.
*
* @param samRecord
* @param referenceSource
* @param samSequenceDictionary
* @param calculateMdTag
* @param calculateNmTag
* @throws IOException
*/
public static void calculateMdAndNmTags(SAMRecord samRecord, ReferenceSource referenceSource,
SAMSequenceDictionary samSequenceDictionary, boolean calculateMdTag, boolean calculateNmTag)
throws IOException {
if (!samRecord.getReadUnmappedFlag() && samRecord.getReferenceIndex() != SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX) {
SAMSequenceRecord sequence = samSequenceDictionary.getSequence(samRecord.getReferenceIndex());
ReferenceRegion region = referenceSource.getRegion(sequence, samRecord.getAlignmentStart(),
samRecord.getAlignmentEnd());
calculateMdAndNmTags(samRecord, region.array, samRecord.getAlignmentStart() - 1, calculateMdTag,
calculateNmTag);
}
}
