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

/**
 * Test that the records are properly paired.
 * @param r1 Should be the first read of the pair
 * @param r2 Should be the second read of the pair
 */
public boolean testPairedRead (SAMRecord r1, SAMRecord r2) {
	if (r1==null || r2==null) {
		return (false);
	}
	
	if (r1.getSecondOfPairFlag() || r2.getFirstOfPairFlag()) {
		log.warn(r1.getReadName()+ " passed as 2nd read of pair, should be the first");
		return false;
	}
	if (r2.getFirstOfPairFlag()) {
		log.warn(r2.getReadName()+ " passed as 1st read of pair, should be the second");
		return false;
	}
	boolean result=testSameReadName(r1, r2);
	if (result==false) {
		log.warn("Read names not the same "+ r1.getReadName()+ " "+ r2.getReadName());
	}
	return result;
}
 

/**
 * This is supposed to check if the mates have valid pairing flags.
 * 
 * @param r1
 * @param r2
 * @return
 */
private boolean checkMateFlags(SAMRecord r1, SAMRecord r2) {
	if (!r1.getReadPairedFlag() || !r2.getReadPairedFlag())
		return false;

	if (r1.getReadUnmappedFlag() != r2.getMateUnmappedFlag())
		return false;
	if (r1.getReadNegativeStrandFlag() != r2.getMateNegativeStrandFlag())
		return false;
	if (r1.getProperPairFlag() != r2.getProperPairFlag())
		return false;
	if (r1.getFirstOfPairFlag() && r2.getFirstOfPairFlag())
		return false;
	if (r1.getSecondOfPairFlag() && r2.getSecondOfPairFlag())
		return false;

	if (r2.getReadUnmappedFlag() != r1.getMateUnmappedFlag())
		return false;
	if (r2.getReadNegativeStrandFlag() != r1.getMateNegativeStrandFlag())
		return false;

	return true;
}
 

void addRecord(final SAMRecord rec) {
    final byte[] quals = (useOriginalQualities ? rec.getOriginalBaseQualities() : rec.getBaseQualities());
    if (quals == null) return;

    final int length = quals.length;
    final boolean rc = rec.getReadNegativeStrandFlag();
    ensureArraysBigEnough(length+1);

    for (int i=0; i<length; ++i) {
        final int cycle = rc ? length-i : i+1;

        if (rec.getReadPairedFlag() && rec.getSecondOfPairFlag()) {
            secondReadTotalsByCycle[cycle] += quals[i];
            secondReadCountsByCycle[cycle] += 1;
        }
        else {
            firstReadTotalsByCycle[cycle] += quals[i];
            firstReadCountsByCycle[cycle] += 1;
        }
    }
}
 

private boolean areReadsMates(final SAMRecord read1, final SAMRecord read2) {
    // must have same name
    return (read1.getReadName().equals(read2.getReadName()) &&
            // must be paired
            read1.getReadPairedFlag() &&
            // one must be first while the other is not
            read1.getFirstOfPairFlag() != read2.getFirstOfPairFlag() &&
            // one must be second while the other is not
            read1.getSecondOfPairFlag() != read2.getSecondOfPairFlag() &&
            // read1 must be mapped
            !read1.getReadUnmappedFlag() &&
            // read2 must be mapped
            !read2.getReadUnmappedFlag() &&
            // read1 must be non-secondary
            !read1.isSecondaryAlignment() &&
            // read2 must be non-secondary
            !read2.isSecondaryAlignment() &&
            // read1 mate position must agree with read2's position
            read1.getMateAlignmentStart() == read2.getAlignmentStart() &&
            // read1 mate reference must agree with read2's reference
            read1.getMateReferenceIndex().equals(read2.getReferenceIndex())
    );
}
 

private List<SAMRecord> getPairedRead (final String name, final int contig, final int start1, final int start2) {
	List<SAMRecord> result = new ArrayList<> ();

	SAMRecordSetBuilder builder = new SAMRecordSetBuilder();
	builder.addPair(name, contig, start1, start2);

	Collection<SAMRecord> recs = builder.getRecords();

	for (SAMRecord r: recs) {
		if (r.getFirstOfPairFlag()) result.add(0, r);
		if (r.getSecondOfPairFlag()) result.add(1, r);
		r.setMappingQuality(10);
	}

	return (result);

}
 

/**
 * @return a paired read, first of pair in the first position of the list, 2nd of pair in the 2nd position.
 */
private List<SAMRecord> getPairedRead () {
	List<SAMRecord> result = new ArrayList<> ();

	SAMRecordSetBuilder builder = new SAMRecordSetBuilder();
	builder.addUnmappedPair("test");
	Collection<SAMRecord> recs = builder.getRecords();

	for (SAMRecord r: recs) {
		if (r.getFirstOfPairFlag()) result.add(0, r);
		if (r.getSecondOfPairFlag()) result.add(1, r);
	}
	return (result);

}
 

public PrimaryAlignmentKey(final SAMRecord rec) {
    if (rec.isSecondaryOrSupplementary()) {
        throw new IllegalArgumentException("PrimaryAligmentKey cannot be a secondary or suplementary alignment");
    }
    this.pairStatus = rec.getReadPairedFlag() ?
            (rec.getSecondOfPairFlag() ? PairStatus.SECOND : PairStatus.FIRST) :
            PairStatus.UNPAIRED;
    this.readName = rec.getReadName();
}
 

private void countRecord(final SAMRecord rec) {
    final boolean isNegativeStrand = rec.getReadNegativeStrandFlag();
    final boolean isReadTwo = rec.getReadPairedFlag() && rec.getSecondOfPairFlag();
    if (isReadTwo) {
        if (isNegativeStrand) this.R2_NEG++;
        else this.R2_POS++;
    } else {
        if (isNegativeStrand) this.R1_NEG++;
        else this.R1_POS++;
    }
}
 

/**
 * Reads share the same name, and one is flagged first of pair, and the other 2nd of pair.
 * Order of reads is unimportant.
 * @param r1
 * @param r2
 * @return
 */
public boolean testReadsArePaired (SAMRecord r1, SAMRecord r2) {
	// test if the pair of reads are flagged first and second read.
	boolean f1=r1.getFirstOfPairFlag() & r2.getSecondOfPairFlag();
	boolean f2=r1.getSecondOfPairFlag() & r2.getFirstOfPairFlag();
	if (f1 | f2) {
		return (testSameReadName(r1,r2));
	}
	return false;
}
 

void add(final SAMRecord record) {
    if (!record.getReadPairedFlag()) throw new PicardException("Record "+record.getReadName()+" is not paired");
    if (record.getFirstOfPairFlag()) { 
        if (mate1 != null) throw new PicardException("Mate 1 already set for record: "+record.getReadName());
        mate1 = record ;
    }
    else if (record.getSecondOfPairFlag()) { 
        if (mate2 != null) throw new PicardException("Mate 2 already set for record: "+record.getReadName());
        mate2 = record ;
    }
    else throw new PicardException("Neither FirstOfPairFlag or SecondOfPairFlag is set for a paired record");
}
 

public int getStreamIndex(SAMRecord record) {
	if (!record.getReadPairedFlag())
		return 0;
	if (record.getFirstOfPairFlag())
		return 1;
	if (record.getSecondOfPairFlag())
		return 2;

	return 0;
}
 

/**
public ReadPair (SAMRecord read1) {
	this.read1=read1;
}
*/

public ReadPair (final SAMRecord read1, final SAMRecord read2) {
	if (read1.getFirstOfPairFlag())
		this.read1 = read1;
	if (read1.getSecondOfPairFlag())
		this.read2=read1;
	if (read2.getFirstOfPairFlag())
		this.read1=read2;
	if (read2.getSecondOfPairFlag())
		this.read2=read2;
}
 

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

    final SamReader in = SamReaderFactory.makeDefault().referenceSequence(REFERENCE_SEQUENCE).open(INPUT);
    final SAMFileHeader.SortOrder order = in.getFileHeader().getSortOrder();
    SAMFileWriter out = null;
    if (OUTPUT != null) {
        IOUtil.assertFileIsWritable(OUTPUT);
        out = new SAMFileWriterFactory().makeSAMOrBAMWriter(in.getFileHeader(), true, OUTPUT);
    }

    final Histogram<Integer> histo = new Histogram<Integer>("clipped_bases", "read_count");

    // Combine any adapters and custom adapter pairs from the command line into an array for use in clipping
    final AdapterPair[] adapters;
    {
        final List<AdapterPair> tmp = new ArrayList<AdapterPair>();
        tmp.addAll(ADAPTERS);
        if (FIVE_PRIME_ADAPTER != null && THREE_PRIME_ADAPTER != null) {
            tmp.add(new CustomAdapterPair(FIVE_PRIME_ADAPTER, THREE_PRIME_ADAPTER));
        }
        adapters = tmp.toArray(new AdapterPair[tmp.size()]);
    }

    ////////////////////////////////////////////////////////////////////////
    // Main loop that consumes reads, clips them and writes them to the output
    ////////////////////////////////////////////////////////////////////////
    final ProgressLogger progress = new ProgressLogger(log, 1000000, "Read");
    final SAMRecordIterator iterator = in.iterator();

    final AdapterMarker adapterMarker = new AdapterMarker(ADAPTER_TRUNCATION_LENGTH, adapters).
            setMaxPairErrorRate(MAX_ERROR_RATE_PE).setMinPairMatchBases(MIN_MATCH_BASES_PE).
            setMaxSingleEndErrorRate(MAX_ERROR_RATE_SE).setMinSingleEndMatchBases(MIN_MATCH_BASES_SE).
            setNumAdaptersToKeep(NUM_ADAPTERS_TO_KEEP).
            setThresholdForSelectingAdaptersToKeep(PRUNE_ADAPTER_LIST_AFTER_THIS_MANY_ADAPTERS_SEEN);

    while (iterator.hasNext()) {
        final SAMRecord rec = iterator.next();
        final SAMRecord rec2 = rec.getReadPairedFlag() && iterator.hasNext() ? iterator.next() : null;
        rec.setAttribute(ReservedTagConstants.XT, null);

        // Do the clipping one way for PE and another for SE reads
        if (rec.getReadPairedFlag()) {
            // Assert that the input file is in query name order only if we see some PE reads
            if (order != SAMFileHeader.SortOrder.queryname) {
                throw new PicardException("Input BAM file must be sorted by queryname");
            }

            if (rec2 == null) throw new PicardException("Missing mate pair for paired read: " + rec.getReadName());
            rec2.setAttribute(ReservedTagConstants.XT, null);

            // Assert that we did in fact just get two mate pairs
            if (!rec.getReadName().equals(rec2.getReadName())) {
                throw new PicardException("Adjacent reads expected to be mate-pairs have different names: " +
                        rec.getReadName() + ", " + rec2.getReadName());
            }

            // establish which of pair is first and which second
            final SAMRecord first, second;

            if (rec.getFirstOfPairFlag() && rec2.getSecondOfPairFlag()) {
                first = rec;
                second = rec2;
            } else if (rec.getSecondOfPairFlag() && rec2.getFirstOfPairFlag()) {
                first = rec2;
                second = rec;
            } else {
                throw new PicardException("Two reads with same name but not correctly marked as 1st/2nd of pair: " + rec.getReadName());
            }

            adapterMarker.adapterTrimIlluminaPairedReads(first, second);
        } else {
            adapterMarker.adapterTrimIlluminaSingleRead(rec);
        }

        // Then output the records, update progress and metrics
        for (final SAMRecord r : new SAMRecord[]{rec, rec2}) {
            if (r != null) {
                progress.record(r);
                if (out != null) out.addAlignment(r);

                final Integer clip = r.getIntegerAttribute(ReservedTagConstants.XT);
                if (clip != null) histo.increment(r.getReadLength() - clip + 1);
            }
        }
    }

    if (out != null) out.close();

    // Lastly output the metrics to file
    final MetricsFile<?, Integer> metricsFile = getMetricsFile();
    metricsFile.setHistogram(histo);
    metricsFile.write(METRICS);

    CloserUtil.close(in);
    return 0;
}
 

protected static void assertPairedMates(final SAMRecord record1, final SAMRecord record2) {
    if (!(record1.getFirstOfPairFlag() && record2.getSecondOfPairFlag() ||
            record2.getFirstOfPairFlag() && record1.getSecondOfPairFlag())) {
        throw new PicardException("Illegal mate state: " + record1.getReadName());
    }
}
 

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

/**
 *
 */
private Counts computeAlleleFraction(final SamLocusIterator.LocusInfo info, final byte refBase) {
    final Counts counts = new Counts();
    final byte altBase = (refBase == 'C') ? (byte) 'A' : (byte) 'T';

    for (final SamLocusIterator.RecordAndOffset rec : info.getRecordAndOffsets()) {
        final byte qual;
        final SAMRecord samrec = rec.getRecord();

        if (USE_OQ) {
            final byte[] oqs = samrec.getOriginalBaseQualities();
            if (oqs != null) qual = oqs[rec.getOffset()];
            else qual = rec.getBaseQuality();
        } else {
            qual = rec.getBaseQuality();
        }

        // Skip if below qual, or if library isn't a match
        if (qual < MINIMUM_QUALITY_SCORE) continue;
        if (!this.library.equals(Optional.ofNullable(samrec.getReadGroup().getLibrary()).orElse(UNKNOWN_LIBRARY))) continue;

        // Get the read base, and get it in "as read" orientation
        final byte base = rec.getReadBase();
        final byte baseAsRead = samrec.getReadNegativeStrandFlag() ? SequenceUtil.complement(base) : base;
        final int read = samrec.getReadPairedFlag() && samrec.getSecondOfPairFlag() ? 2 : 1;

        // Figure out how to count the alternative allele. If the damage is caused by oxidation of G
        // during shearing (in non-rnaseq data), then we know that:
        //     G>T observation is always in read 1
        //     C>A observation is always in read 2
        // But if the substitution is from other causes the distribution of A/T across R1/R2 will be
        // random.
        if (base == refBase) {
            if (baseAsRead == 'G' && read == 1) ++counts.oxidatedC;
            else if (baseAsRead == 'G' && read == 2) ++counts.controlC;
            else if (baseAsRead == 'C' && read == 1) ++counts.controlC;
            else if (baseAsRead == 'C' && read == 2) ++counts.oxidatedC;
        } else if (base == altBase) {
            if (baseAsRead == 'T' && read == 1) ++counts.oxidatedA;
            else if (baseAsRead == 'T' && read == 2) ++counts.controlA;
            else if (baseAsRead == 'A' && read == 1) ++counts.controlA;
            else if (baseAsRead == 'A' && read == 2) ++counts.oxidatedA;
        }
    }

    return counts;
}
 

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

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

public PrimaryAlignmentKey(final SAMRecord rec) {
    this.pairStatus = rec.getReadPairedFlag() ?
            (rec.getSecondOfPairFlag() ? PairStatus.SECOND : PairStatus.FIRST) :
            PairStatus.UNPAIRED;
    this.readName = rec.getReadName();
}