htsjdk.samtools.SAMRecordSetBuilder Java Examples

@Test
public void testFilter() {
	String tagName = "XC";
	Collection<String> cellBarcodesToKeep = Arrays.asList("CELL1", "CELL2");
	Iterator<SAMRecord> underlyingIterator = Collections.emptyIterator();
	BAMTagValueFilter f = new BAMTagValueFilter(underlyingIterator, tagName, cellBarcodesToKeep);

	SAMRecordSetBuilder b = new SAMRecordSetBuilder();
	// second argument is the contig index which is 0 based.  So contig index=0 -> chr1.  index=2 -> chr3, etc.
	b.addFrag("1", 0, 1, false);
	SAMRecord r1 = b.getRecords().iterator().next();

	r1.setAttribute(tagName, "CELL1");
	Assert.assertFalse(f.filterOut(r1));

	r1.setAttribute(tagName, "CELL2");
	Assert.assertFalse(f.filterOut(r1));

	r1.setAttribute(tagName, "FOO");
	Assert.assertTrue(f.filterOut(r1));

}
 

public File build (final List<SAMRecordSetBuilder> setBuilder, final File unsortedSam, final SAMFileHeader header) throws IOException {
    final File sortedSam = VcfTestUtils.createTemporaryIndexedFile("CollectGcBias", ".bam");

    final SAMFileWriter writer = new SAMFileWriterFactory()
            .setCreateIndex(true).makeBAMWriter(header, false, unsortedSam);

    for (final SAMRecordSetBuilder subSetBuilder : setBuilder) {
        for (final SAMRecord record : subSetBuilder) {
            writer.addAlignment(record);
        }
    }
    writer.close();

    final SortSam sorter = new SortSam();
    final String[] args = new String[] {
            "INPUT=" + unsortedSam.getAbsolutePath(),
            "OUTPUT=" + sortedSam.getAbsolutePath(),
            "SORT_ORDER=coordinate"
    };

    sorter.instanceMain(args);

    return sortedSam;
}
 

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

}
 

@Test
public void testBasicPairedRead() {
    SAMRecordSetBuilder builder = new SAMRecordSetBuilder(false, SAMFileHeader.SortOrder.queryname);
    builder.setReadLength(READ_LENGTH);
    builder.addPair("mapped_paired", 1, 1, 31);
    PeekableIterator<SAMRecord> iterator = new PeekableIterator<SAMRecord>(builder.iterator());

    QuerySortedReadPairIteratorUtil.ReadPair pair = QuerySortedReadPairIteratorUtil.getNextReadPair(iterator);
    Assert.assertNotNull(pair);
    Assert.assertNotNull(pair.read1);
    Assert.assertNotNull(pair.read2);
    Assert.assertEquals("mapped_paired", pair.read1.getReadName());
    Assert.assertEquals("mapped_paired", pair.read2.getReadName());

    pair = QuerySortedReadPairIteratorUtil.getNextReadPair(iterator);
    Assert.assertNull(pair);
}
 

@Test
public void testBasicUnmappedReadPair() {
    SAMRecordSetBuilder builder = new SAMRecordSetBuilder(false, SAMFileHeader.SortOrder.queryname);
    builder.setReadLength(READ_LENGTH);
    builder.addUnmappedPair("unmapped_paired");
    PeekableIterator<SAMRecord> iterator = new PeekableIterator<SAMRecord>(builder.iterator());

    QuerySortedReadPairIteratorUtil.ReadPair pair = QuerySortedReadPairIteratorUtil.getNextReadPair(iterator);
    Assert.assertNotNull(pair);
    Assert.assertNotNull(pair.read1);
    Assert.assertNotNull(pair.read2);
    Assert.assertEquals("unmapped_paired", pair.read1.getReadName());
    Assert.assertEquals("unmapped_paired", pair.read2.getReadName());

    pair = QuerySortedReadPairIteratorUtil.getNextReadPair(iterator);
    Assert.assertNull(pair);
}
 

@Test
public void filterOut3() {
	Iterator<SAMRecord> underlyingIterator = Collections.emptyIterator();
	Collection<String> cellBarcodesToKeep = Collections.EMPTY_LIST;
	String cellBCTag ="XC";

	CellBarcodeFilteringIterator f = new CellBarcodeFilteringIterator(underlyingIterator, cellBCTag, cellBarcodesToKeep);
	SAMRecordSetBuilder b = new SAMRecordSetBuilder();
	// second argument is the contig index which is 0 based.  So contig index=0 -> chr1.  index=2 -> chr3, etc.
	b.addFrag("1", 0, 1, false);
	SAMRecord r1 = b.getRecords().iterator().next();

	r1.setAttribute(cellBCTag, "CELL1");
	Assert.assertFalse(f.filterOut(r1));

	r1.setAttribute(cellBCTag, "CELL2");
	Assert.assertFalse(f.filterOut(r1));

}
 

@Test
public void filterOut2() {
	Iterator<SAMRecord> underlyingIterator = Collections.emptyIterator();
	Collection<String> cellBarcodesToKeep = null;
	String cellBCTag ="XC";

	CellBarcodeFilteringIterator f = new CellBarcodeFilteringIterator(underlyingIterator, cellBCTag, cellBarcodesToKeep);
	SAMRecordSetBuilder b = new SAMRecordSetBuilder();
	// second argument is the contig index which is 0 based.  So contig index=0 -> chr1.  index=2 -> chr3, etc.
	b.addFrag("1", 0, 1, false);
	SAMRecord r1 = b.getRecords().iterator().next();

	r1.setAttribute(cellBCTag, "CELL1");
	Assert.assertFalse(f.filterOut(r1));

	r1.setAttribute(cellBCTag, "CELL2");
	Assert.assertFalse(f.filterOut(r1));

}
 

@Test
public void filterOut() {

	Iterator<SAMRecord> underlyingIterator = Collections.emptyIterator();
	Collection<String> cellBarcodesToKeep = Arrays.asList("CELL1", "CELL2");
	String cellBCTag ="XC";

	CellBarcodeFilteringIterator f = new CellBarcodeFilteringIterator(underlyingIterator, cellBCTag, cellBarcodesToKeep);
	SAMRecordSetBuilder b = new SAMRecordSetBuilder();
	// second argument is the contig index which is 0 based.  So contig index=0 -> chr1.  index=2 -> chr3, etc.
	b.addFrag("1", 0, 1, false);
	SAMRecord r1 = b.getRecords().iterator().next();

	r1.setAttribute(cellBCTag, "CELL1");
	Assert.assertFalse(f.filterOut(r1));

	r1.setAttribute(cellBCTag, "CELL2");
	Assert.assertFalse(f.filterOut(r1));

	r1.setAttribute(cellBCTag, "FOO");
	Assert.assertTrue(f.filterOut(r1));

}
 

@Test
public void testFilterEmpty() {
	Iterator<SAMRecord> underlyingIterator = Collections.emptyIterator();
	Collection<String> contigsToFilter = null;

	ChromosomeFilteringIterator f = new ChromosomeFilteringIterator(underlyingIterator, contigsToFilter, true);

	SAMRecordSetBuilder b = new SAMRecordSetBuilder();
	// second argument is the contig index which is 0 based.  So contig index=0 -> chr1.  index=2 -> chr3, etc.
	b.addFrag("1", 0, 1, false);
	b.addFrag("2", 1, 1, false);
	b.addFrag("3", 2, 1, false);
	b.addFrag("4", 3, 1, false);
	Iterator<SAMRecord> recs = b.getRecords().iterator();

	Assert.assertFalse(f.filterOut(recs.next()));
	Assert.assertFalse(f.filterOut(recs.next()));
	// 3 and 4 not on the list.
	Assert.assertFalse(f.filterOut(recs.next()));
	Assert.assertFalse(f.filterOut(recs.next()));


}
 

@Test
public void testFilterInclude() {
	Iterator<SAMRecord> underlyingIterator = Collections.emptyIterator();
	Collection<String> contigsToFilter = Arrays.asList("chr1", "chr2");

	ChromosomeFilteringIterator f = new ChromosomeFilteringIterator(underlyingIterator, contigsToFilter, false);

	SAMRecordSetBuilder b = new SAMRecordSetBuilder();
	// second argument is the contig index which is 0 based.  So contig index=0 -> chr1.  index=2 -> chr3, etc.
	b.addFrag("chr1", 0, 1, false);
	b.addFrag("chr2", 1, 1, false);
	b.addFrag("chr3", 2, 1, false);
	b.addFrag("chr4", 3, 1, false);
	Iterator<SAMRecord> recs = b.getRecords().iterator();

	Assert.assertFalse(f.filterOut(recs.next()));
	Assert.assertFalse(f.filterOut(recs.next()));
	// 3 and 4 not on the list.  Reject.
	Assert.assertTrue(f.filterOut(recs.next()));
	Assert.assertTrue(f.filterOut(recs.next()));


}
 

@Test
public void testFilterExclude() {
	Iterator<SAMRecord> underlyingIterator = Collections.emptyIterator();
	Collection<String> contigsToFilter = Arrays.asList("chr1", "chr2");

	ChromosomeFilteringIterator f = new ChromosomeFilteringIterator(underlyingIterator, contigsToFilter);

	SAMRecordSetBuilder b = new SAMRecordSetBuilder();
	// second argument is the contig index which is 0 based.  So contig index=0 -> chr1.  index=2 -> chr3, etc.
	b.addFrag("chr1", 0, 1, false);
	b.addFrag("chr2", 1, 1, false);
	b.addFrag("chr3", 2, 1, false);
	b.addFrag("chr4", 3, 1, false);
	Iterator<SAMRecord> recs = b.getRecords().iterator();

	Assert.assertTrue(f.filterOut(recs.next()));
	Assert.assertTrue(f.filterOut(recs.next()));
	// 3 and 4 not on the list.
	Assert.assertFalse(f.filterOut(recs.next()));
	Assert.assertFalse(f.filterOut(recs.next()));


}
 

private List<SAMRecord> getRecords () {
	SAMRecordSetBuilder builder  = new SAMRecordSetBuilder();

	// chromosome 2 and chromosome 10
	SAMRecord r1 = builder.addFrag("read1", 1, 1, false);
	Interval i1 = new Interval("2", 1, 10, true, null);
	r1.setAttribute(this.intervalTag, IntervalTagComparator.toString(i1));

	SAMRecord r2 = builder.addFrag("read1", 1, 1, false);
	Interval i2 = new Interval("10", 1, 10, true, null);
	r2.setAttribute(this.intervalTag, IntervalTagComparator.toString(i2));

	List<SAMRecord> result = new ArrayList<>();
	result.add(r1);
	result.add(r2);
	return result;
}
 

@Test
public void testSortIntervalsWithNames () {
	SAMRecordSetBuilder builder  = new SAMRecordSetBuilder();

	// chromosome 2 and chromosome 10
	SAMRecord r1 = builder.addFrag("read1", 1, 1, false);
	Interval i1 = new Interval("2", 1, 10, true, "foo");
	r1.setAttribute(this.intervalTag, IntervalTagComparator.toString(i1));

	SAMRecord r2 = builder.addFrag("read1", 1, 1, false);
	Interval i2 = new Interval("2", 1, 10, true, "bar");
	r2.setAttribute(this.intervalTag, IntervalTagComparator.toString(i2));

	// bar comes before foo.
	int expected = "foo".compareTo("bar");
	IntervalTagComparator c = new IntervalTagComparator(this.intervalTag);
	int comp = c.compare(r1, r2);
	Assert.assertEquals(expected,  comp);
}
 

@Test
public void testBasicHalfmappedReadPair() {
    SAMRecordSetBuilder builder = new SAMRecordSetBuilder(false, SAMFileHeader.SortOrder.queryname);
    builder.setReadLength(READ_LENGTH);
    builder.addPair("halfmapped_paired", 1, 1, 31, false, true, "20M", "20M", true, false, 20);
    PeekableIterator<SAMRecord> iterator = new PeekableIterator<SAMRecord>(builder.iterator());

    QuerySortedReadPairIteratorUtil.ReadPair pair = QuerySortedReadPairIteratorUtil.getNextReadPair(iterator);
    Assert.assertNotNull(pair);
    Assert.assertNotNull(pair.read1);
    Assert.assertNotNull(pair.read2);
    Assert.assertEquals("halfmapped_paired", pair.read1.getReadName());
    Assert.assertEquals("halfmapped_paired", pair.read2.getReadName());

    pair = QuerySortedReadPairIteratorUtil.getNextReadPair(iterator);
    Assert.assertNull(pair);
}
 

@Test
public void testFragmentNoReadPair() {
    SAMRecordSetBuilder builder = new SAMRecordSetBuilder(false, SAMFileHeader.SortOrder.queryname);
    builder.setReadLength(READ_LENGTH);
    builder.addFrag("mapped_frag_a", 1, 1, false);
    builder.addFrag("mapped_frag_b", 1, 1, false);
    PeekableIterator<SAMRecord> iterator = new PeekableIterator<SAMRecord>(builder.iterator());

    QuerySortedReadPairIteratorUtil.ReadPair pair = QuerySortedReadPairIteratorUtil.getNextReadPair(iterator);
    Assert.assertNotNull(pair);
    Assert.assertNotNull(pair.read1);
    Assert.assertNull(pair.read2);
    Assert.assertEquals("mapped_frag_a", pair.read1.getReadName());

    pair = QuerySortedReadPairIteratorUtil.getNextReadPair(iterator);
    Assert.assertNotNull(pair);
    Assert.assertNotNull(pair.read1);
    Assert.assertNull(pair.read2);
    Assert.assertEquals("mapped_frag_b", pair.read1.getReadName());

    pair = QuerySortedReadPairIteratorUtil.getNextReadPair(iterator);
    Assert.assertNull(pair);
}
 

/** Writes the contents of a SAMRecordSetBuilder out to a file. */
private File writeBam(final SAMRecordSetBuilder builder, final File f) {
    try (final SAMFileWriter out = new SAMFileWriterFactory().makeSAMOrBAMWriter(builder.getHeader(), false, f)) {
        builder.forEach(out::addAlignment);
    }
    return f;
}
 

@BeforeTest
void setupBuilder() throws IOException {
    tempSamFileChrM_O = VcfTestUtils.createTemporaryIndexedFile("CollectGcBias", ".bam");
    tempSamFileAllChr = VcfTestUtils.createTemporaryIndexedFile("CollectGcBias", ".bam");

    final File tempSamFileUnsorted = VcfTestUtils.createTemporaryIndexedFile("CollectGcBias", ".bam");

    final SAMFileHeader header = new SAMFileHeader();

    try {
        header.setSequenceDictionary(SAMSequenceDictionaryExtractor.extractDictionary(dict.toPath()));
        header.setSortOrder(SAMFileHeader.SortOrder.unsorted);
    } catch (final SAMException e) {
        e.printStackTrace();
    }

    //build different levels to put into the same bam file for testing multi level collection
    setupTest1(1, readGroupId1, readGroupRecord1, sample1, library1, header, setBuilder1); //Sample 1, Library 1, RG 1
    setupTest1(2, readGroupId2, readGroupRecord2, sample1, library2, header, setBuilder2); //Sample 1, Library 2, RG 2
    setupTest1(3, readGroupId3, readGroupRecord3, sample2, library3, header, setBuilder3); //Sample 2, Library 3, RG 3

    //build one last readgroup for comparing that window count stays the same whether you use all contigs or not
    setupTest2(1, readGroupId1, readGroupRecord1, sample1, library1, header, setBuilder4);

    final List<SAMRecordSetBuilder> test1Builders = new ArrayList<>();
    test1Builders.add(setBuilder1);
    test1Builders.add(setBuilder2);
    test1Builders.add(setBuilder3);

    final List<SAMRecordSetBuilder> test2Builders = new ArrayList<>();
    test2Builders.add(setBuilder4);

    tempSamFileChrM_O = build(test1Builders, tempSamFileUnsorted, header);
    tempSamFileAllChr = build(test2Builders, tempSamFileUnsorted, header);
}
 

@Test(dataProvider = "readPositionIgnoringSoftClips")
public void testGetReadPositionIgnoringSoftClips(final String cigarString, final int startPosition, final int queryPosition, final int expectedReadPosititon) {
    final SAMFileHeader newHeader = SAMRecordSetBuilder.makeDefaultHeader(SAMFileHeader.SortOrder.queryname, 100000,false);
    final SAMRecord rec = new SAMRecord(newHeader);

    rec.setCigarString(cigarString);
    rec.setAlignmentStart(startPosition);

    final int readPosition = AbstractAlignmentMerger.getReadPositionAtReferencePositionIgnoreSoftClips(rec, queryPosition);

    Assert.assertEquals(readPosition, expectedReadPosititon);
}
 

public void setupTest1(final int ID, final String readGroupId, final SAMReadGroupRecord readGroupRecord, final String sample,
                  final String library, final SAMFileHeader header, final SAMRecordSetBuilder setBuilder)
        throws IOException {

    final String separator = ":";
    final int contig1 = 0;
    final int contig2 = 1;
    readGroupRecord.setSample(sample);
    readGroupRecord.setPlatform(platform);
    readGroupRecord.setLibrary(library);
    readGroupRecord.setPlatformUnit(readGroupId);
    header.addReadGroup(readGroupRecord);
    setBuilder.setReadGroup(readGroupRecord);
    setBuilder.setUseNmFlag(true);

    setBuilder.setHeader(header);

    final int max = 800;
    final int min = 1;
    final Random rg = new Random(5);

    //add records that align to chrM and O but not N
    for (int i = 0; i < NUM_READS; i++) {
        final int start = rg.nextInt(max) + min;
        final String newReadName = READ_NAME + separator + ID + separator + i;

        if (i != NUM_READS - 1) {
            setBuilder.addPair(newReadName, contig1, start + ID, start + ID + LENGTH);
        } else {
            setBuilder.addPair(newReadName, contig2, start + ID, start + ID + LENGTH);
        }
    }
}
 

public void setupTest2(final int ID, final String readGroupId, final SAMReadGroupRecord readGroupRecord, final String sample,
                       final String library, final SAMFileHeader header, final SAMRecordSetBuilder setBuilder)
        throws IOException {

    final String separator = ":";
    final int contig1 = 0;
    final int contig2 = 1;
    final int contig3 = 2;
    readGroupRecord.setSample(sample);
    readGroupRecord.setPlatform(platform);
    readGroupRecord.setLibrary(library);
    readGroupRecord.setPlatformUnit(readGroupId);
    setBuilder.setReadGroup(readGroupRecord);
    setBuilder.setUseNmFlag(true);

    setBuilder.setHeader(header);

    final int max = 800;
    final int min = 1;
    final Random rg = new Random(5);

    //add records that align to all 3 chr in reference file
    for (int i = 0; i < NUM_READS; i++) {
        final int start = rg.nextInt(max) + min;
        final String newReadName = READ_NAME + separator + ID + separator + i;

        if (i<=NUM_READS/3) {
            setBuilder.addPair(newReadName, contig1, start + ID, start + ID + LENGTH);
        } else if (i< (NUM_READS - (NUM_READS/3))) {
            setBuilder.addPair(newReadName, contig2, start + ID, start + ID + LENGTH);
        } else {
            setBuilder.addPair(newReadName, contig3, start + ID, start + ID + LENGTH);
        }
    }
}
 

void setup(final int numReads,
           final String readName,
           final int ID,
           final String readGroupId,
           final SAMReadGroupRecord readGroupRecord,
           final String sample,
           final String library,
           final SAMFileHeader header,
           final SAMRecordSetBuilder setBuilder) throws IOException {
    final String separator = ":";
    readGroupRecord.setSample(sample);
    readGroupRecord.setPlatform(platform);
    readGroupRecord.setLibrary(library);
    readGroupRecord.setPlatformUnit(readGroupId);
    header.addReadGroup(readGroupRecord);
    setBuilder.setReadGroup(readGroupRecord);
    setBuilder.setUseNmFlag(true);
    setBuilder.setHeader(header);

    final int max = 15000;
    final int min = 1;
    final Random rg = new Random(5);

    for (int i = 0; i < numReads; i++) {
        final int start = rg.nextInt(max) + min;
        final String newReadName = readName + separator + ID + separator + i;
        setBuilder.addPair(newReadName, 0, start + ID, start + ID + 99);
    }
}
 

public static File writeBamFileWithLargeHeader() throws IOException {
  SAMRecordSetBuilder samRecordSetBuilder =
      new SAMRecordSetBuilder(true, SAMFileHeader.SortOrder.queryname);
  for (int i = 0; i < 1000; i++) {
    int chr = 20;
    int start1 = (i + 1) * 1000;
    int start2 = start1 + 100;
    samRecordSetBuilder.addPair(String.format("test-read-%03d", i), chr, start1,
        start2);
  }

  final File bamFile = File.createTempFile("test", ".bam");
  bamFile.deleteOnExit();
  SAMFileHeader samHeader = samRecordSetBuilder.getHeader();
  StringBuffer sb = new StringBuffer();
  for (int i = 0; i < 1000000; i++) {
    sb.append("0123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789");
  }
  samHeader.addComment(sb.toString());
  final SAMFileWriter bamWriter = new SAMFileWriterFactory()
      .makeSAMOrBAMWriter(samHeader, true, bamFile);
  for (final SAMRecord rec : samRecordSetBuilder.getRecords()) {
    bamWriter.addAlignment(rec);
  }
  bamWriter.close();

  return bamFile;
}
 

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

}
 

@Test
public void filterByReadNumberTest() {
	FilterBamByTag t = new FilterBamByTag();

	// record paired and read is 1st
	List<SAMRecord> recs = getPairedRead ();
	SAMRecord recFirstPaired = recs.get(0);
	SAMRecord recSecondPaired = recs.get(1);

	boolean flag1= t.retainByReadNumber(recFirstPaired, 1);
	boolean flag2= t.retainByReadNumber(recFirstPaired, 2);
	Assert.assertTrue(flag1);
	Assert.assertFalse(flag2);

	// record paired and read is 2st
	recSecondPaired.setProperPairFlag(true);
	recSecondPaired.setSecondOfPairFlag(true);
	flag1= t.retainByReadNumber(recSecondPaired, 1);
	flag2= t.retainByReadNumber(recSecondPaired, 2);
	Assert.assertTrue(flag2);
	Assert.assertFalse(flag1);

	// record unpaired and read is 1st
	SAMRecordSetBuilder builder = new SAMRecordSetBuilder();
	builder.addUnmappedFragment("foo");
	SAMRecord recFirstUnPaired = builder.getRecords().iterator().next();

	flag1= t.retainByReadNumber(recFirstUnPaired, 1);
	flag2= t.retainByReadNumber(recFirstPaired, 2);
	Assert.assertTrue(flag1);
	Assert.assertFalse(flag2);
}
 

@Test
public void testBaseErrorAggregation2() {
    final SAMSequenceRecord samSequenceRecord = new SAMSequenceRecord("chr1", 200);
    final SAMFileHeader samFileHeader = new SAMFileHeader();
    samFileHeader.addSequence(samSequenceRecord);

    final SAMRecordSetBuilder builder = new SAMRecordSetBuilder();

    final List<SAMRecord> samRecords = builder.addPair("Read1234", 0, 1, 1,
            false, false, "16M", "16M", true, false, 20);
    final SAMRecord samRecord1 = samRecords.get(0);
    final SAMRecord samRecord2 = samRecords.get(1);

    samRecord1.setReadBases("CgTGtGGAcAAAgAAA".getBytes());
    samRecord2.setReadBases("CcTGGtGAcAAAgAAA".getBytes());
    final byte[] refBases = "CATGGGGAAAAAAAAA".getBytes();

    BaseErrorAggregation<?> baseErrorAggregation = new BaseErrorAggregation<>(OverlappingReadsErrorCalculator::new, ReadBaseStratification.readDirectionStratifier);

    final int length = getLengthAndAddBases(samSequenceRecord, samRecord1, samRecord2, refBases, baseErrorAggregation);

    final ErrorMetric[] metrics = baseErrorAggregation.getMetrics();
    final OverlappingErrorMetric metric1 = (OverlappingErrorMetric) metrics[0];
    metric1.calculateDerivedFields();
    Assert.assertEquals(metric1.COVARIATE, ReadBaseStratification.ReadDirection.POSITIVE.toString());
    Assert.assertEquals(metric1.TOTAL_BASES, length);
    Assert.assertEquals(metric1.NUM_BASES_WITH_OVERLAPPING_READS, length);
    Assert.assertEquals(metric1.NUM_DISAGREES_WITH_REFERENCE_ONLY, 2L);
    Assert.assertEquals(metric1.NUM_DISAGREES_WITH_REF_AND_MATE, 1L);
    Assert.assertEquals(metric1.NUM_THREE_WAYS_DISAGREEMENT, 1L);

    final OverlappingErrorMetric metric2 = (OverlappingErrorMetric) metrics[1];
    metric2.calculateDerivedFields();
    Assert.assertEquals(metric2.COVARIATE, ReadBaseStratification.ReadDirection.NEGATIVE.toString());
    Assert.assertEquals(metric2.TOTAL_BASES, length);
    Assert.assertEquals(metric2.NUM_BASES_WITH_OVERLAPPING_READS, length);
    Assert.assertEquals(metric2.NUM_DISAGREES_WITH_REFERENCE_ONLY, 2L);
    Assert.assertEquals(metric2.NUM_DISAGREES_WITH_REF_AND_MATE, 1L);
    Assert.assertEquals(metric2.NUM_THREE_WAYS_DISAGREEMENT, 1L);
}
 

public GenerateIntervalTaggedSAMRecords (final SAMSequenceDictionary dict, final String intervalTag, final int readLength, final int maxNumRecords) {
	this.intervalTag=intervalTag;
	this.readLength=readLength;
	this.maxNumRecords=maxNumRecords;
	this.recs = dict.getSequences();
	this.numContigs=recs.size();
	randomGenerator = new Random();
	this.currentNumRecords=0;


	final SAMRecordSetBuilder ret = new SAMRecordSetBuilder(false, SAMFileHeader.SortOrder.coordinate);
	ret.addUnmappedFragment("foo");
	this.samRecordTemplate= ret.getRecords().iterator().next();

}
 

@Test
public void testBaseErrorAggregation() {
    final SAMSequenceRecord samSequenceRecord = new SAMSequenceRecord("chr1", 200);
    final SAMFileHeader samFileHeader = new SAMFileHeader();
    samFileHeader.addSequence(samSequenceRecord);

    final SAMRecordSetBuilder builder = new SAMRecordSetBuilder();

    final List<SAMRecord> samRecords = builder.addPair("Read1234", 0, 1, 1,
            false, false, "16M", "16M", true, false, 20);
    final SAMRecord samRecord1 = samRecords.get(0);
    final SAMRecord samRecord2 = samRecords.get(1);

    samRecord1.setReadBases("CgTGtGGAcAAAgAAA".getBytes());
    samRecord2.setReadBases("CcTGGtGAcAAAgAAA".getBytes());
    final byte[] refBases = "CATGGGGAAAAAAAAA".getBytes();

    BaseErrorAggregation<?> baseErrorAggregation = new BaseErrorAggregation<>(OverlappingReadsErrorCalculator::new, ReadBaseStratification.readOrdinalityStratifier);

    final int length = getLengthAndAddBases(samSequenceRecord, samRecord1, samRecord2, refBases, baseErrorAggregation);

    final ErrorMetric[] metrics = baseErrorAggregation.getMetrics();
    final OverlappingErrorMetric metric1 = (OverlappingErrorMetric) metrics[0];
    metric1.calculateDerivedFields();
    Assert.assertEquals(metric1.COVARIATE, ReadBaseStratification.ReadOrdinality.FIRST.name());
    Assert.assertEquals(metric1.TOTAL_BASES, length);
    Assert.assertEquals(metric1.NUM_BASES_WITH_OVERLAPPING_READS, length);
    Assert.assertEquals(metric1.NUM_DISAGREES_WITH_REFERENCE_ONLY, 2L);
    Assert.assertEquals(metric1.NUM_DISAGREES_WITH_REF_AND_MATE, 1L);
    Assert.assertEquals(metric1.NUM_THREE_WAYS_DISAGREEMENT, 1L);

    final OverlappingErrorMetric metric2 = (OverlappingErrorMetric) metrics[1];
    metric2.calculateDerivedFields();
    Assert.assertEquals(metric2.COVARIATE, ReadBaseStratification.ReadOrdinality.SECOND.name());
    Assert.assertEquals(metric2.TOTAL_BASES, length);
    Assert.assertEquals(metric2.NUM_BASES_WITH_OVERLAPPING_READS, length);
    Assert.assertEquals(metric2.NUM_DISAGREES_WITH_REFERENCE_ONLY, 2L);
    Assert.assertEquals(metric2.NUM_DISAGREES_WITH_REF_AND_MATE, 1L);
    Assert.assertEquals(metric2.NUM_THREE_WAYS_DISAGREEMENT, 1L);
}
 

@Test ()
public void testScoreBaseQualityWrongLength() {
	SAMRecordSetBuilder b = new SAMRecordSetBuilder();
	// This read is plenty long and high quality.
	b.addFrag("read1", 0, 1, false, false, "20M", null, 30);
	// this read is too short.
	b.addFrag("read2", 0, 1, false, false, "5M", null, 30);
	Iterator<SAMRecord> iter = b.getRecords().iterator();
	SAMRecord r = iter.next();

	List <BaseRange> baseRangeList = new ArrayList<>();
	baseRangeList.add(new BaseRange(1, 12));

	BaseQualityFilter bqf = new BaseQualityFilter(baseRangeList, 20);
	int quality = bqf.scoreBaseQuality(r);
	Assert.assertEquals(quality, 0);

	// give a read that's too short.
	// this causes a null pointer exception, let's make it a nicer reported exception.
	r = iter.next();

	try {
		quality = bqf.scoreBaseQuality(r);
	} catch (TranscriptomeException e) {
		System.out.println(e.getMessage());
	}


	Assert.assertEquals(quality, 0);
}
 

@Test
public void filterOut() {
	String strandTag = "GS";
	Iterator<SAMRecord> underlyingIterator = Collections.emptyIterator();
	GeneStrandFilteringIterator f = new GeneStrandFilteringIterator(underlyingIterator, strandTag);

	SAMRecordSetBuilder b = new SAMRecordSetBuilder();
	// positive strand read.
	b.addFrag("1", 1, 1, false);
	b.addFrag("1", 1, 1, true);
	Iterator<SAMRecord> recs = b.getRecords().iterator();

	SAMRecord r1 = recs.next();
	SAMRecord r2 = recs.next();

	// test all 4 possibilities.
	r1.setAttribute(strandTag, "+");
	boolean t1 = f.filterOut(r1);
	Assert.assertFalse(t1);

	r1.setAttribute(strandTag, "-");
	t1 = f.filterOut(r1);
	Assert.assertTrue(t1);


	r2.setAttribute(strandTag, "+");
	t1 = f.filterOut(r2);
	Assert.assertTrue(t1);

	r2.setAttribute(strandTag, "-");
	t1 = f.filterOut(r2);
	Assert.assertFalse(t1);

}
 

@Test
    public void testUnmapBacterialContamination() throws IOException {
        final SAMRecordSetBuilder builder = new SAMRecordSetBuilder(true, SAMFileHeader.SortOrder.queryname);
        final SAMFileHeader header = builder.getHeader();
        final SAMFileHeader.SortOrder sortOrder = header.getSortOrder();
        final SAMFileHeader newHeader = SAMRecordSetBuilder.makeDefaultHeader(sortOrder, 100000,true);
        builder.setHeader(newHeader);

        final File reference = File.createTempFile("reference",".fasta");
        reference.deleteOnExit();

        builder.writeRandomReference(reference.toPath());

        builder.addPair("overlappingpair", 0,500,500, false,false,"20S20M60S","20S20M60M",true,false,45);
        builder.addPair("overlappingpairFirstAligned", 0,500,500, false,true,"20S20M60S",null,true,false,45);
        builder.addPair("overlappingpairSecondAligned", 0,500,500, true,false,null,"20S20M60S",true,false,45);
        builder.addPair("overlappingpairFirstAlignedB", 0,500,500, false,true,"20S20M60S",null,false,true,45);
        builder.addPair("overlappingpairSecondAlignedB", 0,500,500, true,false,null,"20S20M60S",false,true,45);

//        builder.addFrag("frag",1,500,false,false,"20S20M60S",null, 45);
//        builder.addFrag("frag2",1,500,true,false,"20S20M60S",null, 45);
//        builder.addFrag("frag3",1,500,false,false,"20S20M60S",null, 45);
//        builder.addFrag("frag4",1,500,true,false,"20S20M60S",null, 45);

        final File file = newTempSamFile("aligned");

        try (SAMFileWriter writer = new SAMFileWriterFactory().makeWriter(builder.getHeader(), true, file, null)) {
            builder.getRecords().forEach(writer::addAlignment);
        }

        final RevertSam revertSam = new RevertSam();

        revertSam.INPUT = file;
        final File fileUnaligned = newTempSamFile("unaligned");
        revertSam.OUTPUT = fileUnaligned;

        revertSam.SANITIZE = false;
        revertSam.REMOVE_ALIGNMENT_INFORMATION=true;
        revertSam.REMOVE_DUPLICATE_INFORMATION=true;

        revertSam.SORT_ORDER = SAMFileHeader.SortOrder.queryname;

        Assert.assertEquals(revertSam.doWork(),0);

        MergeBamAlignment mergeBamAlignment = new MergeBamAlignment();

        mergeBamAlignment.ALIGNED_BAM = Collections.singletonList(file);
        mergeBamAlignment.UNMAPPED_BAM = fileUnaligned;
        mergeBamAlignment.UNMAP_CONTAMINANT_READS = true;

        //yuck!
        final RequiredReferenceArgumentCollection requiredReferenceArgumentCollection = new RequiredReferenceArgumentCollection();
        requiredReferenceArgumentCollection.REFERENCE_SEQUENCE = reference;
        mergeBamAlignment.referenceSequence = requiredReferenceArgumentCollection;

        final File fileMerged = newTempSamFile("merged");

        mergeBamAlignment.OUTPUT = fileMerged;

        // merge file with itself.
        Assert.assertEquals(mergeBamAlignment.doWork(),0);

        // check that all reads have been unmapped due to bacterial contamination as needed.
        try (SamReader mergedReader = SamReaderFactory.makeDefault().open(fileMerged)) {
            for (SAMRecord mergedRecord : mergedReader) {
                Assert.assertTrue(mergedRecord.getReadUnmappedFlag(), mergedRecord.toString());
                Assert.assertTrue(!mergedRecord.getReadPairedFlag() || mergedRecord.getMateUnmappedFlag(), mergedRecord.toString());
            }
        }
    }