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

private boolean passesFilters(@NotNull final SAMRecord record)
{
    if(record.getMappingQuality() < mMinMappingQuality || record.getReadUnmappedFlag())
        return false;

    if(mDropSupplementaries && (record.getSupplementaryAlignmentFlag() || record.isSecondaryAlignment()))
        return false;

    if(mDropDuplicates && record.getDuplicateReadFlag())
        return false;

    return true;
}
 

private HitsForInsert nextMaybeEmpty() {
    if (!peekIterator.hasNext()) throw new IllegalStateException();
    final String readName = peekIterator.peek().getReadName();
    final HitsForInsert hits = new HitsForInsert();

    Boolean isPaired = null;

    // Accumulate the alignments matching readName.
    do {
        final SAMRecord rec = peekIterator.next();
        replaceHardWithSoftClips(rec);
        // It is critical to do this here, because SamAlignmentMerger uses this exception to determine
        // if the aligned input needs to be sorted.
        if (peekIterator.hasNext() && queryNameComparator.fileOrderCompare(rec, peekIterator.peek()) > 0) {
            throw new IllegalStateException("Underlying iterator is not queryname sorted: " +
            rec + " > " + peekIterator.peek());
        }

        if (isPaired == null) {
            isPaired = rec.getReadPairedFlag();
        } else if (isPaired != rec.getReadPairedFlag()) {
            throw new PicardException("Got a mix of paired and unpaired alignments for read " + readName);
        }

        // Records w/ a supplemental flag are stashed to the side until the primary alignment has
        // been determined, and then re-added into the process later
        if (!rec.getReadPairedFlag() || rec.getFirstOfPairFlag()) {
            if (rec.getSupplementaryAlignmentFlag()) {
                hits.addSupplementalFirstOfPairOrFragment(rec);
            } else {
                hits.addFirstOfPairOrFragment(rec);
            }
        } else if (rec.getSecondOfPairFlag()) {
            if (rec.getSupplementaryAlignmentFlag()) {
                hits.addSupplementalSecondOfPair(rec);
            } else {
                hits.addSecondOfPair(rec);
            }
        } else throw new PicardException("Read is marked as pair but neither first or second: " + readName);
    } while (peekIterator.hasNext() && peekIterator.peek().getReadName().equals(readName));

    // If there is no more than one alignment for each end, no need to do any coordination.
    if (hits.numHits() <= 1) {
        // No HI tags needed if only a single hit
        if (hits.getFirstOfPair(0) != null) {
            hits.getFirstOfPair(0).setAttribute(SAMTag.HI.name(), null);
            hits.getFirstOfPair(0).setNotPrimaryAlignmentFlag(false);
        }
        if (hits.getSecondOfPair(0) != null) {
            hits.getSecondOfPair(0).setAttribute(SAMTag.HI.name(), null);
            hits.getSecondOfPair(0).setNotPrimaryAlignmentFlag(false);
        }
    } else {
        primaryAlignmentSelectionStrategy.pickPrimaryAlignment(hits);
    }

    // Used to check that alignments for first and second were correlated, but this is no longer required.
    return hits;
}
 

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

private void collectReadData(final SAMRecord record) {
    // NB: for read count metrics, do not include supplementary records, but for base count metrics, do include supplementary records.
    if (record.getSupplementaryAlignmentFlag()) return;

    metrics.TOTAL_READS++;
    readLengthHistogram.increment(record.getReadBases().length);

    if (!record.getReadFailsVendorQualityCheckFlag()) {
        metrics.PF_READS++;
        if (isNoiseRead(record)) metrics.PF_NOISE_READS++;

        // See if the read is an adapter sequence
        if (adapterUtility.isAdapter(record)) {
            this.adapterReads++;
        }

        if (!record.getReadUnmappedFlag() && doRefMetrics) {
            metrics.PF_READS_ALIGNED++;
            if (record.getReadPairedFlag() && !record.getProperPairFlag()) metrics.PF_READS_IMPROPER_PAIRS++;
            if (!record.getReadNegativeStrandFlag()) numPositiveStrand++;
            if (record.getReadPairedFlag() && !record.getMateUnmappedFlag()) {
                metrics.READS_ALIGNED_IN_PAIRS++;

                // Check that both ends have mapq > minimum
                final Integer mateMq = record.getIntegerAttribute(SAMTag.MQ.toString());
                if (mateMq == null || mateMq >= MAPPING_QUALITY_THRESHOLD && record.getMappingQuality() >= MAPPING_QUALITY_THRESHOLD) {
                    ++this.chimerasDenominator;

                    // With both reads mapped we can see if this pair is chimeric
                    if (ChimeraUtil.isChimeric(record, maxInsertSize, expectedOrientations)) {
                        ++this.chimeras;
                    }
                }
            } else { // fragment reads or read pairs with one end that maps
                // Consider chimeras that occur *within* the read using the SA tag
                if (record.getMappingQuality() >= MAPPING_QUALITY_THRESHOLD) {
                    ++this.chimerasDenominator;
                    if (record.getAttribute(SAMTag.SA.toString()) != null) ++this.chimeras;
                }
            }
        }
    }
}
 

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

}