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

public void addRecord(final SAMRecord record, final ReferenceSequence ref) {
    if (record.getNotPrimaryAlignmentFlag()) {
        // only want 1 count per read so skip non primary alignments
        return;
    }

    collectReadData(record);
    collectQualityData(record, ref);
}
 

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

public void acceptRecord(final SAMRecord rec, final ReferenceSequence ref) {
    if (!this.includeSecondaryAlignments    && rec.getNotPrimaryAlignmentFlag()) return;
    if (!this.includeSupplementalAlignments && rec.getSupplementaryAlignmentFlag()) return;

    final int length = rec.getReadLength();
    metrics.TOTAL_READS++;
    metrics.TOTAL_BASES += length;

    final boolean isPfRead = !rec.getReadFailsVendorQualityCheckFlag();
    if (isPfRead) {
        metrics.PF_READS++;
        metrics.PF_BASES += length;
    }

    final byte[] quals;
    if (this.useOriginalQualities) {
        byte[] tmp = rec.getOriginalBaseQualities();
        if (tmp == null) tmp = rec.getBaseQualities();
        quals = tmp;
    } else {
        quals = rec.getBaseQualities();
    }

    // add up quals, and quals >= 20
    for (final int qual : quals) {
        metrics.Q20_EQUIVALENT_YIELD += qual;

        if (qual >= 30) {
            metrics.Q20_BASES++;
            metrics.Q30_BASES++;
        } else if (qual >= 20) {
            metrics.Q20_BASES++;
        }

        if (isPfRead) {
            metrics.PF_Q20_EQUIVALENT_YIELD += qual;
            if (qual >= 30) {
                metrics.PF_Q20_BASES++;
                metrics.PF_Q30_BASES++;
            } else if (qual >= 20) {
                metrics.PF_Q20_BASES++;
            }
        }
    }
}
 

@Test
public void testMultiHit() throws IOException {
    final File unmappedSam = new File(TEST_DATA_DIR, "multihit.unmapped.sam");
    final File alignedSam = new File(TEST_DATA_DIR, "multihit.aligned.sam");
    final File merged = File.createTempFile("merged", ".sam");
    merged.deleteOnExit();

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

    // Iterate over the merged output and gather some statistics
    final Map<String, AlignmentAccumulator> accumulatorMap = new HashMap<String, AlignmentAccumulator>();

    final SamReader reader = SamReaderFactory.makeDefault().open(merged);
    for (final SAMRecord rec : reader) {
        final String readName;
        if (!rec.getReadPairedFlag()) readName = rec.getReadName();
        else if (rec.getFirstOfPairFlag()) readName = rec.getReadName() + "/1";
        else readName = rec.getReadName() + "/2";
        AlignmentAccumulator acc = accumulatorMap.get(readName);
        if (acc == null) {
            acc = new AlignmentAccumulator();
            accumulatorMap.put(readName, acc);
        }
        if (rec.getReadUnmappedFlag()) {
            Assert.assertFalse(acc.seenUnaligned, readName);
            acc.seenUnaligned = true;
        } else {
            ++acc.numAlignments;
            if (!rec.getNotPrimaryAlignmentFlag()) {
                Assert.assertNull(acc.primaryAlignmentSequence, readName);
                acc.primaryAlignmentSequence = rec.getReferenceName();
            }
        }
    }

    // Set up expected results
    final Map<String, AlignmentAccumulator> expectedMap = new HashMap<String, AlignmentAccumulator>();
    expectedMap.put("frag_hit", new AlignmentAccumulator(false, 1, "chr1"));
    expectedMap.put("frag_multihit", new AlignmentAccumulator(false, 3, "chr2"));
    expectedMap.put("frag_no_hit", new AlignmentAccumulator(true, 0, null));
    expectedMap.put("pair_both_hit/1", new AlignmentAccumulator(false, 1, "chr7"));
    expectedMap.put("pair_both_hit/2", new AlignmentAccumulator(false, 1, "chr7"));
    expectedMap.put("pair_both_multihit/1", new AlignmentAccumulator(false, 3, "chr8"));
    expectedMap.put("pair_both_multihit/2", new AlignmentAccumulator(false, 3, "chr8"));
    expectedMap.put("pair_first_hit/1", new AlignmentAccumulator(false, 1, "chr1"));
    expectedMap.put("pair_first_hit/2", new AlignmentAccumulator(true, 0, null));
    expectedMap.put("pair_first_multihit/1", new AlignmentAccumulator(false, 3, "chr1"));
    expectedMap.put("pair_first_multihit/2", new AlignmentAccumulator(true, 0, null));
    expectedMap.put("pair_no_hit/1", new AlignmentAccumulator(true, 0, null));
    expectedMap.put("pair_no_hit/2", new AlignmentAccumulator(true, 0, null));
    expectedMap.put("pair_second_hit/1", new AlignmentAccumulator(true, 0, null));
    expectedMap.put("pair_second_hit/2", new AlignmentAccumulator(false, 1, "chr1"));
    expectedMap.put("pair_second_multihit/1", new AlignmentAccumulator(true, 0, null));
    expectedMap.put("pair_second_multihit/2", new AlignmentAccumulator(false, 3, "chr3"));

    Assert.assertEquals(accumulatorMap.size(), expectedMap.size());

    for (final Map.Entry<String, AlignmentAccumulator> entry : expectedMap.entrySet()) {
        final AlignmentAccumulator actual = accumulatorMap.get(entry.getKey());
        Assert.assertEquals(actual, entry.getValue(), entry.getKey());
    }

    assertSamValid(merged);
}
 

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

/**
 * Read a single fragment read from a file, and create one or more aligned records for the read pair based on
 * the lists, merge with the original read, and assert expected results.
 * @param description
 * @param hitSpecs List that describes the aligned SAMRecords to create.
 * @param expectedPrimaryHitIndex Expected value for the HI tag in the primary alignment in the merged output.
 * @param expectedNumReads Expected number of SAMRecords in the merged output.
 * @param expectedPrimaryMapq Mapq of both ends of primary alignment in the merged output.
 * @throws Exception
 */
@Test(dataProvider = "testFragmentMultiHitWithFilteringTestCases")
public void testFragmentMultiHitWithFiltering(final String description, final List<HitSpec> hitSpecs,
                                            final Integer expectedPrimaryHitIndex, final int expectedNumReads,
                                            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.fragment.unmapped.sam");
    final SAMRecordIterator unmappedSamFileIterator = SamReaderFactory.makeDefault().open(unmappedSam).iterator();
    final SAMRecord unmappedRec = 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 < hitSpecs.size(); ++i) {
        final HitSpec hitSpec = hitSpecs.get(i);
        final SAMRecord mappedRec = makeRead(alignedHeader, unmappedRec, hitSpec, i);
        if (mappedRec != null) {
            alignedWriter.addAlignment(mappedRec);
        }
    }
    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",
            false, 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 numReads = 0;
    Integer primaryHitIndex = null;
    int primaryMapq = 0;
    for (final SAMRecord rec : mergedReader) {
        ++numReads;
        if (!rec.getNotPrimaryAlignmentFlag()  && !rec.getReadUnmappedFlag()) {
            final Integer hitIndex = rec.getIntegerAttribute(SAMTag.HI.name());
            final int newHitIndex = (hitIndex == null? -1: hitIndex);
            Assert.assertNull(primaryHitIndex);
            primaryHitIndex = newHitIndex;
            primaryMapq = rec.getMappingQuality();
        }
    }
    Assert.assertEquals(numReads, expectedNumReads);
    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));
    }
}
 

/**
 * Returns true if the input read is primary.
 * i.e. Bit flag not secondary 0x200 or supplemental 0x800
 */
public static boolean isPrimary(SAMRecord read) {
	return ((read.getFlags() & 0x800)  == 0) && (!read.getNotPrimaryAlignmentFlag());
}