Java Code Examples for htsjdk.samtools.SAMRecordSetBuilder#addFrag()

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

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

}
 

@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 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 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 testFilter() {
	Iterator<SAMRecord> underlyingIterator = Collections.emptyIterator();
	String tagName = "XC";
	// filter always returns false.
	BamTagCountingIterator f = new BamTagCountingIterator(underlyingIterator, tagName);


	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, "A");
	// this one counts.
	Assert.assertFalse(f.filterOut(r1));

	// wrong tag set, correct tag is not set.
	r1.setAttribute("XZ", "A");
	r1.setAttribute(tagName, null);
	Assert.assertFalse(f.filterOut(r1));

	r1.setAttribute(tagName, "A");
	// this one counts.
	Assert.assertFalse(f.filterOut(r1));

	r1.setAttribute(tagName, "B");
	// this one counts.
	Assert.assertFalse(f.filterOut(r1));

	r1.setAttribute(tagName, "C");
	// this one counts.
	Assert.assertFalse(f.filterOut(r1));

	ObjectCounter<String> expected = f.getCounts();
	expected.incrementByCount("A", 2);
	expected.incrementByCount("B", 1);
	expected.incrementByCount("C", 1);
	ObjectCounter<String> actual = f.getCounts();

	Assert.assertEquals(actual, expected);



}
 

public SAMRecord createRecordFromCigar(final String cigar) {
    final SAMRecordSetBuilder builder = new SAMRecordSetBuilder();
    builder.addFrag("", 0, 100, false, false, cigar, null, 30);
    return builder.getRecords().stream().findFirst().get();
}
 

@Test(dataProvider = "provideForTestIndelErrors")
public void testIndelErrors(final String[] readCigars, final String errorSubscript, IndelErrorMetric expectedMetric) throws IOException {

    final File temp = File.createTempFile("Indels", ".bam");
    temp.deleteOnExit();

    final int priorQ = 30;

    try (final ReferenceSequenceFileWalker referenceSequenceFileWalker =
                    new ReferenceSequenceFileWalker(CommandLineProgramTest.CHR_M_REFERENCE)) {

        final SAMRecordSetBuilder builder = new SAMRecordSetBuilder();
        builder.getHeader().setSequenceDictionary(referenceSequenceFileWalker.getSequenceDictionary());

        for (int i = 0; i < readCigars.length; i++) {
            // 10M2I3M4D10M5I
            builder.addFrag("Read" + i, 0, 100, false, false, readCigars[i], null, 30);
        }

        try (final SAMFileWriter writer = new SAMFileWriterFactory()
                .setCompressionLevel(2)
                .makeBAMWriter(builder.getHeader(), false, temp)) {
            builder.forEach(writer::addAlignment);
        }
    }

    final File referenceFile = CHR_M_REFERENCE;
    final File vcf = new File(TEST_DIR, "NIST.selected.vcf");

    final File outputBaseFileName = new File(OUTPUT_DATA_PATH, "test");
    final File errorByAll = new File(outputBaseFileName.getAbsolutePath() + errorSubscript);
    errorByAll.deleteOnExit();
    outputBaseFileName.deleteOnExit();

    final String[] args = {
            "INPUT=" + temp,
            "OUTPUT=" + outputBaseFileName,
            "REFERENCE_SEQUENCE=" + referenceFile.getAbsolutePath(),
            "ERROR_METRICS=" + "INDEL_ERROR",
            "ERROR_METRICS=" + "INDEL_ERROR:INDEL_LENGTH",
            "VCF=" + vcf.getAbsolutePath()
    };

    CollectSamErrorMetrics collectSamErrorMetrics = new CollectSamErrorMetrics();
    collectSamErrorMetrics.ERROR_METRICS.clear();

    Assert.assertEquals(collectSamErrorMetrics.instanceMain(args), 0);

    ErrorMetric.setPriorError(QualityUtil.getErrorProbabilityFromPhredScore(priorQ));
    expectedMetric.calculateDerivedFields();

    // Note that soft clipped bases are not counted
    List<IndelErrorMetric> metrics = MetricsFile.readBeans(errorByAll);

    IndelErrorMetric metric = metrics
            .stream()
            .filter(m -> m.COVARIATE.equals(expectedMetric.COVARIATE))
            .findAny()
            .orElseThrow(() -> new AssertionError("didn't find metric with COVARIATE==" + expectedMetric.COVARIATE));

    Assert.assertEquals(metric, expectedMetric);
}
 

@Test
public void testHsMetricsHandlesIndelsAppropriately() throws IOException {
    final SAMRecordSetBuilder withDeletions = new SAMRecordSetBuilder(true, SortOrder.coordinate);
    final SAMRecordSetBuilder withInsertions = new SAMRecordSetBuilder(true, SortOrder.coordinate);
    final IntervalList targets = new IntervalList(withDeletions.getHeader());
    final IntervalList baits   = new IntervalList(withDeletions.getHeader());
    targets.add(new Interval("chr1", 1000, 1199, false, "t1"));
    baits.add(new Interval("chr1", 950,  1049, false, "b1"));
    baits.add(new Interval("chr1", 1050, 1149, false, "b2"));
    baits.add(new Interval("chr1", 1150, 1249, false, "b3"));

    // Generate 100 reads that fully cover the the target in each BAM
    for (int i=0; i<100; ++i) {
        withDeletions.addFrag( "d" + i, 0, 1000, false, false, "100M20D80M", null, 30);
        withInsertions.addFrag("i" + i, 0, 1000, false, false, "100M50I100M", null, 30);
    }

    // Write things out to file
    final File dir = IOUtil.createTempDir("hsmetrics.", ".test");
    final File bs = new File(dir, "baits.interval_list").getAbsoluteFile();
    final File ts = new File(dir, "targets.interval_list").getAbsoluteFile();
    baits.write(bs);
    targets.write(ts);
    final File withDelBam = writeBam(withDeletions,  new File(dir, "with_del.bam"));
    final File withInsBam = writeBam(withInsertions, new File(dir, "with_ins.bam"));

    // Now run CollectHsMetrics four times
    final File out = Files.createTempFile("hsmetrics.", ".txt").toFile();
    runPicardCommandLine(Arrays.asList("INCLUDE_INDELS=false", "SAMPLE_SIZE=0", "TI="+ts.getPath(), "BI="+bs.getPath(), "O="+out.getPath(), "I="+withDelBam.getAbsolutePath()));
    final HsMetrics delsWithoutIndelHandling = readMetrics(out);
    runPicardCommandLine(Arrays.asList("INCLUDE_INDELS=true", "SAMPLE_SIZE=0", "TI="+ts.getPath(), "BI="+bs.getPath(), "O="+out.getPath(), "I="+withDelBam.getAbsolutePath()));
    final HsMetrics delsWithIndelHandling = readMetrics(out);
    runPicardCommandLine(Arrays.asList("INCLUDE_INDELS=false", "SAMPLE_SIZE=0", "TI="+ts.getPath(), "BI="+bs.getPath(), "O="+out.getPath(), "I="+withInsBam.getAbsolutePath()));
    final HsMetrics insWithoutIndelHandling = readMetrics(out);
    runPicardCommandLine(Arrays.asList("INCLUDE_INDELS=true", "SAMPLE_SIZE=0", "TI="+ts.getPath(), "BI="+bs.getPath(), "O="+out.getPath(), "I="+withInsBam.getAbsolutePath()));
    final HsMetrics insWithIndelHandling = readMetrics(out);

    IOUtil.deleteDirectoryTree(dir);

    Assert.assertEquals(delsWithoutIndelHandling.MEAN_TARGET_COVERAGE, 90.0);  // 100X over 180/200 bases due to deletion
    Assert.assertEquals(delsWithIndelHandling.MEAN_TARGET_COVERAGE, 100.0);    // 100X with counting the deletion

    Assert.assertEquals(insWithoutIndelHandling.PCT_USABLE_BASES_ON_TARGET, 200/250d); // 50/250 inserted bases are not counted as on target
    Assert.assertEquals(insWithIndelHandling.PCT_USABLE_BASES_ON_TARGET,   1.0d);      // inserted bases are counted as on target
}
 

public static File writeBamFile(int numPairs, SAMFileHeader.SortOrder sortOrder)
    throws IOException {
  // file will be both queryname and coordinate sorted, so use one or the other
  SAMRecordSetBuilder samRecordSetBuilder = new SAMRecordSetBuilder(true, sortOrder);
  for (int i = 0; i < numPairs; i++) {
    int chr = 20;
    int start1 = (i + 1) * 1000;
    int start2 = start1 + 100;
    if (i == 5) { // add two unmapped fragments instead of a mapped pair
      samRecordSetBuilder.addFrag(String.format("test-read-%03d-1", i), chr, start1,
          false, true, null,
          null,
          -1, false);
      samRecordSetBuilder.addFrag(String.format("test-read-%03d-2", i), chr, start2,
          false, true, null,
          null,
          -1, false);
    } else {
      samRecordSetBuilder.addPair(String.format("test-read-%03d", i), chr, start1,
          start2);
    }
  }
  if (numPairs > 0) { // add two unplaced unmapped fragments if non-empty
    samRecordSetBuilder.addUnmappedFragment(String.format
        ("test-read-%03d-unplaced-unmapped", numPairs++));
    samRecordSetBuilder.addUnmappedFragment(String.format
        ("test-read-%03d-unplaced-unmapped", numPairs++));
  }

  final File bamFile = File.createTempFile("test", ".bam");
  bamFile.deleteOnExit();
  SAMFileHeader samHeader = samRecordSetBuilder.getHeader();
  final SAMFileWriter bamWriter = new SAMFileWriterFactory()
      .makeSAMOrBAMWriter(samHeader, true, bamFile);
  for (final SAMRecord rec : samRecordSetBuilder.getRecords()) {
    bamWriter.addAlignment(rec);
  }
  bamWriter.close();

  // create BAM index
  if (sortOrder.equals(SAMFileHeader.SortOrder.coordinate)) {
    SamReader samReader = SamReaderFactory.makeDefault()
        .enable(SamReaderFactory.Option.INCLUDE_SOURCE_IN_RECORDS)
        .open(bamFile);
    BAMIndexer.createIndex(samReader, new File(bamFile.getAbsolutePath()
        .replaceFirst("\\.bam$", BAMIndex.BAMIndexSuffix)));
  }

  return bamFile;
}