htsjdk.samtools.SamPairUtil Java Examples

public PerUnitInsertSizeMetricsCollector(final String sample, final String library, final String readGroup) {
    this.sample = sample;
    this.library = library;
    this.readGroup = readGroup;
    String prefix = null;
    if (this.readGroup != null) {
        prefix = this.readGroup + ".";
    }
    else if (this.library != null) {
        prefix = this.library + ".";
    }
    else if (this.sample != null) {
        prefix = this.sample + ".";
    }
    else {
        prefix = "All_Reads.";
    }
    histograms.put(SamPairUtil.PairOrientation.FR,     new Histogram<Integer>("insert_size", prefix + "fr_count"));
    histograms.put(SamPairUtil.PairOrientation.TANDEM, new Histogram<Integer>("insert_size", prefix + "tandem_count"));
    histograms.put(SamPairUtil.PairOrientation.RF,     new Histogram<Integer>("insert_size", prefix + "rf_count"));
}
 

public PerUnitInsertSizeMetricsCollector(
        final String sample,
        final String library,
        final String readGroup,
        final double minimumPct,
        final double deviations,
        final Integer histogramWidth) {

    this.sample = sample;
    this.library = library;
    this.readGroup = readGroup;
    this.minimumPct = minimumPct;
    this.histogramWidth = histogramWidth;
    this.deviations = deviations;

    String prefix = createHistogramValuePrefix();

    histograms.put(SamPairUtil.PairOrientation.FR,     new Histogram<>("insert_size", prefix + "fr_count"));
    histograms.put(SamPairUtil.PairOrientation.TANDEM, new Histogram<>("insert_size", prefix + "tandem_count"));
    histograms.put(SamPairUtil.PairOrientation.RF,     new Histogram<>("insert_size", prefix + "rf_count"));
}
 

@Test(dataProvider="data")
public void testSortingOnSamAlignmentMerger(final File unmapped, final File aligned, final boolean sorted, final boolean coordinateSorted, final String testName)
    throws IOException {

    final File target = File.createTempFile("target", "bam");
    target.deleteOnExit();
    final SamAlignmentMerger merger = new SamAlignmentMerger(unmapped,  target, fasta, null, true, false,
            false, Collections.singletonList(aligned), 1, null, null, null, null, null, null,
            Collections.singletonList(SamPairUtil.PairOrientation.FR),
            coordinateSorted ? SAMFileHeader.SortOrder.coordinate : SAMFileHeader.SortOrder.queryname,
            new BestMapqPrimaryAlignmentSelectionStrategy(), false, false, 30);

    merger.mergeAlignment(Defaults.REFERENCE_FASTA);
    Assert.assertEquals(sorted, !merger.getForceSort());
    final SAMRecordIterator it = SamReaderFactory.makeDefault().open(target).iterator();
    int aln = 0;
    while (it.hasNext()) {
        final SAMRecord rec = it.next();
        if (!rec.getReadUnmappedFlag()) {
            aln++;
        }
    }
    Assert.assertEquals(aln, 6, "Incorrect number of aligned reads in merged BAM file");
}
 

private void doMergeAlignment(final File unmappedBam, final List<File> alignedBams,
                              final List<File> read1AlignedBams, final List<File> read2AlignedBams, final Integer read1Trim, final Integer read2Trim,
                              final boolean alignReadsOnly, final boolean clipAdapters, final boolean isBisulfiteSequence, final int maxInsOrDels,
                              final String progRecordId, final String progGroupVersion, final String progGroupCommandLine, final String progGroupName,
                              final boolean pairedRun, final File refSeq, final File output,
                              final SamPairUtil.PairOrientation expectedOrientation, final MergeBamAlignment.PrimaryAlignmentStrategy primaryAlignmentStrategy,
                              final String attributesToRetain,
                              final Boolean includeSecondary,
                              final Boolean unmapContaminantReads,
                              final SAMFileHeader.SortOrder sortOrder) {
    doMergeAlignment(unmappedBam, alignedBams, read1AlignedBams, read2AlignedBams, read1Trim, read2Trim,
            alignReadsOnly, clipAdapters, isBisulfiteSequence, maxInsOrDels,
            progRecordId, progGroupVersion, progGroupCommandLine, progGroupName, pairedRun, refSeq, output,
            expectedOrientation, primaryAlignmentStrategy, attributesToRetain, includeSecondary, unmapContaminantReads,
            sortOrder, null);
}
 

@Override
protected InsertSizeCollectorArgs makeArg(SAMRecord samRecord, ReferenceSequence refSeq) {
    final int insertSize = Math.abs(samRecord.getInferredInsertSize());
    final SamPairUtil.PairOrientation orientation = SamPairUtil.getPairOrientation(samRecord);

    return new InsertSizeCollectorArgs(insertSize, orientation);
}
 

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

@Override
protected InsertSizeMetricsCollectorArgs makeArg(SAMRecord samRecord, ReferenceSequence refSeq) {
    // inferred insert size is negative if the mate maps to lower position than the read, so use abs
    final int insertSize = Math.abs(samRecord.getInferredInsertSize());
    final SamPairUtil.PairOrientation orientation = SamPairUtil.getPairOrientation(samRecord);

    return new InsertSizeMetricsCollectorArgs(insertSize, orientation);
}
 

@Test
public void testRemoveNmMdAndUqOnOverlappingReads() throws IOException {
    final File output = File.createTempFile("testRemoveNmMdAndUqOnOverlappingReads", ".sam");
    output.deleteOnExit();
    doMergeAlignment(new File(TEST_DATA_DIR, "removetags.unmapped.sam"),
            Collections.singletonList(new File(TEST_DATA_DIR, "removetags.aligned.sam")),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true,
            new File(TEST_DATA_DIR, "removetags.fasta"), output,
            SamPairUtil.PairOrientation.FR, null,
            null, null, null, SAMFileHeader.SortOrder.queryname);

    final SamReader result = SamReaderFactory.makeDefault().open(output);
    for (final SAMRecord rec : result) {
        boolean hasTags = false;
        if (rec.getReadName().startsWith("CLIPPED")) {
            final String[] readNameFields = rec.getReadName().split(":");
            final int index = rec.getFirstOfPairFlag() ? 1 : 2;
            hasTags = Integer.parseInt(readNameFields[index]) == 1;
        }
        if (hasTags) {
            Assert.assertNull(rec.getAttribute("MD"));
            Assert.assertNull(rec.getAttribute("NM"));
        } else {
            Assert.assertNotNull(rec.getAttribute("MD"));
            Assert.assertNotNull(rec.getAttribute("NM"));
        }
    }
    result.close();
}
 

@Test(expectedExceptions = {IllegalStateException.class, PicardException.class})
public void testOldQuerynameSortFails() throws IOException {

    final File merged = File.createTempFile("merged", BamFileIoUtils.BAM_FILE_EXTENSION);
    merged.deleteOnExit();

    doMergeAlignment(badorderUnmappedBam, Collections.singletonList(badorderAlignedBam),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true, fasta, merged,
            SamPairUtil.PairOrientation.FR, null,
            null, null, null, null);
    Assert.fail("Merger should have failed because unmapped reads are not in queryname order but didn't");
}
 

/**
 * Test that clipping of FR reads for fragments shorter than read length happens only when it should.
 */
@Test
public void testShortFragmentClipping() throws Exception {
    final File output = File.createTempFile("testShortFragmentClipping", ".sam");
    output.deleteOnExit();
    doMergeAlignment(new File(TEST_DATA_DIR, "cliptest.unmapped.sam"),
            Collections.singletonList(new File(TEST_DATA_DIR, "cliptest.aligned.sam")),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true,
            new File(TEST_DATA_DIR, "cliptest.fasta"), output,
            SamPairUtil.PairOrientation.FR, null,
            null, null, null, null);

    final SamReader result = SamReaderFactory.makeDefault().open(output);
    final Map<String, SAMRecord> firstReadEncountered = new HashMap<String, SAMRecord>();

    for (final SAMRecord rec : result) {
        final SAMRecord otherEnd = firstReadEncountered.get(rec.getReadName());
        if (otherEnd == null) {
            firstReadEncountered.put(rec.getReadName(), rec);
        } else {
            final int fragmentStart = Math.min(rec.getAlignmentStart(), otherEnd.getAlignmentStart());
            final int fragmentEnd = Math.max(rec.getAlignmentEnd(), otherEnd.getAlignmentEnd());
            final String[] readNameFields = rec.getReadName().split(":");
            // Read name of each pair includes the expected fragment start and fragment end positions.
            final int expectedFragmentStart = Integer.parseInt(readNameFields[1]);
            final int expectedFragmentEnd = Integer.parseInt(readNameFields[2]);
            Assert.assertEquals(fragmentStart, expectedFragmentStart, rec.getReadName());
            Assert.assertEquals(fragmentEnd, expectedFragmentEnd, rec.getReadName());
        }
    }
    result.close();
}
 

/**
 * Various scenarios for EarliestFragmentStrategy.  Confirms that one of the expected ones is selected.
 * Note that there may be an arbitrary selection due to a tie.
 */
@Test(dataProvider = "testEarliestFragmentStrategyDataProvider")
public void testEarliestFragmentStrategy(final String testName, final MultipleAlignmentSpec[] specs) throws IOException {

    final File output = File.createTempFile(testName, ".sam");
    output.deleteOnExit();
    final File[] sams = createSamFilesToBeMerged(specs);

    doMergeAlignment(sams[0], Collections.singletonList(sams[1]),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true, fasta, output,
            SamPairUtil.PairOrientation.FR, MergeBamAlignment.PrimaryAlignmentStrategy.EarliestFragment,
            ONE_OF_THE_BEST_TAG,
            null, false, null);


    final SamReader mergedReader = SamReaderFactory.makeDefault().open(output);
    boolean seenPrimary = false;
    for (final SAMRecord rec : mergedReader) {
        if (!rec.getNotPrimaryAlignmentFlag()) {
            seenPrimary = true;
            final Integer oneOfTheBest = rec.getIntegerAttribute(ONE_OF_THE_BEST_TAG);
            Assert.assertEquals(oneOfTheBest, new Integer(1), "Read not marked as one of the best is primary: " + rec);
        }
    }
    CloserUtil.close(mergedReader);
    Assert.assertTrue(seenPrimary, "Never saw primary alignment");
}
 

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

}
 

public InsertSizeMetricsCollectorArgs(final int insertSize, final SamPairUtil.PairOrientation po) {
    this.insertSize = insertSize;
    this.po = po;
}
 

public SamPairUtil.PairOrientation getPairOrientation() {
    return po;
}
 

@Test
public void testMappedToMultipleStrands() throws Exception {
    final File outputMappedToMultipleStands = File.createTempFile("mappedToMultipleStrands", ".sam");
    outputMappedToMultipleStands.deleteOnExit();

    doMergeAlignment(mergingUnmappedBam,
            Collections.singletonList(multipleStrandsAlignedBam),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true, fasta, outputMappedToMultipleStands,
            SamPairUtil.PairOrientation.FR, null, null, null, null, null);

    final SamReader result = SamReaderFactory.makeDefault().open(outputMappedToMultipleStands);

    for (final SAMRecord sam : result) {
        if (sam.getReadName().equals("test:1") && !sam.getReadUnmappedFlag()) {
            if (sam.getReadNegativeStrandFlag() && sam.getFirstOfPairFlag()) {
                Assert.assertEquals(sam.getReadString(), "TTTACTGATGTTATGACCATTACTCCGAAAGTGCCAAGATCATGAAGGGCAAGGAGAGAGTGGGATCCCCGGGTAC", "Read aligned to negative strand has unexpected bases.");
            } else {
                Assert.assertEquals(sam.getReadString(), "GTACCCGGGGATCCCACTCTCTCCTTGCCCTTCATGATCTTGGCACTTTCGGAGTAATGGTCATAACATCAGTAAA", "Read aligned to positive strand has unexpected bases.");
            }
        }

        if (sam.getReadName().equals("test:2") && !sam.getReadUnmappedFlag()) {
            if (sam.getReadNegativeStrandFlag() && sam.getSecondOfPairFlag()) {
                Assert.assertEquals(sam.getReadString(), "TTATTCACTTAGTGTGTTTTTCCTGAGAACTTGCTATGTGTTAGGTCCTAGGCTGGGTGGGATCCTCTAGAGTCGA", "Read aligned to negative strand has unexpected bases.");
            } else {
                Assert.assertEquals(sam.getReadString(), "TCGACTCTAGAGGATCCCACCCAGCCTAGGACCTAACACATAGCAAGTTCTCAGGAAAAACACACTAAGTGAATAA", "Read aligned to positive strand has unexpected bases.");
            }
        }

        if (sam.getReadName().equals("test:5") && !sam.getReadUnmappedFlag()) {
            if (sam.getReadNegativeStrandFlag()) {
                Assert.assertEquals(sam.getReadString(), "AGTTTTGGTTTGTCAGACCCAGCCCTGGGCACAGATGAGGAATTCTGGCTTCTCCTCCTGTGGGATCCCCGGGTAC", "Read aligned to negative strand has unexpected bases.");
            } else {
                Assert.assertEquals(sam.getReadString(), "GTACCCGGGGATCCCACAGGAGGAGAAGCCAGAATTCCTCATCTGTGCCCAGGGCTGGGTCTGACAAACCAAAACT", "Read aligned to positive strand has unexpected bases.");
            }
        }
    }
}
 

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

/**
 * Confirm that paired reads are rejected by PrimaryAlignmentStrategy.EarliestFragment.
 */
@Test(expectedExceptions = UnsupportedOperationException.class)
public void testEarliestFragmentStrategyPaired() throws Exception {
    final File output = File.createTempFile("mergeTest", ".sam");
    output.deleteOnExit();

    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 cigar = "16M";

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

    final SAMRecord secondOfPair = new SAMRecord(header);
    secondOfPair.setReadName("theRead");
    secondOfPair.setReadString("ACGTACGTACGTACGT");
    secondOfPair.setBaseQualityString("5555555555555555");
    secondOfPair.setReadUnmappedFlag(true);
    secondOfPair.setReadPairedFlag(true);
    secondOfPair.setSecondOfPairFlag(true);
    SamPairUtil.setMateInfo(firstOfPair, secondOfPair);

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

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

    // Populate the header with SAMSequenceRecords
    header.setSequenceDictionary(SAMSequenceDictionaryExtractor.extractDictionary(sequenceDict2.toPath()));

    // Create 2 alignments for each end of pair
    final SAMFileWriter alignedWriter = factory.makeSAMWriter(header, false, alignedSam);
    for (int i = 1; i <= 2; ++i) {
        final SAMRecord firstOfPairAligned = new SAMRecord(header);
        firstOfPairAligned.setReadName(firstOfPair.getReadName());
        firstOfPairAligned.setReadBases(firstOfPair.getReadBases());
        firstOfPairAligned.setBaseQualities(firstOfPair.getBaseQualities());
        firstOfPairAligned.setReferenceName("chr1");
        firstOfPairAligned.setAlignmentStart(i);
        firstOfPairAligned.setCigarString(cigar);
        firstOfPairAligned.setMappingQuality(100);
        firstOfPairAligned.setReadPairedFlag(true);
        firstOfPairAligned.setFirstOfPairFlag(true);
        firstOfPairAligned.setAttribute(SAMTag.HI.name(), i);

        final SAMRecord secondOfPairAligned = new SAMRecord(header);
        secondOfPairAligned.setReadName(secondOfPair.getReadName());
        secondOfPairAligned.setReadBases(secondOfPair.getReadBases());
        secondOfPairAligned.setBaseQualities(secondOfPair.getBaseQualities());
        secondOfPairAligned.setReferenceName("chr1");
        secondOfPairAligned.setAlignmentStart(i + 10);
        secondOfPairAligned.setCigarString(cigar);
        secondOfPairAligned.setMappingQuality(100);
        secondOfPairAligned.setReadPairedFlag(true);
        secondOfPairAligned.setSecondOfPairFlag(true);
        secondOfPairAligned.setAttribute(SAMTag.HI.name(), i);

        SamPairUtil.setMateInfo(firstOfPairAligned, secondOfPairAligned);

        alignedWriter.addAlignment(firstOfPairAligned);
        alignedWriter.addAlignment(secondOfPairAligned);
    }
    alignedWriter.close();

    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.EarliestFragment,
            null, null, null, null);

    Assert.fail("Exception was not thrown");
}
 

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

/**
 * Minimal test of merging data from separate read 1 and read 2 alignments
 */
@Test(dataProvider="separateTrimmed")
public void testMergingFromSeparatedReadTrimmedAlignments(final File unmapped, final List<File> r1Align, final List<File> r2Align, final int r1Trim, final int r2Trim, final String testName) throws Exception {
     final File output = File.createTempFile("mergeMultipleAlignmentsTest", ".sam");
     output.deleteOnExit();

    doMergeAlignment(unmapped, null, r1Align, r2Align, r1Trim, r2Trim,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true, fasta, output,
            SamPairUtil.PairOrientation.FR, null, null, null, null, null);

    SamReaderFactory factory = SamReaderFactory.makeDefault();
    final SamReader result = factory.open(output);
     for (final SAMRecord sam : result) {
        // Get the alignment record
        final List<File> rFiles = sam.getFirstOfPairFlag() ? r1Align : r2Align;
        SAMRecord alignment = null;
        for (final File f : rFiles) {
            for (final SAMRecord tmp : factory.open(f)) {
                if (tmp.getReadName().equals(sam.getReadName())) {
                    alignment = tmp;
                    break;
                }
            }
            if (alignment != null) break;
        }

        final int trim = sam.getFirstOfPairFlag() ? r1Trim : r2Trim;
        final int notWrittenToFastq = sam.getReadLength() - (trim + alignment.getReadLength());
        final int beginning = sam.getReadNegativeStrandFlag() ? notWrittenToFastq : trim;
        final int end = sam.getReadNegativeStrandFlag() ? trim : notWrittenToFastq;

        if (!sam.getReadUnmappedFlag()) {
            final CigarElement firstMergedCigarElement = sam.getCigar().getCigarElement(0);
            final CigarElement lastMergedCigarElement = sam.getCigar().getCigarElement(sam.getCigar().getCigarElements().size() - 1);
            final CigarElement firstAlignedCigarElement = alignment.getCigar().getCigarElement(0);
            final CigarElement lastAlignedCigarElement = alignment.getCigar().getCigarElement(alignment.getCigar().getCigarElements().size() - 1);

            if (beginning > 0) {
                Assert.assertEquals(firstMergedCigarElement.getOperator(), CigarOperator.S, "First element is not a soft clip");
                Assert.assertEquals(firstMergedCigarElement.getLength(), beginning + ((firstAlignedCigarElement.getOperator() == CigarOperator.S) ? firstAlignedCigarElement.getLength() : 0));
            }
            if (end > 0) {
                Assert.assertEquals(lastMergedCigarElement.getOperator(), CigarOperator.S, "Last element is not a soft clip");
                Assert.assertEquals(lastMergedCigarElement.getLength(), end + ((lastAlignedCigarElement.getOperator() == CigarOperator.S) ? lastAlignedCigarElement.getLength() : 0));
            }
        }
    }

}
 

@Test
public void testMergerFromMultipleFiles() throws Exception {
    final File output = File.createTempFile("mergeTest", ".sam");
    output.deleteOnExit();

    doMergeAlignment(unmappedBam, Arrays.asList(oneHalfAlignedBam, otherHalfAlignedBam),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true, fasta, output,
            SamPairUtil.PairOrientation.FR, null, null, null, null, null);

    final SamReader result = SamReaderFactory.makeDefault().open(output);

    for (final SAMRecord sam : result) {
        // This tests that we clip both (a) when the adapter is marked in the unmapped BAM file and
        // (b) when the insert size is less than the read length
        if (sam.getReadName().equals("both_reads_align_clip_adapter") ||
                sam.getReadName().equals("both_reads_align_clip_marked")) {
            Assert.assertEquals(sam.getReferenceName(), "chr7");
            if (sam.getReadNegativeStrandFlag()) {
                Assert.assertEquals(sam.getCigarString(), "5S96M", "Incorrect CIGAR string for " +
                        sam.getReadName());
            } else {
                Assert.assertEquals(sam.getCigarString(), "96M5S", "Incorrect CIGAR string for " +
                        sam.getReadName());
            }
        }
        // This tests that we DON'T clip when we run off the end if there are equal to or more than
        // MIN_ADAPTER_BASES hanging off the end
        else if (sam.getReadName().equals("both_reads_align_min_adapter_bases_exceeded")) {
            Assert.assertEquals(sam.getReferenceName(), "chr7");
            Assert.assertTrue(!sam.getCigarString().contains("S"),
                    "Read was clipped when it should not be.");
        } else if (sam.getReadName().equals("neither_read_aligns_or_present")) {
            Assert.assertTrue(sam.getReadUnmappedFlag(), "Read should be unmapped but isn't");
        }
        // Two pairs in which only the first read should align
        else if (sam.getReadName().equals("both_reads_present_only_first_aligns") ||
                sam.getReadName().equals("read_2_too_many_gaps")) {
            if (sam.getFirstOfPairFlag()) {
                Assert.assertEquals(sam.getReferenceName(), "chr7", "Read should be mapped but isn't");
            } else {
                Assert.assertTrue(sam.getReadUnmappedFlag(), "Read should not be mapped but is");
            }
        } else {
            throw new Exception("Unexpected read name: " + sam.getReadName());
        }
    }
}
 

@Test
public void testMerger() throws Exception {
    final File output = File.createTempFile("mergeTest", ".sam");
    output.deleteOnExit();

    doMergeAlignment(unmappedBam, Collections.singletonList(alignedBam),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true, fasta, output,
            SamPairUtil.PairOrientation.FR, null, null, null, null, null);

    SamReader result = SamReaderFactory.makeDefault().open(output);
    Assert.assertEquals(result.getFileHeader().getSequenceDictionary().getSequences().size(), 8,
            "Number of sequences did not match");
    SAMProgramRecord pg = result.getFileHeader().getProgramRecords().get(0);
    Assert.assertEquals(pg.getProgramGroupId(), "0");
    Assert.assertEquals(pg.getProgramVersion(), "1.0");
    Assert.assertEquals(pg.getCommandLine(), "align!");
    Assert.assertEquals(pg.getProgramName(), "myAligner");

    final SAMReadGroupRecord rg = result.getFileHeader().getReadGroups().get(0);
    Assert.assertEquals(rg.getReadGroupId(), "0");
    Assert.assertEquals(rg.getSample(), "Hi,Mom!");

    Assert.assertEquals(result.getFileHeader().getSortOrder(), SAMFileHeader.SortOrder.coordinate);

    for (final SAMRecord sam : result) {
        // This tests that we clip both (a) when the adapter is marked in the unmapped BAM file and
        // (b) when the insert size is less than the read length
        if (sam.getReadName().equals("both_reads_align_clip_adapter") ||
                sam.getReadName().equals("both_reads_align_clip_marked")) {
            Assert.assertEquals(sam.getReferenceName(), "chr7");
            if (sam.getReadNegativeStrandFlag()) {
                Assert.assertEquals(sam.getCigarString(), "5S96M", "Incorrect CIGAR string for " +
                        sam.getReadName());
            } else {
                Assert.assertEquals(sam.getCigarString(), "96M5S", "Incorrect CIGAR string for " +
                        sam.getReadName());
            }
        }
        // This tests that we DON'T clip when we run off the end if there are equal to or more than
        // MIN_ADAPTER_BASES hanging off the end
        else if (sam.getReadName().equals("both_reads_align_min_adapter_bases_exceeded")) {
            Assert.assertEquals(sam.getReferenceName(), "chr7");
            Assert.assertTrue(sam.getCigarString().indexOf("S") == -1,
                    "Read was clipped when it should not be.");
        } else if (sam.getReadName().equals("neither_read_aligns_or_present")) {
            Assert.assertTrue(sam.getReadUnmappedFlag(), "Read should be unmapped but isn't");
        }
        // Two pairs in which only the first read should align
        else if (sam.getReadName().equals("both_reads_present_only_first_aligns") ||
                sam.getReadName().equals("read_2_too_many_gaps")) {
            if (sam.getFirstOfPairFlag()) {
                Assert.assertEquals(sam.getReferenceName(), "chr7", "Read should be mapped but isn't");
            } else {
                Assert.assertTrue(sam.getReadUnmappedFlag(), "Read should not be mapped but is");
            }
        } else {
            throw new Exception("Unexpected read name: " + sam.getReadName());
        }

        final boolean pos = !sam.getReadNegativeStrandFlag();
        Assert.assertEquals(sam.getAttribute("aa"), pos ? "Hello" : "olleH");
        Assert.assertEquals(sam.getAttribute("ab"), pos ? "ATTCGG" : "CCGAAT");
        Assert.assertEquals(sam.getAttribute("ac"), pos ? new byte[] {1,2,3} : new byte[] {3,2,1});
        Assert.assertEquals(sam.getAttribute("as"), pos ? new short[]{1,2,3} : new short[]{3,2,1});
        Assert.assertEquals(sam.getAttribute("ai"), pos ? new int[]  {1,2,3} : new int[]  {3,2,1});
        Assert.assertEquals(sam.getAttribute("af"), pos ? new float[]{1,2,3} : new float[]{3,2,1});
    }

    // Quick test to make sure the program record gets picked up from the file if not specified
    // on the command line.
    doMergeAlignment(unmappedBam, Collections.singletonList(alignedBam),
            null, null, null, null,
            false, true, false, 1,
            null, null, null, null,
            true, fasta, output,
            SamPairUtil.PairOrientation.FR, null, null, null, null, null);

    CloserUtil.close(result);

    result = SamReaderFactory.makeDefault().open(output);
    pg = result.getFileHeader().getProgramRecords().get(0);
    Assert.assertEquals(pg.getProgramGroupId(), "1",
            "Program group ID not picked up correctly from aligned BAM");
    Assert.assertEquals(pg.getProgramVersion(), "2.0",
            "Program version not picked up correctly from aligned BAM");
    Assert.assertNull(pg.getCommandLine(),
            "Program command line not picked up correctly from aligned BAM");
    Assert.assertEquals(pg.getProgramName(), "Hey!",
            "Program name not picked up correctly from aligned BAM");

    CloserUtil.close(result);
}
 

@Test
public void testMergerWithSupplemental() throws Exception {
    final File outputWithSupplemental = File.createTempFile("mergeWithSupplementalTest", ".sam");
    outputWithSupplemental.deleteOnExit();

    doMergeAlignment(unmappedBam,
            Collections.singletonList(supplementalReadAlignedBam),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true, fasta, outputWithSupplemental,
            SamPairUtil.PairOrientation.FR, null, null, null, null, null);

    final SamReader result = SamReaderFactory.makeDefault().open(outputWithSupplemental);

    final List<Integer> clipAdapterFlags = new ArrayList<Integer>(Arrays.asList(99, 2147, 147, 2195));
    final List<Integer> foundClipAdapterFlags = new ArrayList<Integer>();

    for (final SAMRecord sam : result) {
        if (sam.getReadName().equals("both_reads_align_clip_adapter")) {
            foundClipAdapterFlags.add(sam.getFlags());
        }

        // This tests that we clip both (a) when the adapter is marked in the unmapped BAM file and
        // (b) when the insert size is less than the read length
        if (sam.getReadName().equals("both_reads_align_clip_marked")) {
            Assert.assertEquals(sam.getReferenceName(), "chr7");
            if (sam.getReadNegativeStrandFlag()) {
                Assert.assertEquals(sam.getCigarString(), "5S96M", "Incorrect CIGAR string for " + sam.getReadName());
            } else {
                Assert.assertEquals(sam.getCigarString(), "96M5S", "Incorrect CIGAR string for " + sam.getReadName());
            }
        }
        else if (sam.getReadName().equals("both_reads_align_clip_adapter")) {
            Assert.assertEquals(sam.getReferenceName(), "chr7");
            if (!sam.getSupplementaryAlignmentFlag()) {
                if (sam.getReadNegativeStrandFlag()) Assert.assertEquals(sam.getCigarString(), "5S96M", "Incorrect CIGAR string for " + sam.getReadName());
                else Assert.assertEquals(sam.getCigarString(), "96M5S", "Incorrect CIGAR string for " + sam.getReadName());
            }
        }
        // This tests that we DON'T clip when we run off the end if there are equal to or more than
        // MIN_ADAPTER_BASES hanging off the end
        else if (sam.getReadName().equals("both_reads_align_min_adapter_bases_exceeded")) {
            Assert.assertEquals(sam.getReferenceName(), "chr7");
            Assert.assertTrue(!sam.getCigarString().contains("S"),
                    "Read was clipped when it should not be.");
        } else if (sam.getReadName().equals("neither_read_aligns_or_present")) {
            Assert.assertTrue(sam.getReadUnmappedFlag(), "Read should be unmapped but isn't");
        }
        // Two pairs in which only the first read should align
        else if (sam.getReadName().equals("both_reads_present_only_first_aligns") ||
                sam.getReadName().equals("read_2_too_many_gaps")) {
            if (sam.getFirstOfPairFlag()) {
                Assert.assertEquals(sam.getReferenceName(), "chr7", "Read should be mapped but isn't");
            } else {
                Assert.assertTrue(sam.getReadUnmappedFlag(), "Read should not be mapped but is");
            }
        } else {
            throw new Exception("Unexpected read name: " + sam.getReadName());
        }
    }

    // Make sure that we have the appropriate primary and supplementary reads in the new file
    Assert.assertEquals(foundClipAdapterFlags.size(), clipAdapterFlags.size());
    Collections.sort(clipAdapterFlags);
    Collections.sort(foundClipAdapterFlags);
    for (int i = 0; i < clipAdapterFlags.size(); i++) {
        Assert.assertEquals(clipAdapterFlags.get(i), foundClipAdapterFlags.get(i));
    }

}
 

private static boolean matchesExpectedOrientations(final SAMRecord rec, final Set<PairOrientation> expectedOrientations) {
    return expectedOrientations.contains(SamPairUtil.getPairOrientation(rec)) && rec.getAttribute(SAMTag.SA.toString()) == null;
}
 

public InsertSizeCollectorArgs(final int insertSize, final SamPairUtil.PairOrientation po) {
    this.insertSize = insertSize;
    this.po = po;
}
 

public SamPairUtil.PairOrientation getPairOrientation() {
    return po;
}
 

public void addMetricsToFile(final MetricsFile<InsertSizeMetrics,Integer> file) {
    // get the number of inserts, and the maximum and minimum keys across, across all orientations
    for (final Histogram<Integer> h : this.histograms.values()) {
        totalInserts += h.getCount();
    }
    if (0 == totalInserts) return; // nothing to store

    for(final Map.Entry<SamPairUtil.PairOrientation, Histogram<Integer>> entry : histograms.entrySet()) {
        final SamPairUtil.PairOrientation pairOrientation = entry.getKey();
        final Histogram<Integer> histogram = entry.getValue();
        final double total = histogram.getCount();

        // Only include a category if it has a sufficient percentage of the data in it
        if( total >= totalInserts * minimumPct ) {
            final InsertSizeMetrics metrics = new InsertSizeMetrics();
            metrics.SAMPLE             = this.sample;
            metrics.LIBRARY            = this.library;
            metrics.READ_GROUP         = this.readGroup;
            metrics.PAIR_ORIENTATION   = pairOrientation;
            if (!histogram.isEmpty()) {
                metrics.READ_PAIRS = (long) total;
                metrics.MAX_INSERT_SIZE = (int) histogram.getMax();
                metrics.MIN_INSERT_SIZE = (int) histogram.getMin();
                metrics.MEDIAN_INSERT_SIZE = histogram.getMedian();
                metrics.MODE_INSERT_SIZE = histogram.getMode();
                metrics.MEDIAN_ABSOLUTE_DEVIATION = histogram.getMedianAbsoluteDeviation();

                final double median = histogram.getMedian();
                double covered = 0;
                double low = median;
                double high = median;

                while (low >= histogram.getMin()-1 || high <= histogram.getMax()+1) {
                    final Histogram.Bin<Integer> lowBin = histogram.get((int) low);
                    if (lowBin != null) covered += lowBin.getValue();

                    if (low != high) {
                        final Histogram.Bin<Integer> highBin = histogram.get((int) high);
                        if (highBin != null) covered += highBin.getValue();
                    }

                    final double percentCovered = covered / total;
                    final int distance = (int) (high - low) + 1;
                    if (percentCovered >= 0.1 && metrics.WIDTH_OF_10_PERCENT == 0) metrics.WIDTH_OF_10_PERCENT = distance;
                    if (percentCovered >= 0.2 && metrics.WIDTH_OF_20_PERCENT == 0) metrics.WIDTH_OF_20_PERCENT = distance;
                    if (percentCovered >= 0.3 && metrics.WIDTH_OF_30_PERCENT == 0) metrics.WIDTH_OF_30_PERCENT = distance;
                    if (percentCovered >= 0.4 && metrics.WIDTH_OF_40_PERCENT == 0) metrics.WIDTH_OF_40_PERCENT = distance;
                    if (percentCovered >= 0.5 && metrics.WIDTH_OF_50_PERCENT == 0) metrics.WIDTH_OF_50_PERCENT = distance;
                    if (percentCovered >= 0.6 && metrics.WIDTH_OF_60_PERCENT == 0) metrics.WIDTH_OF_60_PERCENT = distance;
                    if (percentCovered >= 0.7 && metrics.WIDTH_OF_70_PERCENT == 0) metrics.WIDTH_OF_70_PERCENT = distance;
                    if (percentCovered >= 0.8 && metrics.WIDTH_OF_80_PERCENT == 0) metrics.WIDTH_OF_80_PERCENT = distance;
                    if (percentCovered >= 0.9 && metrics.WIDTH_OF_90_PERCENT == 0) metrics.WIDTH_OF_90_PERCENT = distance;
                    if (percentCovered >= 0.95 && metrics.WIDTH_OF_95_PERCENT == 0) metrics.WIDTH_OF_95_PERCENT = distance;
                    if (percentCovered >= 0.99 && metrics.WIDTH_OF_99_PERCENT == 0) metrics.WIDTH_OF_99_PERCENT = distance;

                    --low;
                    ++high;
                }
            }

            // Trim the Histogram down to get rid of outliers that would make the chart useless.
            final Histogram<Integer> trimmedHistogram = histogram; // alias it
            trimmedHistogram.trimByWidth(getWidthToTrimTo(metrics));

            if (!trimmedHistogram.isEmpty()) {
                metrics.MEAN_INSERT_SIZE = trimmedHistogram.getMean();
                metrics.STANDARD_DEVIATION = trimmedHistogram.getStandardDeviation();
            }

            file.addHistogram(trimmedHistogram);
            file.addMetric(metrics);
        }
    }
}
 

SAMRecord[] toSAMRecord(SAMFileHeader header) {

		SAMRecord first = firstSAMRecord(header);
		SAMRecord second = secondSAMRecord(header);

		SamPairUtil.setMateInfo(first, second, header);

		first.setProperPairFlag(true);
		second.setProperPairFlag(true);

		return new SAMRecord[] { first, second };
	}