htsjdk.samtools.SamPairUtil Java Examples

The following examples show how to use htsjdk.samtools.SamPairUtil. You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example #1
Source File: InsertSizeMetricsCollector.java    From picard with MIT License 6 votes vote down vote up
public PerUnitInsertSizeMetricsCollector(final String sample, final String library, final String readGroup) {
    this.sample = sample;
    this.library = library;
    this.readGroup = readGroup;
    String prefix = null;
    if (this.readGroup != null) {
        prefix = this.readGroup + ".";
    }
    else if (this.library != null) {
        prefix = this.library + ".";
    }
    else if (this.sample != null) {
        prefix = this.sample + ".";
    }
    else {
        prefix = "All_Reads.";
    }
    histograms.put(SamPairUtil.PairOrientation.FR,     new Histogram<Integer>("insert_size", prefix + "fr_count"));
    histograms.put(SamPairUtil.PairOrientation.TANDEM, new Histogram<Integer>("insert_size", prefix + "tandem_count"));
    histograms.put(SamPairUtil.PairOrientation.RF,     new Histogram<Integer>("insert_size", prefix + "rf_count"));
}
 
Example #2
Source File: PerUnitInsertSizeMetricsCollector.java    From gatk with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
public PerUnitInsertSizeMetricsCollector(
        final String sample,
        final String library,
        final String readGroup,
        final double minimumPct,
        final double deviations,
        final Integer histogramWidth) {

    this.sample = sample;
    this.library = library;
    this.readGroup = readGroup;
    this.minimumPct = minimumPct;
    this.histogramWidth = histogramWidth;
    this.deviations = deviations;

    String prefix = createHistogramValuePrefix();

    histograms.put(SamPairUtil.PairOrientation.FR,     new Histogram<>("insert_size", prefix + "fr_count"));
    histograms.put(SamPairUtil.PairOrientation.TANDEM, new Histogram<>("insert_size", prefix + "tandem_count"));
    histograms.put(SamPairUtil.PairOrientation.RF,     new Histogram<>("insert_size", prefix + "rf_count"));
}
 
Example #3
Source File: MergeBamAlignmentTest.java    From picard with MIT License 6 votes vote down vote up
@Test(dataProvider="data")
public void testSortingOnSamAlignmentMerger(final File unmapped, final File aligned, final boolean sorted, final boolean coordinateSorted, final String testName)
    throws IOException {

    final File target = File.createTempFile("target", "bam");
    target.deleteOnExit();
    final SamAlignmentMerger merger = new SamAlignmentMerger(unmapped,  target, fasta, null, true, false,
            false, Collections.singletonList(aligned), 1, null, null, null, null, null, null,
            Collections.singletonList(SamPairUtil.PairOrientation.FR),
            coordinateSorted ? SAMFileHeader.SortOrder.coordinate : SAMFileHeader.SortOrder.queryname,
            new BestMapqPrimaryAlignmentSelectionStrategy(), false, false, 30);

    merger.mergeAlignment(Defaults.REFERENCE_FASTA);
    Assert.assertEquals(sorted, !merger.getForceSort());
    final SAMRecordIterator it = SamReaderFactory.makeDefault().open(target).iterator();
    int aln = 0;
    while (it.hasNext()) {
        final SAMRecord rec = it.next();
        if (!rec.getReadUnmappedFlag()) {
            aln++;
        }
    }
    Assert.assertEquals(aln, 6, "Incorrect number of aligned reads in merged BAM file");
}
 
Example #4
Source File: MergeBamAlignmentTest.java    From picard with MIT License 6 votes vote down vote up
private void doMergeAlignment(final File unmappedBam, final List<File> alignedBams,
                              final List<File> read1AlignedBams, final List<File> read2AlignedBams, final Integer read1Trim, final Integer read2Trim,
                              final boolean alignReadsOnly, final boolean clipAdapters, final boolean isBisulfiteSequence, final int maxInsOrDels,
                              final String progRecordId, final String progGroupVersion, final String progGroupCommandLine, final String progGroupName,
                              final boolean pairedRun, final File refSeq, final File output,
                              final SamPairUtil.PairOrientation expectedOrientation, final MergeBamAlignment.PrimaryAlignmentStrategy primaryAlignmentStrategy,
                              final String attributesToRetain,
                              final Boolean includeSecondary,
                              final Boolean unmapContaminantReads,
                              final SAMFileHeader.SortOrder sortOrder) {
    doMergeAlignment(unmappedBam, alignedBams, read1AlignedBams, read2AlignedBams, read1Trim, read2Trim,
            alignReadsOnly, clipAdapters, isBisulfiteSequence, maxInsOrDels,
            progRecordId, progGroupVersion, progGroupCommandLine, progGroupName, pairedRun, refSeq, output,
            expectedOrientation, primaryAlignmentStrategy, attributesToRetain, includeSecondary, unmapContaminantReads,
            sortOrder, null);
}
 
Example #5
Source File: InsertSizeMetricsCollector.java    From picard with MIT License 5 votes vote down vote up
@Override
protected InsertSizeCollectorArgs makeArg(SAMRecord samRecord, ReferenceSequence refSeq) {
    final int insertSize = Math.abs(samRecord.getInferredInsertSize());
    final SamPairUtil.PairOrientation orientation = SamPairUtil.getPairOrientation(samRecord);

    return new InsertSizeCollectorArgs(insertSize, orientation);
}
 
Example #6
Source File: CollectJumpingLibraryMetrics.java    From picard with MIT License 5 votes vote down vote up
/**
 * 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;
}
 
Example #7
Source File: InsertSizeMetricsCollector.java    From gatk with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
@Override
protected InsertSizeMetricsCollectorArgs makeArg(SAMRecord samRecord, ReferenceSequence refSeq) {
    // inferred insert size is negative if the mate maps to lower position than the read, so use abs
    final int insertSize = Math.abs(samRecord.getInferredInsertSize());
    final SamPairUtil.PairOrientation orientation = SamPairUtil.getPairOrientation(samRecord);

    return new InsertSizeMetricsCollectorArgs(insertSize, orientation);
}
 
Example #8
Source File: MergeBamAlignmentTest.java    From picard with MIT License 5 votes vote down vote up
@Test
public void testRemoveNmMdAndUqOnOverlappingReads() throws IOException {
    final File output = File.createTempFile("testRemoveNmMdAndUqOnOverlappingReads", ".sam");
    output.deleteOnExit();
    doMergeAlignment(new File(TEST_DATA_DIR, "removetags.unmapped.sam"),
            Collections.singletonList(new File(TEST_DATA_DIR, "removetags.aligned.sam")),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true,
            new File(TEST_DATA_DIR, "removetags.fasta"), output,
            SamPairUtil.PairOrientation.FR, null,
            null, null, null, SAMFileHeader.SortOrder.queryname);

    final SamReader result = SamReaderFactory.makeDefault().open(output);
    for (final SAMRecord rec : result) {
        boolean hasTags = false;
        if (rec.getReadName().startsWith("CLIPPED")) {
            final String[] readNameFields = rec.getReadName().split(":");
            final int index = rec.getFirstOfPairFlag() ? 1 : 2;
            hasTags = Integer.parseInt(readNameFields[index]) == 1;
        }
        if (hasTags) {
            Assert.assertNull(rec.getAttribute("MD"));
            Assert.assertNull(rec.getAttribute("NM"));
        } else {
            Assert.assertNotNull(rec.getAttribute("MD"));
            Assert.assertNotNull(rec.getAttribute("NM"));
        }
    }
    result.close();
}
 
Example #9
Source File: MergeBamAlignmentTest.java    From picard with MIT License 5 votes vote down vote up
@Test(expectedExceptions = {IllegalStateException.class, PicardException.class})
public void testOldQuerynameSortFails() throws IOException {

    final File merged = File.createTempFile("merged", BamFileIoUtils.BAM_FILE_EXTENSION);
    merged.deleteOnExit();

    doMergeAlignment(badorderUnmappedBam, Collections.singletonList(badorderAlignedBam),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true, fasta, merged,
            SamPairUtil.PairOrientation.FR, null,
            null, null, null, null);
    Assert.fail("Merger should have failed because unmapped reads are not in queryname order but didn't");
}
 
Example #10
Source File: MergeBamAlignmentTest.java    From picard with MIT License 5 votes vote down vote up
/**
 * Test that clipping of FR reads for fragments shorter than read length happens only when it should.
 */
@Test
public void testShortFragmentClipping() throws Exception {
    final File output = File.createTempFile("testShortFragmentClipping", ".sam");
    output.deleteOnExit();
    doMergeAlignment(new File(TEST_DATA_DIR, "cliptest.unmapped.sam"),
            Collections.singletonList(new File(TEST_DATA_DIR, "cliptest.aligned.sam")),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true,
            new File(TEST_DATA_DIR, "cliptest.fasta"), output,
            SamPairUtil.PairOrientation.FR, null,
            null, null, null, null);

    final SamReader result = SamReaderFactory.makeDefault().open(output);
    final Map<String, SAMRecord> firstReadEncountered = new HashMap<String, SAMRecord>();

    for (final SAMRecord rec : result) {
        final SAMRecord otherEnd = firstReadEncountered.get(rec.getReadName());
        if (otherEnd == null) {
            firstReadEncountered.put(rec.getReadName(), rec);
        } else {
            final int fragmentStart = Math.min(rec.getAlignmentStart(), otherEnd.getAlignmentStart());
            final int fragmentEnd = Math.max(rec.getAlignmentEnd(), otherEnd.getAlignmentEnd());
            final String[] readNameFields = rec.getReadName().split(":");
            // Read name of each pair includes the expected fragment start and fragment end positions.
            final int expectedFragmentStart = Integer.parseInt(readNameFields[1]);
            final int expectedFragmentEnd = Integer.parseInt(readNameFields[2]);
            Assert.assertEquals(fragmentStart, expectedFragmentStart, rec.getReadName());
            Assert.assertEquals(fragmentEnd, expectedFragmentEnd, rec.getReadName());
        }
    }
    result.close();
}
 
Example #11
Source File: MergeBamAlignmentTest.java    From picard with MIT License 5 votes vote down vote up
/**
 * 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");
}
 
Example #12
Source File: Utils.java    From cramtools with Apache License 2.0 4 votes vote down vote up
/**
 * Copied from net.sf.picard.sam.SamPairUtil. This is a more permissive
 * version of the method, which does not reset alignment start and reference
 * for unmapped reads.
 * 
 * @param rec1
 * @param rec2
 * @param header
 */
public static void setLooseMateInfo(final SAMRecord rec1, final SAMRecord rec2, final SAMFileHeader header) {
	if (rec1.getReferenceName() != SAMRecord.NO_ALIGNMENT_REFERENCE_NAME
			&& rec1.getReferenceIndex() == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX)
		rec1.setReferenceIndex(header.getSequenceIndex(rec1.getReferenceName()));
	if (rec2.getReferenceName() != SAMRecord.NO_ALIGNMENT_REFERENCE_NAME
			&& rec2.getReferenceIndex() == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX)
		rec2.setReferenceIndex(header.getSequenceIndex(rec2.getReferenceName()));

	// If neither read is unmapped just set their mate info
	if (!rec1.getReadUnmappedFlag() && !rec2.getReadUnmappedFlag()) {

		rec1.setMateReferenceIndex(rec2.getReferenceIndex());
		rec1.setMateAlignmentStart(rec2.getAlignmentStart());
		rec1.setMateNegativeStrandFlag(rec2.getReadNegativeStrandFlag());
		rec1.setMateUnmappedFlag(false);
		rec1.setAttribute(SAMTag.MQ.name(), rec2.getMappingQuality());

		rec2.setMateReferenceIndex(rec1.getReferenceIndex());
		rec2.setMateAlignmentStart(rec1.getAlignmentStart());
		rec2.setMateNegativeStrandFlag(rec1.getReadNegativeStrandFlag());
		rec2.setMateUnmappedFlag(false);
		rec2.setAttribute(SAMTag.MQ.name(), rec1.getMappingQuality());
	}
	// Else if they're both unmapped set that straight
	else if (rec1.getReadUnmappedFlag() && rec2.getReadUnmappedFlag()) {
		rec1.setMateNegativeStrandFlag(rec2.getReadNegativeStrandFlag());
		rec1.setMateUnmappedFlag(true);
		rec1.setAttribute(SAMTag.MQ.name(), null);
		rec1.setInferredInsertSize(0);

		rec2.setMateNegativeStrandFlag(rec1.getReadNegativeStrandFlag());
		rec2.setMateUnmappedFlag(true);
		rec2.setAttribute(SAMTag.MQ.name(), null);
		rec2.setInferredInsertSize(0);
	}
	// And if only one is mapped copy it's coordinate information to the
	// mate
	else {
		final SAMRecord mapped = rec1.getReadUnmappedFlag() ? rec2 : rec1;
		final SAMRecord unmapped = rec1.getReadUnmappedFlag() ? rec1 : rec2;

		mapped.setMateReferenceIndex(unmapped.getReferenceIndex());
		mapped.setMateAlignmentStart(unmapped.getAlignmentStart());
		mapped.setMateNegativeStrandFlag(unmapped.getReadNegativeStrandFlag());
		mapped.setMateUnmappedFlag(true);
		mapped.setInferredInsertSize(0);

		unmapped.setMateReferenceIndex(mapped.getReferenceIndex());
		unmapped.setMateAlignmentStart(mapped.getAlignmentStart());
		unmapped.setMateNegativeStrandFlag(mapped.getReadNegativeStrandFlag());
		unmapped.setMateUnmappedFlag(false);
		unmapped.setInferredInsertSize(0);
	}

	boolean firstIsFirst = rec1.getAlignmentStart() < rec2.getAlignmentStart();
	int insertSize = firstIsFirst ? SamPairUtil.computeInsertSize(rec1, rec2) : SamPairUtil.computeInsertSize(rec2,
			rec1);

	rec1.setInferredInsertSize(firstIsFirst ? insertSize : -insertSize);
	rec2.setInferredInsertSize(firstIsFirst ? -insertSize : insertSize);

}
 
Example #13
Source File: InsertSizeMetricsCollectorArgs.java    From gatk with BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
public InsertSizeMetricsCollectorArgs(final int insertSize, final SamPairUtil.PairOrientation po) {
    this.insertSize = insertSize;
    this.po = po;
}
 
Example #14
Source File: InsertSizeMetricsCollectorArgs.java    From gatk with BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
public SamPairUtil.PairOrientation getPairOrientation() {
    return po;
}
 
Example #15
Source File: MergeBamAlignmentTest.java    From picard with MIT License 4 votes vote down vote up
@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.");
            }
        }
    }
}
 
Example #16
Source File: MergeBamAlignmentTest.java    From picard with MIT License 4 votes vote down vote up
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));
    }
}
 
Example #17
Source File: MergeBamAlignmentTest.java    From picard with MIT License 4 votes vote down vote up
/**
 * Confirm that paired reads are rejected by PrimaryAlignmentStrategy.EarliestFragment.
 */
@Test(expectedExceptions = UnsupportedOperationException.class)
public void testEarliestFragmentStrategyPaired() throws Exception {
    final File output = File.createTempFile("mergeTest", ".sam");
    output.deleteOnExit();

    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 cigar = "16M";

    final SAMRecord firstOfPair = new SAMRecord(header);
    firstOfPair.setReadName("theRead");
    firstOfPair.setReadString("ACGTACGTACGTACGT");
    firstOfPair.setBaseQualityString("5555555555555555");
    firstOfPair.setReadUnmappedFlag(true);
    firstOfPair.setReadPairedFlag(true);
    firstOfPair.setFirstOfPairFlag(true);

    final SAMRecord secondOfPair = new SAMRecord(header);
    secondOfPair.setReadName("theRead");
    secondOfPair.setReadString("ACGTACGTACGTACGT");
    secondOfPair.setBaseQualityString("5555555555555555");
    secondOfPair.setReadUnmappedFlag(true);
    secondOfPair.setReadPairedFlag(true);
    secondOfPair.setSecondOfPairFlag(true);
    SamPairUtil.setMateInfo(firstOfPair, secondOfPair);

    final SAMFileWriter unmappedWriter = factory.makeSAMWriter(header, false, unmappedSam);
    unmappedWriter.addAlignment(firstOfPair);
    unmappedWriter.addAlignment(secondOfPair);
    unmappedWriter.close();

    final File alignedSam = File.createTempFile("aligned.", ".sam");
    alignedSam.deleteOnExit();

    // Populate the header with SAMSequenceRecords
    header.setSequenceDictionary(SAMSequenceDictionaryExtractor.extractDictionary(sequenceDict2.toPath()));

    // Create 2 alignments for each end of pair
    final SAMFileWriter alignedWriter = factory.makeSAMWriter(header, false, alignedSam);
    for (int i = 1; i <= 2; ++i) {
        final SAMRecord firstOfPairAligned = new SAMRecord(header);
        firstOfPairAligned.setReadName(firstOfPair.getReadName());
        firstOfPairAligned.setReadBases(firstOfPair.getReadBases());
        firstOfPairAligned.setBaseQualities(firstOfPair.getBaseQualities());
        firstOfPairAligned.setReferenceName("chr1");
        firstOfPairAligned.setAlignmentStart(i);
        firstOfPairAligned.setCigarString(cigar);
        firstOfPairAligned.setMappingQuality(100);
        firstOfPairAligned.setReadPairedFlag(true);
        firstOfPairAligned.setFirstOfPairFlag(true);
        firstOfPairAligned.setAttribute(SAMTag.HI.name(), i);

        final SAMRecord secondOfPairAligned = new SAMRecord(header);
        secondOfPairAligned.setReadName(secondOfPair.getReadName());
        secondOfPairAligned.setReadBases(secondOfPair.getReadBases());
        secondOfPairAligned.setBaseQualities(secondOfPair.getBaseQualities());
        secondOfPairAligned.setReferenceName("chr1");
        secondOfPairAligned.setAlignmentStart(i + 10);
        secondOfPairAligned.setCigarString(cigar);
        secondOfPairAligned.setMappingQuality(100);
        secondOfPairAligned.setReadPairedFlag(true);
        secondOfPairAligned.setSecondOfPairFlag(true);
        secondOfPairAligned.setAttribute(SAMTag.HI.name(), i);

        SamPairUtil.setMateInfo(firstOfPairAligned, secondOfPairAligned);

        alignedWriter.addAlignment(firstOfPairAligned);
        alignedWriter.addAlignment(secondOfPairAligned);
    }
    alignedWriter.close();

    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.EarliestFragment,
            null, null, null, null);

    Assert.fail("Exception was not thrown");
}
 
Example #18
Source File: MergeBamAlignmentTest.java    From picard with MIT License 4 votes vote down vote up
/**
 * 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);
}
 
Example #19
Source File: MergeBamAlignmentTest.java    From picard with MIT License 4 votes vote down vote up
/**
 * Minimal test of merging data from separate read 1 and read 2 alignments
 */
@Test(dataProvider="separateTrimmed")
public void testMergingFromSeparatedReadTrimmedAlignments(final File unmapped, final List<File> r1Align, final List<File> r2Align, final int r1Trim, final int r2Trim, final String testName) throws Exception {
     final File output = File.createTempFile("mergeMultipleAlignmentsTest", ".sam");
     output.deleteOnExit();

    doMergeAlignment(unmapped, null, r1Align, r2Align, r1Trim, r2Trim,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true, fasta, output,
            SamPairUtil.PairOrientation.FR, null, null, null, null, null);

    SamReaderFactory factory = SamReaderFactory.makeDefault();
    final SamReader result = factory.open(output);
     for (final SAMRecord sam : result) {
        // Get the alignment record
        final List<File> rFiles = sam.getFirstOfPairFlag() ? r1Align : r2Align;
        SAMRecord alignment = null;
        for (final File f : rFiles) {
            for (final SAMRecord tmp : factory.open(f)) {
                if (tmp.getReadName().equals(sam.getReadName())) {
                    alignment = tmp;
                    break;
                }
            }
            if (alignment != null) break;
        }

        final int trim = sam.getFirstOfPairFlag() ? r1Trim : r2Trim;
        final int notWrittenToFastq = sam.getReadLength() - (trim + alignment.getReadLength());
        final int beginning = sam.getReadNegativeStrandFlag() ? notWrittenToFastq : trim;
        final int end = sam.getReadNegativeStrandFlag() ? trim : notWrittenToFastq;

        if (!sam.getReadUnmappedFlag()) {
            final CigarElement firstMergedCigarElement = sam.getCigar().getCigarElement(0);
            final CigarElement lastMergedCigarElement = sam.getCigar().getCigarElement(sam.getCigar().getCigarElements().size() - 1);
            final CigarElement firstAlignedCigarElement = alignment.getCigar().getCigarElement(0);
            final CigarElement lastAlignedCigarElement = alignment.getCigar().getCigarElement(alignment.getCigar().getCigarElements().size() - 1);

            if (beginning > 0) {
                Assert.assertEquals(firstMergedCigarElement.getOperator(), CigarOperator.S, "First element is not a soft clip");
                Assert.assertEquals(firstMergedCigarElement.getLength(), beginning + ((firstAlignedCigarElement.getOperator() == CigarOperator.S) ? firstAlignedCigarElement.getLength() : 0));
            }
            if (end > 0) {
                Assert.assertEquals(lastMergedCigarElement.getOperator(), CigarOperator.S, "Last element is not a soft clip");
                Assert.assertEquals(lastMergedCigarElement.getLength(), end + ((lastAlignedCigarElement.getOperator() == CigarOperator.S) ? lastAlignedCigarElement.getLength() : 0));
            }
        }
    }

}
 
Example #20
Source File: MergeBamAlignmentTest.java    From picard with MIT License 4 votes vote down vote up
@Test
public void testMergerFromMultipleFiles() throws Exception {
    final File output = File.createTempFile("mergeTest", ".sam");
    output.deleteOnExit();

    doMergeAlignment(unmappedBam, Arrays.asList(oneHalfAlignedBam, otherHalfAlignedBam),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true, fasta, output,
            SamPairUtil.PairOrientation.FR, null, null, null, null, null);

    final SamReader result = SamReaderFactory.makeDefault().open(output);

    for (final SAMRecord sam : result) {
        // This tests that we clip both (a) when the adapter is marked in the unmapped BAM file and
        // (b) when the insert size is less than the read length
        if (sam.getReadName().equals("both_reads_align_clip_adapter") ||
                sam.getReadName().equals("both_reads_align_clip_marked")) {
            Assert.assertEquals(sam.getReferenceName(), "chr7");
            if (sam.getReadNegativeStrandFlag()) {
                Assert.assertEquals(sam.getCigarString(), "5S96M", "Incorrect CIGAR string for " +
                        sam.getReadName());
            } else {
                Assert.assertEquals(sam.getCigarString(), "96M5S", "Incorrect CIGAR string for " +
                        sam.getReadName());
            }
        }
        // This tests that we DON'T clip when we run off the end if there are equal to or more than
        // MIN_ADAPTER_BASES hanging off the end
        else if (sam.getReadName().equals("both_reads_align_min_adapter_bases_exceeded")) {
            Assert.assertEquals(sam.getReferenceName(), "chr7");
            Assert.assertTrue(!sam.getCigarString().contains("S"),
                    "Read was clipped when it should not be.");
        } else if (sam.getReadName().equals("neither_read_aligns_or_present")) {
            Assert.assertTrue(sam.getReadUnmappedFlag(), "Read should be unmapped but isn't");
        }
        // Two pairs in which only the first read should align
        else if (sam.getReadName().equals("both_reads_present_only_first_aligns") ||
                sam.getReadName().equals("read_2_too_many_gaps")) {
            if (sam.getFirstOfPairFlag()) {
                Assert.assertEquals(sam.getReferenceName(), "chr7", "Read should be mapped but isn't");
            } else {
                Assert.assertTrue(sam.getReadUnmappedFlag(), "Read should not be mapped but is");
            }
        } else {
            throw new Exception("Unexpected read name: " + sam.getReadName());
        }
    }
}
 
Example #21
Source File: MergeBamAlignmentTest.java    From picard with MIT License 4 votes vote down vote up
@Test
public void testMerger() throws Exception {
    final File output = File.createTempFile("mergeTest", ".sam");
    output.deleteOnExit();

    doMergeAlignment(unmappedBam, Collections.singletonList(alignedBam),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true, fasta, output,
            SamPairUtil.PairOrientation.FR, null, null, null, null, null);

    SamReader result = SamReaderFactory.makeDefault().open(output);
    Assert.assertEquals(result.getFileHeader().getSequenceDictionary().getSequences().size(), 8,
            "Number of sequences did not match");
    SAMProgramRecord pg = result.getFileHeader().getProgramRecords().get(0);
    Assert.assertEquals(pg.getProgramGroupId(), "0");
    Assert.assertEquals(pg.getProgramVersion(), "1.0");
    Assert.assertEquals(pg.getCommandLine(), "align!");
    Assert.assertEquals(pg.getProgramName(), "myAligner");

    final SAMReadGroupRecord rg = result.getFileHeader().getReadGroups().get(0);
    Assert.assertEquals(rg.getReadGroupId(), "0");
    Assert.assertEquals(rg.getSample(), "Hi,Mom!");

    Assert.assertEquals(result.getFileHeader().getSortOrder(), SAMFileHeader.SortOrder.coordinate);

    for (final SAMRecord sam : result) {
        // This tests that we clip both (a) when the adapter is marked in the unmapped BAM file and
        // (b) when the insert size is less than the read length
        if (sam.getReadName().equals("both_reads_align_clip_adapter") ||
                sam.getReadName().equals("both_reads_align_clip_marked")) {
            Assert.assertEquals(sam.getReferenceName(), "chr7");
            if (sam.getReadNegativeStrandFlag()) {
                Assert.assertEquals(sam.getCigarString(), "5S96M", "Incorrect CIGAR string for " +
                        sam.getReadName());
            } else {
                Assert.assertEquals(sam.getCigarString(), "96M5S", "Incorrect CIGAR string for " +
                        sam.getReadName());
            }
        }
        // This tests that we DON'T clip when we run off the end if there are equal to or more than
        // MIN_ADAPTER_BASES hanging off the end
        else if (sam.getReadName().equals("both_reads_align_min_adapter_bases_exceeded")) {
            Assert.assertEquals(sam.getReferenceName(), "chr7");
            Assert.assertTrue(sam.getCigarString().indexOf("S") == -1,
                    "Read was clipped when it should not be.");
        } else if (sam.getReadName().equals("neither_read_aligns_or_present")) {
            Assert.assertTrue(sam.getReadUnmappedFlag(), "Read should be unmapped but isn't");
        }
        // Two pairs in which only the first read should align
        else if (sam.getReadName().equals("both_reads_present_only_first_aligns") ||
                sam.getReadName().equals("read_2_too_many_gaps")) {
            if (sam.getFirstOfPairFlag()) {
                Assert.assertEquals(sam.getReferenceName(), "chr7", "Read should be mapped but isn't");
            } else {
                Assert.assertTrue(sam.getReadUnmappedFlag(), "Read should not be mapped but is");
            }
        } else {
            throw new Exception("Unexpected read name: " + sam.getReadName());
        }

        final boolean pos = !sam.getReadNegativeStrandFlag();
        Assert.assertEquals(sam.getAttribute("aa"), pos ? "Hello" : "olleH");
        Assert.assertEquals(sam.getAttribute("ab"), pos ? "ATTCGG" : "CCGAAT");
        Assert.assertEquals(sam.getAttribute("ac"), pos ? new byte[] {1,2,3} : new byte[] {3,2,1});
        Assert.assertEquals(sam.getAttribute("as"), pos ? new short[]{1,2,3} : new short[]{3,2,1});
        Assert.assertEquals(sam.getAttribute("ai"), pos ? new int[]  {1,2,3} : new int[]  {3,2,1});
        Assert.assertEquals(sam.getAttribute("af"), pos ? new float[]{1,2,3} : new float[]{3,2,1});
    }

    // Quick test to make sure the program record gets picked up from the file if not specified
    // on the command line.
    doMergeAlignment(unmappedBam, Collections.singletonList(alignedBam),
            null, null, null, null,
            false, true, false, 1,
            null, null, null, null,
            true, fasta, output,
            SamPairUtil.PairOrientation.FR, null, null, null, null, null);

    CloserUtil.close(result);

    result = SamReaderFactory.makeDefault().open(output);
    pg = result.getFileHeader().getProgramRecords().get(0);
    Assert.assertEquals(pg.getProgramGroupId(), "1",
            "Program group ID not picked up correctly from aligned BAM");
    Assert.assertEquals(pg.getProgramVersion(), "2.0",
            "Program version not picked up correctly from aligned BAM");
    Assert.assertNull(pg.getCommandLine(),
            "Program command line not picked up correctly from aligned BAM");
    Assert.assertEquals(pg.getProgramName(), "Hey!",
            "Program name not picked up correctly from aligned BAM");

    CloserUtil.close(result);
}
 
Example #22
Source File: MergeBamAlignmentTest.java    From picard with MIT License 4 votes vote down vote up
@Test
public void testMergerWithSupplemental() throws Exception {
    final File outputWithSupplemental = File.createTempFile("mergeWithSupplementalTest", ".sam");
    outputWithSupplemental.deleteOnExit();

    doMergeAlignment(unmappedBam,
            Collections.singletonList(supplementalReadAlignedBam),
            null, null, null, null,
            false, true, false, 1,
            "0", "1.0", "align!", "myAligner",
            true, fasta, outputWithSupplemental,
            SamPairUtil.PairOrientation.FR, null, null, null, null, null);

    final SamReader result = SamReaderFactory.makeDefault().open(outputWithSupplemental);

    final List<Integer> clipAdapterFlags = new ArrayList<Integer>(Arrays.asList(99, 2147, 147, 2195));
    final List<Integer> foundClipAdapterFlags = new ArrayList<Integer>();

    for (final SAMRecord sam : result) {
        if (sam.getReadName().equals("both_reads_align_clip_adapter")) {
            foundClipAdapterFlags.add(sam.getFlags());
        }

        // This tests that we clip both (a) when the adapter is marked in the unmapped BAM file and
        // (b) when the insert size is less than the read length
        if (sam.getReadName().equals("both_reads_align_clip_marked")) {
            Assert.assertEquals(sam.getReferenceName(), "chr7");
            if (sam.getReadNegativeStrandFlag()) {
                Assert.assertEquals(sam.getCigarString(), "5S96M", "Incorrect CIGAR string for " + sam.getReadName());
            } else {
                Assert.assertEquals(sam.getCigarString(), "96M5S", "Incorrect CIGAR string for " + sam.getReadName());
            }
        }
        else if (sam.getReadName().equals("both_reads_align_clip_adapter")) {
            Assert.assertEquals(sam.getReferenceName(), "chr7");
            if (!sam.getSupplementaryAlignmentFlag()) {
                if (sam.getReadNegativeStrandFlag()) Assert.assertEquals(sam.getCigarString(), "5S96M", "Incorrect CIGAR string for " + sam.getReadName());
                else Assert.assertEquals(sam.getCigarString(), "96M5S", "Incorrect CIGAR string for " + sam.getReadName());
            }
        }
        // This tests that we DON'T clip when we run off the end if there are equal to or more than
        // MIN_ADAPTER_BASES hanging off the end
        else if (sam.getReadName().equals("both_reads_align_min_adapter_bases_exceeded")) {
            Assert.assertEquals(sam.getReferenceName(), "chr7");
            Assert.assertTrue(!sam.getCigarString().contains("S"),
                    "Read was clipped when it should not be.");
        } else if (sam.getReadName().equals("neither_read_aligns_or_present")) {
            Assert.assertTrue(sam.getReadUnmappedFlag(), "Read should be unmapped but isn't");
        }
        // Two pairs in which only the first read should align
        else if (sam.getReadName().equals("both_reads_present_only_first_aligns") ||
                sam.getReadName().equals("read_2_too_many_gaps")) {
            if (sam.getFirstOfPairFlag()) {
                Assert.assertEquals(sam.getReferenceName(), "chr7", "Read should be mapped but isn't");
            } else {
                Assert.assertTrue(sam.getReadUnmappedFlag(), "Read should not be mapped but is");
            }
        } else {
            throw new Exception("Unexpected read name: " + sam.getReadName());
        }
    }

    // Make sure that we have the appropriate primary and supplementary reads in the new file
    Assert.assertEquals(foundClipAdapterFlags.size(), clipAdapterFlags.size());
    Collections.sort(clipAdapterFlags);
    Collections.sort(foundClipAdapterFlags);
    for (int i = 0; i < clipAdapterFlags.size(); i++) {
        Assert.assertEquals(clipAdapterFlags.get(i), foundClipAdapterFlags.get(i));
    }

}
 
Example #23
Source File: ChimeraUtil.java    From picard with MIT License 4 votes vote down vote up
private static boolean matchesExpectedOrientations(final SAMRecord rec, final Set<PairOrientation> expectedOrientations) {
    return expectedOrientations.contains(SamPairUtil.getPairOrientation(rec)) && rec.getAttribute(SAMTag.SA.toString()) == null;
}
 
Example #24
Source File: InsertSizeMetricsCollector.java    From picard with MIT License 4 votes vote down vote up
public InsertSizeCollectorArgs(final int insertSize, final SamPairUtil.PairOrientation po) {
    this.insertSize = insertSize;
    this.po = po;
}
 
Example #25
Source File: InsertSizeMetricsCollector.java    From picard with MIT License 4 votes vote down vote up
public SamPairUtil.PairOrientation getPairOrientation() {
    return po;
}
 
Example #26
Source File: InsertSizeMetricsCollector.java    From picard with MIT License 4 votes vote down vote up
public void addMetricsToFile(final MetricsFile<InsertSizeMetrics,Integer> file) {
    // get the number of inserts, and the maximum and minimum keys across, across all orientations
    for (final Histogram<Integer> h : this.histograms.values()) {
        totalInserts += h.getCount();
    }
    if (0 == totalInserts) return; // nothing to store

    for(final Map.Entry<SamPairUtil.PairOrientation, Histogram<Integer>> entry : histograms.entrySet()) {
        final SamPairUtil.PairOrientation pairOrientation = entry.getKey();
        final Histogram<Integer> histogram = entry.getValue();
        final double total = histogram.getCount();

        // Only include a category if it has a sufficient percentage of the data in it
        if( total >= totalInserts * minimumPct ) {
            final InsertSizeMetrics metrics = new InsertSizeMetrics();
            metrics.SAMPLE             = this.sample;
            metrics.LIBRARY            = this.library;
            metrics.READ_GROUP         = this.readGroup;
            metrics.PAIR_ORIENTATION   = pairOrientation;
            if (!histogram.isEmpty()) {
                metrics.READ_PAIRS = (long) total;
                metrics.MAX_INSERT_SIZE = (int) histogram.getMax();
                metrics.MIN_INSERT_SIZE = (int) histogram.getMin();
                metrics.MEDIAN_INSERT_SIZE = histogram.getMedian();
                metrics.MODE_INSERT_SIZE = histogram.getMode();
                metrics.MEDIAN_ABSOLUTE_DEVIATION = histogram.getMedianAbsoluteDeviation();

                final double median = histogram.getMedian();
                double covered = 0;
                double low = median;
                double high = median;

                while (low >= histogram.getMin()-1 || high <= histogram.getMax()+1) {
                    final Histogram.Bin<Integer> lowBin = histogram.get((int) low);
                    if (lowBin != null) covered += lowBin.getValue();

                    if (low != high) {
                        final Histogram.Bin<Integer> highBin = histogram.get((int) high);
                        if (highBin != null) covered += highBin.getValue();
                    }

                    final double percentCovered = covered / total;
                    final int distance = (int) (high - low) + 1;
                    if (percentCovered >= 0.1 && metrics.WIDTH_OF_10_PERCENT == 0) metrics.WIDTH_OF_10_PERCENT = distance;
                    if (percentCovered >= 0.2 && metrics.WIDTH_OF_20_PERCENT == 0) metrics.WIDTH_OF_20_PERCENT = distance;
                    if (percentCovered >= 0.3 && metrics.WIDTH_OF_30_PERCENT == 0) metrics.WIDTH_OF_30_PERCENT = distance;
                    if (percentCovered >= 0.4 && metrics.WIDTH_OF_40_PERCENT == 0) metrics.WIDTH_OF_40_PERCENT = distance;
                    if (percentCovered >= 0.5 && metrics.WIDTH_OF_50_PERCENT == 0) metrics.WIDTH_OF_50_PERCENT = distance;
                    if (percentCovered >= 0.6 && metrics.WIDTH_OF_60_PERCENT == 0) metrics.WIDTH_OF_60_PERCENT = distance;
                    if (percentCovered >= 0.7 && metrics.WIDTH_OF_70_PERCENT == 0) metrics.WIDTH_OF_70_PERCENT = distance;
                    if (percentCovered >= 0.8 && metrics.WIDTH_OF_80_PERCENT == 0) metrics.WIDTH_OF_80_PERCENT = distance;
                    if (percentCovered >= 0.9 && metrics.WIDTH_OF_90_PERCENT == 0) metrics.WIDTH_OF_90_PERCENT = distance;
                    if (percentCovered >= 0.95 && metrics.WIDTH_OF_95_PERCENT == 0) metrics.WIDTH_OF_95_PERCENT = distance;
                    if (percentCovered >= 0.99 && metrics.WIDTH_OF_99_PERCENT == 0) metrics.WIDTH_OF_99_PERCENT = distance;

                    --low;
                    ++high;
                }
            }

            // Trim the Histogram down to get rid of outliers that would make the chart useless.
            final Histogram<Integer> trimmedHistogram = histogram; // alias it
            trimmedHistogram.trimByWidth(getWidthToTrimTo(metrics));

            if (!trimmedHistogram.isEmpty()) {
                metrics.MEAN_INSERT_SIZE = trimmedHistogram.getMean();
                metrics.STANDARD_DEVIATION = trimmedHistogram.getStandardDeviation();
            }

            file.addHistogram(trimmedHistogram);
            file.addMetric(metrics);
        }
    }
}
 
Example #27
Source File: Read.java    From cramtools with Apache License 2.0 3 votes vote down vote up
SAMRecord[] toSAMRecord(SAMFileHeader header) {

		SAMRecord first = firstSAMRecord(header);
		SAMRecord second = secondSAMRecord(header);

		SamPairUtil.setMateInfo(first, second, header);

		first.setProperPairFlag(true);
		second.setProperPairFlag(true);

		return new SAMRecord[] { first, second };
	}