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

void writeSequence(SequencesReader reader, long seqId, byte[] dataBuffer, byte[] qualityBuffer, int flags) throws IOException {
  final SAMRecord rec = new SAMRecord(mWriter.getFileHeader());

  final int length = reader.read(seqId, dataBuffer);
  rec.setReadName(mHasNames ? reader.name(seqId) : String.valueOf(seqId));
  rec.setReferenceName("*");
  rec.setFlags(flags);
  rec.setReadBases(DnaUtils.bytesToSequenceIncCG(dataBuffer, 0, length).getBytes());
  if (mReader.hasQualityData()) {
    reader.readQuality(seqId, qualityBuffer);
    rec.setBaseQualities(Arrays.copyOf(qualityBuffer, length));
  }
  rec.setInferredInsertSize(0);
  rec.setMateAlignmentStart(0);
  rec.setMateReferenceName("*");
  if (mReadGroupRecord != null) {
    rec.setAttribute(ReadGroupUtils.RG_ATTRIBUTE, mReadGroupRecord.getReadGroupId());
  }

  mWriter.addAlignment(rec);
}
 

private static GATKRead makeRead(final SAMFileHeader header, final String contig, final int readStart, final int fragmentLength, final byte[] bases, final byte[] qual, final String cigar) {
    final SAMRecord read = new SAMRecord(header);
    read.setReferenceName(contig);
    read.setAlignmentStart(readStart);
    read.setReadPairedFlag(true);
    read.setReadUnmappedFlag(false);
    read.setMateUnmappedFlag(false);
    read.setMateReferenceName("Mate");
    read.setReadNegativeStrandFlag(false);
    read.setMateNegativeStrandFlag(true);
    read.setReadBases(bases);
    read.setBaseQualities(qual);


    final int mateStart = readStart + fragmentLength - bases.length;
    read.setMateAlignmentStart(mateStart);
    read.setInferredInsertSize(fragmentLength);
    read.setProperPairFlag(true);
    read.setCigarString(cigar);


    return new SAMRecordToGATKReadAdapter(read);
}
 

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

public static SAMRecord makeTestRecord(SAMFileHeader sfh, int flags, String sequence, String mateSequence, int pos, int matePos, String readGroup) {
  final SAMRecord sr = new SAMRecord(sfh);
  sr.setFlags(flags);
  sr.setReferenceName(sequence);
  sr.setMateReferenceName(mateSequence);
  sr.setAlignmentStart(pos);
  sr.setMateAlignmentStart(matePos);
  sr.setAttribute("RG", readGroup);
  if (pos < matePos) {
    sr.setInferredInsertSize(matePos - pos + 5);
  } else {
    sr.setInferredInsertSize(matePos - pos  - 5); //-1 * (pos - matePos + 5));
  }
  return sr;
}
 

@Test(dataProvider = "topStrandDataProvider")
public void testIsTopStrand(final int referenceIndex, final int alignmentStart, final int mateReferenceIndex, final int mateAlignmentStart,
                            final boolean firstOfPairFlag, final boolean negativeStrandFlag, final boolean mateNegativeStrandFlag,
                            final boolean mapped, final boolean mateMapped,
                            final UmiUtil.ReadStrand strand) {

    final int readLength = 15;
    final int contigLength = 500;
    final SAMFileHeader header = new SAMFileHeader();
    final SAMSequenceDictionary sequenceDictionary = new SAMSequenceDictionary();

    sequenceDictionary.addSequence(new SAMSequenceRecord("chr1", contigLength));
    sequenceDictionary.addSequence(new SAMSequenceRecord("chr2", contigLength));

    header.setSequenceDictionary(sequenceDictionary);

    final SAMRecord rec = new SAMRecord(header);

    rec.setReadUnmappedFlag(!mapped);
    rec.setMateUnmappedFlag(!mateMapped);
    rec.setReadPairedFlag(true);

    rec.setCigarString(readLength + "M");
    rec.setAttribute("MC", readLength + "M");

    rec.setReferenceIndex(referenceIndex);
    rec.setAlignmentStart(alignmentStart);

    rec.setMateReferenceIndex(mateReferenceIndex);
    rec.setMateAlignmentStart(mateAlignmentStart);

    rec.setFirstOfPairFlag(firstOfPairFlag);
    rec.setReadNegativeStrandFlag(negativeStrandFlag);
    rec.setMateNegativeStrandFlag(mateNegativeStrandFlag);

    Assert.assertEquals(UmiUtil.getStrand(rec), strand);
}
 

@Test
public void test() {
	SAMFileHeader header = new SAMFileHeader();
	SAMRecord r1 = new SAMRecord(header);
	r1.setReadName("readName");
	r1.setFlags(4);
	r1.setReferenceIndex(SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX);
	r1.setAlignmentStart(SAMRecord.NO_ALIGNMENT_START);
	r1.setMappingQuality(SAMRecord.NO_MAPPING_QUALITY);
	r1.setCigar(new Cigar());
	r1.setMateReferenceIndex(SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX);
	r1.setMateAlignmentStart(SAMRecord.NO_ALIGNMENT_START);
	r1.setReadBases("A".getBytes());
	r1.setBaseQualityString("!");

	BAMRecordView view = new BAMRecordView(new byte[1024]);
	translate(r1, view);
	r1.setReadName("2");
	translate(r1, view);

	List<SAMRecord> list = toSAMRecord(view, header);
	assertEquals(2, list.size());

	Iterator<SAMRecord> iterator = list.iterator();
	SAMRecord r2 = iterator.next();
	r1.setReadName("readName");
	compare(r1, r2);

	r1.setReadName("2");
	r2 = iterator.next();
	compare(r1, r2);
}
 

@Override
protected int doWork() {
	IOUtil.assertFileIsReadable(INPUT);
	IOUtil.assertFileIsWritable(OUTPUT);
	buildPatterns();

	SamReader in = SamReaderFactory.makeDefault().open(INPUT);

	SAMFileHeader fileHeader = editSequenceDictionary(in.getFileHeader().clone());
	SamHeaderUtil.addPgRecord(fileHeader, this);
	SAMFileWriter out = new SAMFileWriterFactory().makeSAMOrBAMWriter(fileHeader, true, OUTPUT);
	ProgressLogger progLog=new ProgressLogger(log);

	final boolean sequencesRemoved = fileHeader.getSequenceDictionary().getSequences().size() != in.getFileHeader().getSequenceDictionary().getSequences().size();
	
	FilteredReadsMetric m = new FilteredReadsMetric();
	
	for (final SAMRecord r : in) {
		progLog.record(r);
		if (!filterRead(r)) {
               String sequenceName = stripReferencePrefix(r.getReferenceName());
               String mateSequenceName = null;
               if (r.getMateReferenceIndex() != -1)
				mateSequenceName = stripReferencePrefix(r.getMateReferenceName());
		    if (sequencesRemoved || sequenceName != null) {
		        if (sequenceName == null)
					sequenceName = r.getReferenceName();
                   // Even if sequence name has not been edited, if sequences have been removed, need to set
                   // reference name again to invalidate reference index cache.
                   r.setReferenceName(sequenceName);
               }
               if (r.getMateReferenceIndex() != -1 && (sequencesRemoved || mateSequenceName != null)) {
                   if (mateSequenceName == null)
					mateSequenceName = r.getMateReferenceName();
                   // It's possible that the mate was mapped to a reference sequence that has been dropped from
                   // the sequence dictionary.  If so, set the mate to be unmapped.
                   if (fileHeader.getSequenceDictionary().getSequence(mateSequenceName) != null)
					r.setMateReferenceName(mateSequenceName);
				else {
                       r.setMateUnmappedFlag(true);
                       r.setMateReferenceIndex(SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX);
                       r.setMateAlignmentStart(0);
                   }
               }
			out.addAlignment(r);
			m.READS_ACCEPTED++;
		} else {
			m.READS_REJECTED++;
		}
	}
	// write the summary if the summary file is not null.
	writeSummary(this.SUMMARY, m);
       CloserUtil.close(in);
       out.close();
       FilterProgramUtils.reportAndCheckFilterResults("reads", m.READS_ACCEPTED, m.READS_REJECTED,
			PASSING_READ_THRESHOLD, log);
	return (0);
}
 

/**
 * Takes an individual SAMRecord and applies the set of changes/reversions to it that
 * have been requested by program level options.
 */
public void revertSamRecord(final SAMRecord rec) {
    if (RESTORE_ORIGINAL_QUALITIES) {
        final byte[] oq = rec.getOriginalBaseQualities();
        if (oq != null) {
            rec.setBaseQualities(oq);
            rec.setOriginalBaseQualities(null);
        }
    }

    if (REMOVE_DUPLICATE_INFORMATION) {
        rec.setDuplicateReadFlag(false);
    }

    if (REMOVE_ALIGNMENT_INFORMATION) {
        if (rec.getReadNegativeStrandFlag()) {
            rec.reverseComplement(true);
            rec.setReadNegativeStrandFlag(false);
        }

        // Remove all alignment based information about the read itself
        rec.setReferenceIndex(SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX);
        rec.setAlignmentStart(SAMRecord.NO_ALIGNMENT_START);
        rec.setCigarString(SAMRecord.NO_ALIGNMENT_CIGAR);
        rec.setMappingQuality(SAMRecord.NO_MAPPING_QUALITY);

        rec.setInferredInsertSize(0);
        rec.setNotPrimaryAlignmentFlag(false);
        rec.setProperPairFlag(false);
        rec.setReadUnmappedFlag(true);

        // Then remove any mate flags and info related to alignment
        rec.setMateAlignmentStart(SAMRecord.NO_ALIGNMENT_START);
        rec.setMateNegativeStrandFlag(false);
        rec.setMateReferenceIndex(SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX);
        rec.setMateUnmappedFlag(rec.getReadPairedFlag());

        if (RESTORE_HARDCLIPS) {
            String hardClippedBases = rec.getStringAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASES_TAG);
            String hardClippedQualities = rec.getStringAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASE_QUALITIES_TAG);
            if (hardClippedBases != null && hardClippedQualities != null) {
                // Record has already been reverse complemented if this was on the negative strand
                rec.setReadString(rec.getReadString() + hardClippedBases);
                rec.setBaseQualities(SAMUtils.fastqToPhred(SAMUtils.phredToFastq(rec.getBaseQualities()) + hardClippedQualities));

                // Remove hard clipping storage tags
                rec.setAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASES_TAG, null);
                rec.setAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASE_QUALITIES_TAG, null);
            }
        }

        // And then remove any tags that are calculated from the alignment
        ATTRIBUTE_TO_CLEAR.forEach(tag -> rec.setAttribute(tag, null));
    }
}
 

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

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

}
 

@Test
public void test() throws IOException {
	byte[] buf = new byte[1024];
	BAMRecordView view = new BAMRecordView(buf);
	view.setRefID(0);
	view.setAlignmentStart(77);
	view.setMappingScore(44);
	view.setIndexBin(99);
	view.setFlags(555);
	view.setMateRefID(0);
	view.setMateAlStart(78);
	view.setInsertSize(133);

	view.setReadName("name1");
	view.setCigar(TextCigarCodec.decode("10M"));
	view.setBases("AAAAAAAAAA".getBytes());
	view.setQualityScores("BBBBBBBBBB".getBytes());

	int id = 'A' << 16 | 'M' << 8 | 'A';
	view.addTag(id, "Q".getBytes(), 0, 1);

	int len = view.finish();

	System.out.println(Arrays.toString(Arrays.copyOf(buf, len)));

	ByteArrayOutputStream baos = new ByteArrayOutputStream();

	SAMFileHeader header = new SAMFileHeader();
	header.addSequence(new SAMSequenceRecord("14", 14));

	ByteArrayOutputStream baos2 = new ByteArrayOutputStream();
	SAMFileWriter writer = new SAMFileWriterFactory().makeBAMWriter(header, true, baos2);
	SAMRecord record = new SAMRecord(header);
	record.setReferenceIndex(0);
	record.setAlignmentStart(1);
	record.setCigarString("10M");
	record.setFlags(555);
	record.setMappingQuality(44);
	record.setMateReferenceIndex(0);
	record.setMateAlignmentStart(0);
	record.setInferredInsertSize(133);
	record.setReadName("name1");
	record.setReadBases("AAAAAAAAAA".getBytes());
	record.setBaseQualities("BBBBBBBBBB".getBytes());
	record.setAttribute("AM", 'Q');

	System.out.println("BAMFileWriter.addAlignment():");
	writer.addAlignment(record);
	System.out.println(".");
	writer.close();

	System.out.println("------------------------------------------");
	System.out.println();
	System.out.println(new String(baos2.toByteArray()));
	System.out.println();

	SAMFileReader.setDefaultValidationStringency(ValidationStringency.SILENT);
	SAMFileReader reader2 = new SAMFileReader(new ByteArrayInputStream(baos2.toByteArray()));
	SAMRecordIterator iterator = reader2.iterator();
	while (iterator.hasNext()) {
		record = iterator.next();
		System.out.println(record.getSAMString());
	}
	System.out.println("------------------------------------------");

	BlockCompressedOutputStream bcos = new BlockCompressedOutputStream(baos, null);
	bcos.write("BAM\1".getBytes());
	bcos.write(toByteArray(header));
	CramInt.writeInt32(header.getSequenceDictionary().size(), bcos);
	for (final SAMSequenceRecord sequenceRecord : header.getSequenceDictionary().getSequences()) {
		byte[] bytes = sequenceRecord.getSequenceName().getBytes();
		CramInt.writeInt32(bytes.length + 1, bcos);
		bcos.write(sequenceRecord.getSequenceName().getBytes());
		bcos.write(0);
		CramInt.writeInt32(sequenceRecord.getSequenceLength(), bcos);
	}
	bcos.write(buf, 0, len);
	bcos.close();

	System.out.println(new String(baos.toByteArray()));

	SAMFileReader reader = new SAMFileReader(new ByteArrayInputStream(baos.toByteArray()));
	iterator = reader.iterator();
	while (iterator.hasNext()) {
		record = iterator.next();
		System.out.println(record.getSAMString());
	}
	reader.close();

}