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

private SAMRecord createSamRecord(final SAMFileHeader header, final String baseName, final FastqRecord frec, final boolean paired) {
    final SAMRecord srec = new SAMRecord(header);
    srec.setReadName(baseName);
    srec.setReadString(frec.getReadString());
    srec.setReadUnmappedFlag(true);
    srec.setAttribute(ReservedTagConstants.READ_GROUP_ID, READ_GROUP_NAME);
    final byte[] quals = StringUtil.stringToBytes(frec.getBaseQualityString());
    convertQuality(quals, QUALITY_FORMAT);
    for (final byte qual : quals) {
        final int uQual = qual & 0xff;
        if (uQual < MIN_Q || uQual > MAX_Q) {
            throw new PicardException("Base quality " + uQual + " is not in the range " + MIN_Q + ".." +
            MAX_Q + " for read " + frec.getReadHeader());
        }
    }
    srec.setBaseQualities(quals);

    if (paired) {
        srec.setReadPairedFlag(true);
        srec.setMateUnmappedFlag(true);
    }
    return srec ;
}
 

private void addAlignmentsForBestFragmentMapqStrategy(
        final SAMFileWriter writer, final SAMRecord unmappedRecord, final String sequence, final int[] mapqs) {
    boolean reverse = false;
    int alignmentStart = 1;
    for (final int mapq : mapqs) {
        final SAMRecord alignedRecord = new SAMRecord(writer.getFileHeader());
        alignedRecord.setReadName(unmappedRecord.getReadName());
        alignedRecord.setReadBases(unmappedRecord.getReadBases());
        alignedRecord.setBaseQualities(unmappedRecord.getBaseQualities());
        alignedRecord.setReferenceName(sequence);
        alignedRecord.setAlignmentStart(alignmentStart);
        alignmentStart += 10; // Any old position will do
        alignedRecord.setReadNegativeStrandFlag(reverse);
        reverse = !reverse;
        alignedRecord.setCigarString(unmappedRecord.getReadBases().length + "M");
        alignedRecord.setMappingQuality(mapq);
        alignedRecord.setReadPairedFlag(unmappedRecord.getReadPairedFlag());
        alignedRecord.setFirstOfPairFlag(unmappedRecord.getFirstOfPairFlag());
        alignedRecord.setSecondOfPairFlag(unmappedRecord.getSecondOfPairFlag());
        alignedRecord.setMateUnmappedFlag(true);
        writer.addAlignment(alignedRecord);
    }
}
 

private void addAlignmentForMostStrategy(
        final SAMFileWriter writer, final SAMRecord unmappedRecord, final MostDistantStrategyAlignmentSpec spec,
        final boolean reverse) {
    final SAMRecord alignedRecord = new SAMRecord(writer.getFileHeader());
    alignedRecord.setReadName(unmappedRecord.getReadName());
    alignedRecord.setReadBases(unmappedRecord.getReadBases());
    alignedRecord.setBaseQualities(unmappedRecord.getBaseQualities());
    alignedRecord.setReferenceName(spec.sequence);
    alignedRecord.setAlignmentStart(spec.alignmentStart);
    alignedRecord.setReadNegativeStrandFlag(reverse);
    alignedRecord.setCigarString(unmappedRecord.getReadBases().length + "M");
    alignedRecord.setMappingQuality(spec.mapQ);
    alignedRecord.setReadPairedFlag(unmappedRecord.getReadPairedFlag());
    alignedRecord.setFirstOfPairFlag(unmappedRecord.getFirstOfPairFlag());
    alignedRecord.setSecondOfPairFlag(unmappedRecord.getSecondOfPairFlag());
    alignedRecord.setMateUnmappedFlag(true);
    writer.addAlignment(alignedRecord);
}
 

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

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

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

/**
 * Creates a new SAM record from the basecall data
 */
private SAMRecord createSamRecord(final ReadData readData, final String readName, final boolean isPf, final boolean firstOfPair,
                                  final String unmatchedBarcode, final String barcodeQuality,
                                  final List<String> molecularIndexes, final List<String> molecularIndexQualities) {
    final SAMRecord sam = new SAMRecord(null);
    sam.setReadName(readName);
    sam.setReadBases(readData.getBases());
    sam.setBaseQualities(readData.getQualities());

    // Flag values
    sam.setReadPairedFlag(isPairedEnd);
    sam.setReadUnmappedFlag(true);
    sam.setReadFailsVendorQualityCheckFlag(!isPf);
    if (isPairedEnd) {
        sam.setMateUnmappedFlag(true);
        sam.setFirstOfPairFlag(firstOfPair);
        sam.setSecondOfPairFlag(!firstOfPair);
    }

    if (filters.filterOut(sam)) {
        sam.setAttribute(ReservedTagConstants.XN, 1);
    }

    if (this.readGroupId != null) {
        sam.setAttribute(SAMTag.RG.name(), readGroupId);
    }

    // If it's a barcoded run and it has been decided that the original BC value should be added to the record, do it
    if (unmatchedBarcode != null) {
        sam.setAttribute(SAMTag.BC.name(), unmatchedBarcode);
        if (barcodeQuality != null ) {
            sam.setAttribute(SAMTag.QT.name(), barcodeQuality);
        }
    }

    if (!molecularIndexes.isEmpty()) {
        if (!this.MOLECULAR_INDEX_TAG.isEmpty()) {
            sam.setAttribute(this.MOLECULAR_INDEX_TAG, String.join(MOLECULAR_INDEX_, molecularIndexes));
        }
        if (!this.MOLECULAR_INDEX_QUALITY_TAG.isEmpty()) {
            sam.setAttribute(this.MOLECULAR_INDEX_QUALITY_TAG, String.join(MOLECULAR_INDEX_, molecularIndexQualities));
        }
        if (!this.tagPerMolecularIndex.isEmpty()) {
            if (tagPerMolecularIndex.size() != molecularIndexes.size()) {
                throw new PicardException("Found " + molecularIndexes.size() + " molecular indexes but only " + tagPerMolecularIndex.size() + " SAM tags given.");
            }
            for (int i = 0; i < this.tagPerMolecularIndex.size(); i++) {
                sam.setAttribute(this.tagPerMolecularIndex.get(i), molecularIndexes.get(i));
            }
        }
    }

    return sam;
}
 

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

/**
 * Read a single paired-end read from a file, and create one or more aligned records for the read pair based on
 * the lists, merge with the original paired-end read, and assert expected results.
 * @param description
 * @param firstOfPair List that describes the aligned SAMRecords to create for the first end.
 * @param secondOfPair List that describes the aligned SAMRecords to create for the second end.
 * @param expectedPrimaryHitIndex Expected value for the HI tag in the primary alignment in the merged output.
 * @param expectedNumFirst Expected number of first-of-pair SAMRecords in the merged output.
 * @param expectedNumSecond Expected number of second-of-pair SAMRecords in the merged output.
 * @param expectedPrimaryMapq Sum of mapqs of both ends of primary alignment in the merged output.
 * @throws Exception
 */
@Test(dataProvider = "testPairedMultiHitWithFilteringTestCases")
public void testPairedMultiHitWithFiltering(final String description, final List<HitSpec> firstOfPair, final List<HitSpec> secondOfPair,
                                            final Integer expectedPrimaryHitIndex, final int expectedNumFirst,
                                            final int expectedNumSecond, final int expectedPrimaryMapq) throws Exception {

    // Create the aligned file by copying bases, quals, readname from the unmapped read, and conforming to each HitSpec.
    final File unmappedSam = new File(TEST_DATA_DIR, "multihit.filter.unmapped.sam");
    final SAMRecordIterator unmappedSamFileIterator = SamReaderFactory.makeDefault().open(unmappedSam).iterator();
    final SAMRecord firstUnmappedRec = unmappedSamFileIterator.next();
    final SAMRecord secondUnmappedRec = unmappedSamFileIterator.next();
    unmappedSamFileIterator.close();
    final File alignedSam = File.createTempFile("aligned.", ".sam");
    alignedSam.deleteOnExit();
    final SAMFileHeader alignedHeader = new SAMFileHeader();
    alignedHeader.setSortOrder(SAMFileHeader.SortOrder.queryname);
    alignedHeader.setSequenceDictionary(SamReaderFactory.makeDefault().getFileHeader(sequenceDict).getSequenceDictionary());
    final SAMFileWriter alignedWriter = new SAMFileWriterFactory().makeSAMWriter(alignedHeader, true, alignedSam);
    for (int i = 0; i < Math.max(firstOfPair.size(), secondOfPair.size()); ++i) {
        final HitSpec firstHitSpec = firstOfPair.isEmpty()? null: firstOfPair.get(i);
        final HitSpec secondHitSpec = secondOfPair.isEmpty()? null: secondOfPair.get(i);
        final SAMRecord first = makeRead(alignedHeader, firstUnmappedRec, firstHitSpec, true, i);
        final SAMRecord second = makeRead(alignedHeader, secondUnmappedRec, secondHitSpec, false, i);
        if (first != null && second != null) SamPairUtil.setMateInfo(first, second);
        if (first != null) {
            if (second == null) first.setMateUnmappedFlag(true);
            alignedWriter.addAlignment(first);
        }
        if (second != null) {
            if (first == null) second.setMateUnmappedFlag(true);
            alignedWriter.addAlignment(second);
        }
    }
    alignedWriter.close();

    // Merge aligned file with original unmapped file.
    final File mergedSam = File.createTempFile("merged.", ".sam");
    mergedSam.deleteOnExit();

    doMergeAlignment(unmappedSam, Collections.singletonList(alignedSam),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true, fasta, mergedSam,
            null, null, null, null, null, null);

    assertSamValid(mergedSam);

    // Tally metrics and check for agreement with expected.
    final SamReader mergedReader = SamReaderFactory.makeDefault().open(mergedSam);
    int numFirst = 0;
    int numSecond = 0;
    Integer primaryHitIndex = null;
    int primaryMapq = 0;
    for (final SAMRecord rec : mergedReader) {
        if (rec.getFirstOfPairFlag()) ++numFirst;
        if (rec.getSecondOfPairFlag()) ++numSecond;
        if (!rec.getNotPrimaryAlignmentFlag()  && !rec.getReadUnmappedFlag()) {
            final Integer hitIndex = rec.getIntegerAttribute(SAMTag.HI.name());
            final int newHitIndex = (hitIndex == null? -1: hitIndex);
            if (primaryHitIndex == null) primaryHitIndex = newHitIndex;
            else Assert.assertEquals(newHitIndex, primaryHitIndex.intValue());
            primaryMapq += rec.getMappingQuality();
        }
    }
    Assert.assertEquals(numFirst, expectedNumFirst);
    Assert.assertEquals(numSecond, expectedNumSecond);
    Assert.assertEquals(primaryHitIndex, expectedPrimaryHitIndex);
    Assert.assertEquals(primaryMapq, expectedPrimaryMapq);
}
 

private void testBestFragmentMapqStrategy(final String testName, final int[] firstMapQs, final int[] secondMapQs,
                                          final boolean includeSecondary, final int expectedFirstMapq,
                                          final int expectedSecondMapq) throws Exception {
    final File unmappedSam = File.createTempFile("unmapped.", ".sam");
    unmappedSam.deleteOnExit();
    final SAMFileWriterFactory factory = new SAMFileWriterFactory();
    final SAMFileHeader header = new SAMFileHeader();
    header.setSortOrder(SAMFileHeader.SortOrder.queryname);

    final String readName = "theRead";
    final SAMRecord firstUnmappedRead = new SAMRecord(header);
    firstUnmappedRead.setReadName(readName);
    firstUnmappedRead.setReadString("ACGTACGTACGTACGT");
    firstUnmappedRead.setBaseQualityString("5555555555555555");
    firstUnmappedRead.setReadUnmappedFlag(true);
    firstUnmappedRead.setMateUnmappedFlag(true);
    firstUnmappedRead.setReadPairedFlag(true);
    firstUnmappedRead.setFirstOfPairFlag(true);

    final SAMRecord secondUnmappedRead = new SAMRecord(header);
    secondUnmappedRead.setReadName(readName);
    secondUnmappedRead.setReadString("TCGAACGTTCGAACTG");
    secondUnmappedRead.setBaseQualityString("6666666666666666");
    secondUnmappedRead.setReadUnmappedFlag(true);
    secondUnmappedRead.setMateUnmappedFlag(true);
    secondUnmappedRead.setReadPairedFlag(true);
    secondUnmappedRead.setSecondOfPairFlag(true);

    final SAMFileWriter unmappedWriter = factory.makeSAMWriter(header, false, unmappedSam);
    unmappedWriter.addAlignment(firstUnmappedRead);
    unmappedWriter.addAlignment(secondUnmappedRead);
    unmappedWriter.close();

    final File alignedSam = File.createTempFile("aligned.", ".sam");
    alignedSam.deleteOnExit();

    final String sequence = "chr1";
    // Populate the header with SAMSequenceRecords
    header.setSequenceDictionary(SAMSequenceDictionaryExtractor.extractDictionary(sequenceDict2.toPath()));

    final SAMFileWriter alignedWriter = factory.makeSAMWriter(header, false, alignedSam);

    addAlignmentsForBestFragmentMapqStrategy(alignedWriter, firstUnmappedRead, sequence, firstMapQs);
    addAlignmentsForBestFragmentMapqStrategy(alignedWriter, secondUnmappedRead, sequence, secondMapQs);
    alignedWriter.close();

    final File output = File.createTempFile("testBestFragmentMapqStrategy." + testName, ".sam");
    output.deleteOnExit();
    doMergeAlignment(unmappedSam, Collections.singletonList(alignedSam),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true,
            fasta, output,
            SamPairUtil.PairOrientation.FR,
            MergeBamAlignment.PrimaryAlignmentStrategy.BestEndMapq,
            null, includeSecondary, null, null);
    final SamReader reader = SamReaderFactory.makeDefault().open(output);

    int numFirstRecords = 0;
    int numSecondRecords = 0;
    int firstPrimaryMapq = -1;
    int secondPrimaryMapq = -1;
    for (final SAMRecord rec: reader) {
        Assert.assertTrue(rec.getReadPairedFlag());
        if (rec.getFirstOfPairFlag()) ++numFirstRecords;
        else if (rec.getSecondOfPairFlag()) ++ numSecondRecords;
        else Assert.fail("unpossible!");
        if (!rec.getReadUnmappedFlag() && !rec.getNotPrimaryAlignmentFlag()) {
            if (rec.getFirstOfPairFlag()) {
                Assert.assertEquals(firstPrimaryMapq, -1);
                firstPrimaryMapq = rec.getMappingQuality();
            } else {
                Assert.assertEquals(secondPrimaryMapq, -1);
                secondPrimaryMapq = rec.getMappingQuality();
            }
        } else if (rec.getNotPrimaryAlignmentFlag()) {
            Assert.assertTrue(rec.getMateUnmappedFlag());
        }
    }
    reader.close();
    Assert.assertEquals(firstPrimaryMapq, expectedFirstMapq);
    Assert.assertEquals(secondPrimaryMapq, expectedSecondMapq);
    if (!includeSecondary) {
        Assert.assertEquals(numFirstRecords, 1);
        Assert.assertEquals(numSecondRecords, 1);
    } else {
        // If no alignments for an end, there will be a single unmapped record
        Assert.assertEquals(numFirstRecords, Math.max(1, firstMapQs.length));
        Assert.assertEquals(numSecondRecords, Math.max(1, secondMapQs.length));
    }
}
 

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

}