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

/**
 * Calculates edit distance for the input read.
 * If the input c2r is not null, compare to the actual reference.
 * If c2r is null, check the input NM tag.
 */
public static int getEditDistance(SAMRecord read, CompareToReference2 c2r, boolean includeSoftClipping) {
	
	Integer distance = null;
	
	if (read.getReadUnmappedFlag()) {
		distance = read.getReadLength();
	} else if (c2r != null) {
		//distance = c2r.numMismatches(read) + getNumIndelBases(read);
		distance = c2r.numMismatches(read, includeSoftClipping) + getNumIndelBases(read);
	} else {
		distance = read.getIntegerAttribute("NM");
		
		if (distance == null) {
			// STAR format
			distance = read.getIntegerAttribute("nM");
		}
		
		if (distance == null) {
			distance = read.getReadLength();
		}
	}
			
	return distance;
}
 

/**
 *  Returns original edit distance as set in YX tag.
 */
public static int getOrigEditDistance(SAMRecord read) {
	
	Integer distance = null;
	
	if (read.getReadUnmappedFlag()) {
		distance = read.getReadLength();
	} else {
		distance = read.getIntegerAttribute("YX");
		
		if (distance == null) {
			distance = read.getReadLength();
		}
	}
			
	return distance;
}
 

/**
 * Convenience method for retrieving int attribute
 */
public static int getIntAttribute(SAMRecord read, String attribute) {
	Integer num = read.getIntegerAttribute(attribute);
	
	if (num == null) {
		return 0;
	} else {
		return num;
	}
}
 

/**
 * Get the NH or IH value
 * @param sam the sam record
 * @return the value contained in the NH attribute, or if that doesn't exist,
 * the IH attribute. If neither exist, returns null.
 */
public static Integer getNHOrIH(SAMRecord sam) {
  final Integer nh = sam.getIntegerAttribute(ATTRIBUTE_NH);
  if (nh != null) {
    return nh;
  }
  return sam.getIntegerAttribute(ATTRIBUTE_IH);
}
 

/**
 * Writes out a SAMRecord in Maq fastq format
 *
 * @param codec the code to write to
 * @param rec   the SAMRecord to write
 */
private void writeFastqRecord(final BinaryCodec codec, final SAMRecord rec) {

    // Trim the run barcode off the read name
    String readName = rec.getReadName();
    if (namePrefix != null && readName.startsWith(namePrefix)) {
        readName = readName.substring(nameTrim);
    }
    // Writes the length of the read name and then the name (null-terminated)
    codec.writeString(readName, true, true);

    final char[] seqs = rec.getReadString().toCharArray();
    final char[] quals = rec.getBaseQualityString().toCharArray();

    int retainedLength = seqs.length;
    if (clipAdapters){
        // adjust to a shorter length iff clipping tag exists
        Integer trimPoint = rec.getIntegerAttribute(ReservedTagConstants.XT);
        if (trimPoint != null) {
            ValidationUtils.validateArg(rec.getReadLength() == seqs.length, () -> "length of read and seqs differ. Found " + rec.getReadLength() + " and '" + seqs.length + ".");

            retainedLength = Math.min(seqs.length, Math.max(SEED_REGION_LENGTH, trimPoint -1));
        }
    }

    // Write the length of the sequence
    codec.writeInt(basesToWrite != null ? basesToWrite : seqs.length);

    // Calculate and write the sequence and qualities
    final byte[] seqsAndQuals = encodeSeqsAndQuals(seqs, quals, retainedLength);
    codec.writeBytes(seqsAndQuals);
}
 

public int compare(final SAMRecord rec1, final SAMRecord rec2) {
    final Integer hi1 = rec1.getIntegerAttribute(SAMTag.HI.name());
    final Integer hi2 = rec2.getIntegerAttribute(SAMTag.HI.name());
    if (hi1 == null) {
        if (hi2 == null) return 0;
        else return 1;
    } else if (hi2 == null) {
        return -1;
    } else {
        return hi1.compareTo(hi2);
    }
}
 

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

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

@Override
public int compare(final SAMRecord rec1, final SAMRecord rec2) {
    final Integer hi1 = rec1.getIntegerAttribute(SAMTag.HI.name());
    final Integer hi2 = rec2.getIntegerAttribute(SAMTag.HI.name());
    if (hi1 == null) {
        if (hi2 == null) return 0;
        else return 1;
    } else if (hi2 == null) {
        return -1;
    } else {
        return hi1.compareTo(hi2);
    }
}
 

public boolean qcCheck(SAMRecord sr){
Cigar cigar = sr.getCigar();
int rLen = sr.getReadLength();
int effectiveLen = 0;
if(cigar==null) 
    return false;
else{
    for(final CigarElement ce : cigar.getCigarElements()){
	CigarOperator op = ce.getOperator();
	int cigarLen = ce.getLength();
	switch(op)
	    {
	    case M:
		{
		    effectiveLen += cigarLen;
		    break;
		}
	    case I:
		{
		    effectiveLen += cigarLen;
		    break;
		}
	    default: 
		break;
	    }
    }
}
boolean readdebug = false;
if(readdebug){
    HLA.log.appendln(sr.getSAMString());
    HLA.log.appendln("EffectiveLen:\t" + effectiveLen);
    HLA.log.appendln("ReadLen:\t" + rLen);
}
Integer i = sr.getIntegerAttribute("NM");
int nm = 0;
if(i!=null)
    nm = i.intValue();
if(readdebug)
    HLA.log.appendln("NM=\t" + nm);
if(nm < 16){
    if(readdebug)
	HLA.log.appendln("PASSWED QC");
    return true;
}
if(readdebug){
    HLA.log.appendln("FAILED QC");
    HLA.log.appendln(sr.getSAMString());
}
return false;
   }
 

/**
 * Filter a SAM record based on specified filtering parameters.
 * Does not consider position based filters, or inversion status.
 *
 * @param params parameters
 * @param rec record
 * @return true if record should be accepted for processing
 */
private static boolean checkRecordProperties(final SamFilterParams params, final SAMRecord rec) {
  final int flags = rec.getFlags();
  if ((flags & params.requireUnsetFlags()) != 0) {
    return false;
  }
  if ((flags & params.requireSetFlags()) != params.requireSetFlags()) {
    return false;
  }
  if (rec.getAlignmentStart() == 0 && params.excludeUnplaced()) {
    return false;
  }

  final boolean mated = (flags & SamBamConstants.SAM_READ_IS_MAPPED_IN_PROPER_PAIR) != 0;
  final boolean unmapped = rec.getReadUnmappedFlag();
  if (!mated && !unmapped && params.excludeUnmated()) {
    return false;
  }

  final int minMapQ = params.minMapQ();
  if (minMapQ >= 0 && rec.getMappingQuality() < minMapQ) {
    return false;
  }

  final Integer nh = SamUtils.getNHOrIH(rec);
  final int maxNH = params.maxAlignmentCount();
  if (maxNH >= 0 && nh != null && nh > maxNH) {
      return false;
  }
  if (params.excludeVariantInvalid()) {
    if (nh != null && nh == 0) {
      return false;
    }
    if (!rec.getReadUnmappedFlag() && rec.getAlignmentStart() <= 0) {
      return false;
    }
  }

  final IntegerOrPercentage maxAS = mated ? params.maxMatedAlignmentScore() : params.maxUnmatedAlignmentScore();
  if (maxAS != null) {
    final Integer as = rec.getIntegerAttribute(SamUtils.ATTRIBUTE_ALIGNMENT_SCORE);
    if (as != null && as > maxAS.getValue(rec.getReadLength())) {
      return false;
    }
  }
  final int minReadLength = params.minReadLength();
  if (minReadLength >= 0 && rec.getReadLength() < minReadLength) {
    return false;
  }

  if (params.subsampleFraction() != null) {
    final double sFrac;
    if (params.subsampleRampFraction() != null) { // Use subsample ramping
      if (rec.getAlignmentStart() == 0) {
        sFrac = (params.subsampleFraction() + params.subsampleRampFraction()) / 2;
      } else {
        final int pos = rec.getAlignmentStart();
        final int refLength = rec.getHeader().getSequence(rec.getReferenceIndex()).getSequenceLength();
        final double lengthFrac = Math.max(0, Math.min(1.0, (double) pos / refLength));
        sFrac = params.subsampleFraction() + lengthFrac * (params.subsampleRampFraction() - params.subsampleFraction());
      }
    } else {
      sFrac = params.subsampleFraction();
    }
    // Subsample using hash of read name, ensures pairs are kept together
    return (internalHash(rec.getReadName(), params.subsampleSeed()) & SUBSAMPLE_MASK) < sFrac * SUBSAMPLE_MAX;
  }

  if (params.selectReadGroups() != null) {
    final SAMReadGroupRecord srg = rec.getReadGroup();
    if (srg == null || !params.selectReadGroups().contains(srg.getReadGroupId())) {
      return false;
    }
  }

  return true;
}
 

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

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

/**
 * Run IlluminaBasecallsToSam on a few test cases, and verify that results agree with hand-checked expectation.
 */
@Test(dataProvider="data")
public void testBasic(final String LANE, final String readStructure,
                      final String fivePrimerAdapter, final String threePrimerAdapter,
                      final int expectedNumClippedRecords) throws Exception {
    // Create the SAM file from Gerald output
    final File samFile = File.createTempFile("." + LANE + ".illuminaBasecallsToSam", ".sam");
    samFile.deleteOnExit();
    final List<String> illuminaArgv = CollectionUtil.makeList("BASECALLS_DIR=" + TEST_DATA_DIR,
            "LANE=" + LANE,
            "RUN_BARCODE=" + RUN_BARCODE,
            "READ_STRUCTURE=" + readStructure,
            "OUTPUT=" + samFile,
            "SEQUENCING_CENTER=BI",
            "ALIAS=" + ALIAS
    );
    if (fivePrimerAdapter != null) {
        illuminaArgv.addAll(CollectionUtil.makeList(
            "ADAPTERS_TO_CHECK=null",
            "FIVE_PRIME_ADAPTER=" + fivePrimerAdapter,
            "THREE_PRIME_ADAPTER=" + threePrimerAdapter
        ));
    }
    Assert.assertEquals(runPicardCommandLine(illuminaArgv), 0);

    // Read the file and confirm it contains what is expected
    final SamReader samReader = SamReaderFactory.makeDefault().open(samFile);

    // look for clipped adaptor attribute in lane 3 PE (2) and in lane 6 (1) non-PE
    int count = 0;
    for (final SAMRecord record : samReader) {
        if (record.getIntegerAttribute(ReservedTagConstants.XT) != null) {
            count++;
            if ((count == 1 || count == 2) && LANE.equals("2")){
                Assert.assertEquals (114, (int)record.getIntegerAttribute(ReservedTagConstants.XT));
            } else if (count == 1 || count == 2 && LANE.equals("1")) {
                Assert.assertEquals(68, (int) record.getIntegerAttribute(ReservedTagConstants.XT));
            }
        }
    }

    // Check the total number of clipped records
    Assert.assertEquals(count, expectedNumClippedRecords);

    samReader.close();
    samFile.delete();
}
 

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