Java Code Examples for htsjdk.samtools.SAMRecord#NO_ALIGNMENT_REFERENCE_INDEX

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

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

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

public static byte[] getRefBasesOfFail(CramContext cramContext, int refId) {
	if (refId == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX)
		return new byte[0];

	SAMSequenceRecord sequenceRecord = cramContext.samFileHeader.getSequence(refId);
	if (sequenceRecord == null)
		throw new RuntimeException("Reference sequence not found in the header: " + refId);
	byte[] ref = cramContext.referenceSource.getReferenceBases(sequenceRecord, true);
	if (ref == null)
		throw new RuntimeException("Failed to fetch reference bases for sequence " + refId + ", name: "
				+ sequenceRecord.getSequenceName());

	return ref;
}
 

private static void setNextMate(final CramCompressionRecord record, final CramCompressionRecord next) {
	record.mateAlignmentStart = next.alignmentStart;
	record.setMateUnmapped(next.isSegmentUnmapped());
	record.setMateNegativeStrand(next.isNegativeStrand());
	record.mateSequenceID = next.sequenceId;
	if (record.mateSequenceID == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX)
		record.mateAlignmentStart = SAMRecord.NO_ALIGNMENT_START;
}
 

private void advance() {
	int candidateIndex = getIndexOfMinAlignment();
	if (candidateIndex < 0) {
		next = null;
	} else {
		next = records[candidateIndex];
		SAMSequenceRecord sequence = header.getSequence(next.getReferenceName());

		next.setHeader(header);

		next.setReferenceIndex(sequence.getSequenceIndex());

		next.setReadName(sources[candidateIndex].id + delim + next.getReadName());

		if (next.getMateReferenceIndex() == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX) {
			next.setMateAlignmentStart(SAMRecord.NO_ALIGNMENT_START);
		} else {
			SAMSequenceRecord mateSequence = header.getSequence(next.getMateReferenceName());
			next.setMateReferenceIndex(mateSequence.getSequenceIndex());
		}

		if (sources[candidateIndex].it.hasNext())
			records[candidateIndex] = sources[candidateIndex].it.next();
		else
			records[candidateIndex] = null;
	}
}
 

private static void setNextMate(CramCompressionRecord r, CramCompressionRecord next) {
	r.mateAlignmentStart = next.alignmentStart;
	r.setMateUnmapped(next.isSegmentUnmapped());
	r.setMateNegativeStrand(next.isNegativeStrand());
	r.mateSequenceID = next.sequenceId;
	if (r.mateSequenceID == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX)
		r.mateAlignmentStart = SAMRecord.NO_ALIGNMENT_START;
}
 

private static Collection<Entry> getMutliRefEntries(CompressionHeader header, Slice slice, long containerOffset,
		int[] landmarks) throws IllegalArgumentException, IllegalAccessException, IOException {
	final DataReaderFactory dataReaderFactory = new DataReaderFactory();
	final Map<Integer, InputStream> inputMap = new HashMap<Integer, InputStream>();
	for (final Integer exId : slice.external.keySet()) {
		inputMap.put(exId, new ByteArrayInputStream(slice.external.get(exId).getRawContent()));
	}

	final RefSeqIdReader reader = new RefSeqIdReader(slice.sequenceId, slice.alignmentStart,
			ValidationStringency.SILENT);
	dataReaderFactory.buildReader(reader,
			new DefaultBitInputStream(new ByteArrayInputStream(slice.coreBlock.getRawContent())), inputMap, header,
			slice.sequenceId);
	reader.APDelta = header.APDelta;

	for (int i = 0; i < slice.nofRecords; i++) {
		final CramCompressionRecord record = new CramCompressionRecord();
		record.sliceIndex = slice.index;
		record.index = i;

		reader.read();

		if (record.sequenceId == slice.sequenceId) {
			record.sequenceId = slice.sequenceId;
		} else {
			if (record.sequenceId == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX)
				record.sequenceName = SAMRecord.NO_ALIGNMENT_REFERENCE_NAME;
		}
	}

	Map<Integer, AlignmentSpan> spans = reader.getReferenceSpans();
	List<Entry> entries = new ArrayList<CramIndex.Entry>(spans.size());
	for (int seqId : spans.keySet()) {
		Entry e = new Entry();
		e.sequenceId = seqId;
		AlignmentSpan span = spans.get(seqId);
		e.alignmentStart = span.getStart();
		e.alignmentSpan = span.getSpan();
		e.sliceSize = slice.size;
		e.sliceIndex = slice.index;

		e.containerStartOffset = containerOffset;
		e.sliceOffset = landmarks[slice.index];

		entries.add(e);
	}

	return entries;
}
 

protected void doRun() throws IOException {
	cramHeader = CramIO.readCramHeader(cramIS);
	FixBAMFileHeader fix = new FixBAMFileHeader(referenceSource);

	reader = newReader();
	reader.reverseNegativeReads = reverse;
	MAIN_LOOP: while (!brokenPipe.get()
			&& (container = ContainerIO.readContainer(cramHeader.getVersion(), cramIS)) != null) {
		if (container.isEOF())
			break;
		DataReaderFactory f = new DataReaderFactory();

		for (Slice s : container.slices) {
			if (s.sequenceId != SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX && s.sequenceId != -2) {
				SAMSequenceRecord sequence = cramHeader.getSamFileHeader().getSequence(s.sequenceId);

				if (sequence == null)
					throw new RuntimeException("Null sequence for id: " + s.sequenceId);

				ref = referenceSource.getReferenceBases(sequence, true);

				if (!s.validateRefMD5(ref)) {
					log.error(String
							.format("Reference sequence MD5 mismatch for slice: seq id %d, start %d, span %d, expected MD5 %s",
									s.sequenceId, s.alignmentStart, s.alignmentSpan,
									String.format("%032x", new BigInteger(1, s.refMD5))));
					throw new RuntimeException("Reference checksum mismatch.");
				}
			} else
				ref = new byte[0];

			Map<Integer, InputStream> inputMap = new HashMap<Integer, InputStream>();
			for (Integer exId : s.external.keySet()) {
				inputMap.put(exId, new ByteArrayInputStream(s.external.get(exId).getRawContent()));
			}

			reader.referenceSequence = ref;
			reader.prevAlStart = s.alignmentStart;
			reader.substitutionMatrix = container.header.substitutionMatrix;
			reader.recordCounter = 0;
			try {
				f.buildReader(reader,
						new DefaultBitInputStream(new ByteArrayInputStream(s.coreBlock.getRawContent())),
						inputMap, container.header, s.sequenceId);
			} catch (IllegalArgumentException e) {
				throw new RuntimeException(e);
			}

			for (int i = 0; i < s.nofRecords; i++) {
				reader.read();
				if (maxRecords > -1) {
					if (maxRecords == 0)
						break MAIN_LOOP;
					maxRecords--;
				}
			}

			containerHasBeenRead();
		}
	}
	if (!brokenPipe.get())
		reader.finish();
}
 

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

}
 

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