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

The following examples show how to use htsjdk.samtools.SAMRecord#setReadName() . 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: Read.java    From cramtools with Apache License 2.0 6 votes vote down vote up
SAMRecord secondSAMRecord(SAMFileHeader header) {
	SAMRecord r = new SAMRecord(header);
	r.setReadName(evidenceRecord.getReadName());
	r.setReferenceName(evidenceRecord.Chromosome);
	r.setAlignmentStart(Integer.valueOf(evidenceRecord.MateOffsetInReference) + 1);
	r.setMappingQuality(Integer.valueOf(evidenceRecord.ScoreAllele0));
	r.setReadPairedFlag(true);
	r.setReadUnmappedFlag(false);
	r.setReadNegativeStrandFlag(negative);
	r.setFirstOfPairFlag(evidenceRecord.side == 1);
	r.setSecondOfPairFlag(!r.getFirstOfPairFlag());

	r.setCigar(new Cigar(Utils.toCigarOperatorList(secondHalf.samCigarElements)));

	r.setReadBases(Utils.toByteArray(secondHalf.readBasesBuf));
	r.setBaseQualities(Utils.toByteArray(secondHalf.readScoresBuf));

	r.setAttribute("GC", Utils.toString(secondHalf.gcList));
	r.setAttribute("GS", Utils.toString(secondHalf.gsBuf));
	r.setAttribute("GQ", SAMUtils.phredToFastq(Utils.toByteArray(secondHalf.gqBuf)));

	return r;
}
 
Example 2
Source File: MergeBamAlignmentTest.java    From picard with MIT License 6 votes vote down vote up
private void addAlignmentForMostStrategy(
        final SAMFileWriter writer, final SAMRecord unmappedRecord, final MostDistantStrategyAlignmentSpec spec,
        final boolean reverse) {
    final SAMRecord alignedRecord = new SAMRecord(writer.getFileHeader());
    alignedRecord.setReadName(unmappedRecord.getReadName());
    alignedRecord.setReadBases(unmappedRecord.getReadBases());
    alignedRecord.setBaseQualities(unmappedRecord.getBaseQualities());
    alignedRecord.setReferenceName(spec.sequence);
    alignedRecord.setAlignmentStart(spec.alignmentStart);
    alignedRecord.setReadNegativeStrandFlag(reverse);
    alignedRecord.setCigarString(unmappedRecord.getReadBases().length + "M");
    alignedRecord.setMappingQuality(spec.mapQ);
    alignedRecord.setReadPairedFlag(unmappedRecord.getReadPairedFlag());
    alignedRecord.setFirstOfPairFlag(unmappedRecord.getFirstOfPairFlag());
    alignedRecord.setSecondOfPairFlag(unmappedRecord.getSecondOfPairFlag());
    alignedRecord.setMateUnmappedFlag(true);
    writer.addAlignment(alignedRecord);
}
 
Example 3
Source File: MergeBamAlignmentTest.java    From picard with MIT License 6 votes vote down vote up
private void addAlignmentsForBestFragmentMapqStrategy(
        final SAMFileWriter writer, final SAMRecord unmappedRecord, final String sequence, final int[] mapqs) {
    boolean reverse = false;
    int alignmentStart = 1;
    for (final int mapq : mapqs) {
        final SAMRecord alignedRecord = new SAMRecord(writer.getFileHeader());
        alignedRecord.setReadName(unmappedRecord.getReadName());
        alignedRecord.setReadBases(unmappedRecord.getReadBases());
        alignedRecord.setBaseQualities(unmappedRecord.getBaseQualities());
        alignedRecord.setReferenceName(sequence);
        alignedRecord.setAlignmentStart(alignmentStart);
        alignmentStart += 10; // Any old position will do
        alignedRecord.setReadNegativeStrandFlag(reverse);
        reverse = !reverse;
        alignedRecord.setCigarString(unmappedRecord.getReadBases().length + "M");
        alignedRecord.setMappingQuality(mapq);
        alignedRecord.setReadPairedFlag(unmappedRecord.getReadPairedFlag());
        alignedRecord.setFirstOfPairFlag(unmappedRecord.getFirstOfPairFlag());
        alignedRecord.setSecondOfPairFlag(unmappedRecord.getSecondOfPairFlag());
        alignedRecord.setMateUnmappedFlag(true);
        writer.addAlignment(alignedRecord);
    }
}
 
Example 4
Source File: FastqToSam.java    From picard with MIT License 6 votes vote down vote up
private SAMRecord createSamRecord(final SAMFileHeader header, final String baseName, final FastqRecord frec, final boolean paired) {
    final SAMRecord srec = new SAMRecord(header);
    srec.setReadName(baseName);
    srec.setReadString(frec.getReadString());
    srec.setReadUnmappedFlag(true);
    srec.setAttribute(ReservedTagConstants.READ_GROUP_ID, READ_GROUP_NAME);
    final byte[] quals = StringUtil.stringToBytes(frec.getBaseQualityString());
    convertQuality(quals, QUALITY_FORMAT);
    for (final byte qual : quals) {
        final int uQual = qual & 0xff;
        if (uQual < MIN_Q || uQual > MAX_Q) {
            throw new PicardException("Base quality " + uQual + " is not in the range " + MIN_Q + ".." +
            MAX_Q + " for read " + frec.getReadHeader());
        }
    }
    srec.setBaseQualities(quals);

    if (paired) {
        srec.setReadPairedFlag(true);
        srec.setMateUnmappedFlag(true);
    }
    return srec ;
}
 
Example 5
Source File: Read.java    From cramtools with Apache License 2.0 6 votes vote down vote up
SAMRecord firstSAMRecord(SAMFileHeader header) {
	SAMRecord r = new SAMRecord(header);
	r.setReadName(evidenceRecord.getReadName());
	r.setReferenceName(evidenceRecord.Chromosome);
	r.setAlignmentStart(Integer.valueOf(evidenceRecord.OffsetInReference) + 1);
	r.setMappingQuality(Integer.valueOf(evidenceRecord.ScoreAllele0));
	r.setReadPairedFlag(true);
	r.setReadUnmappedFlag(false);
	r.setReadNegativeStrandFlag(negative);
	r.setFirstOfPairFlag(evidenceRecord.side == 0);
	r.setSecondOfPairFlag(!r.getFirstOfPairFlag());

	r.setCigar(new Cigar(Utils.toCigarOperatorList(firstHalf.samCigarElements)));

	r.setReadBases(Utils.toByteArray(firstHalf.readBasesBuf));
	r.setBaseQualities(Utils.toByteArray(firstHalf.readScoresBuf));

	r.setAttribute("GC", Utils.toString(firstHalf.gcList));
	r.setAttribute("GS", Utils.toString(firstHalf.gsBuf));
	r.setAttribute("GQ", SAMUtils.phredToFastq(Utils.toByteArray(firstHalf.gqBuf)));

	return r;
}
 
Example 6
Source File: SamWriterWrapper.java    From rtg-tools with BSD 2-Clause "Simplified" License 6 votes vote down vote up
void writeSequence(SequencesReader reader, long seqId, byte[] dataBuffer, byte[] qualityBuffer, int flags) throws IOException {
  final SAMRecord rec = new SAMRecord(mWriter.getFileHeader());

  final int length = reader.read(seqId, dataBuffer);
  rec.setReadName(mHasNames ? reader.name(seqId) : String.valueOf(seqId));
  rec.setReferenceName("*");
  rec.setFlags(flags);
  rec.setReadBases(DnaUtils.bytesToSequenceIncCG(dataBuffer, 0, length).getBytes());
  if (mReader.hasQualityData()) {
    reader.readQuality(seqId, qualityBuffer);
    rec.setBaseQualities(Arrays.copyOf(qualityBuffer, length));
  }
  rec.setInferredInsertSize(0);
  rec.setMateAlignmentStart(0);
  rec.setMateReferenceName("*");
  if (mReadGroupRecord != null) {
    rec.setAttribute(ReadGroupUtils.RG_ATTRIBUTE, mReadGroupRecord.getReadGroupId());
  }

  mWriter.addAlignment(rec);
}
 
Example 7
Source File: BamToBfqWriter.java    From picard with MIT License 5 votes vote down vote up
/**
 * Path for writing bfqs for single-end reads
 *
 * @param iterator  the iterator with he SAM Records to write
 * @param filters   the list of filters to be applied
 */
private void writeSingleEndBfqs(final Iterator<SAMRecord> iterator, final List<SamRecordFilter> filters) {

    // Open the codecs for writing
    int fileIndex = 0;
    initializeNextBfqFiles(fileIndex++);

    int records = 0;

    final FilteringSamIterator it = new FilteringSamIterator(iterator, new AggregateFilter(filters));
    while (it.hasNext()) {
        final SAMRecord record = it.next();
        records++;
        if (records % increment == 0) {

            record.setReadName(record.getReadName() + "/1");
            writeFastqRecord(codec1, record);
            wrote++;
            if (wrote % 1000000 == 0) {
                log.info(wrote + " records processed.");
            }
            if (chunk > 0 && wrote % chunk == 0) {
                initializeNextBfqFiles(fileIndex++);
            }
        }
    }
}
 
Example 8
Source File: BAMRecordViewTest.java    From cramtools with Apache License 2.0 5 votes vote down vote up
@Test
public void test() {
	SAMFileHeader header = new SAMFileHeader();
	SAMRecord r1 = new SAMRecord(header);
	r1.setReadName("readName");
	r1.setFlags(4);
	r1.setReferenceIndex(SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX);
	r1.setAlignmentStart(SAMRecord.NO_ALIGNMENT_START);
	r1.setMappingQuality(SAMRecord.NO_MAPPING_QUALITY);
	r1.setCigar(new Cigar());
	r1.setMateReferenceIndex(SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX);
	r1.setMateAlignmentStart(SAMRecord.NO_ALIGNMENT_START);
	r1.setReadBases("A".getBytes());
	r1.setBaseQualityString("!");

	BAMRecordView view = new BAMRecordView(new byte[1024]);
	translate(r1, view);
	r1.setReadName("2");
	translate(r1, view);

	List<SAMRecord> list = toSAMRecord(view, header);
	assertEquals(2, list.size());

	Iterator<SAMRecord> iterator = list.iterator();
	SAMRecord r2 = iterator.next();
	r1.setReadName("readName");
	compare(r1, r2);

	r1.setReadName("2");
	r2 = iterator.next();
	compare(r1, r2);
}
 
Example 9
Source File: SamOutputTest.java    From rtg-tools with BSD 2-Clause "Simplified" License 5 votes vote down vote up
public void test() throws IOException {
  final SAMFileHeader header = new SAMFileHeader();
  header.setSortOrder(SAMFileHeader.SortOrder.coordinate);
  header.getSequenceDictionary().addSequence(new SAMSequenceRecord("seq", 1000000));
  final SAMRecord rec = new SAMRecord(header);
  rec.setAlignmentStart(500000);
  rec.setReferenceName("seq");
  rec.setReadName("read");
  rec.setReadString("acgt");
  rec.setBaseQualityString("####");
  rec.setCigarString("4=");
  try (TestDirectory dir = new TestDirectory()) {
    final MemoryPrintStream mps = new MemoryPrintStream();
    try (SamOutput out = SamOutput.getSamOutput(new File(dir, "foo"), mps.outputStream(), header, true, true, null)) {
      out.getWriter().addAlignment(rec);
    }
    mNano.check("samoutput_expected_1.sam", MainResult.run(new ExtractCli(), "--header", new File(dir, "foo.bam").getPath(), "seq:500000+1").out());
    assertTrue(new File(dir, "foo.bam.bai").exists());
    assertEquals("", mps.toString());

    try (SamOutput out = SamOutput.getSamOutput(new File(dir, "foo.sam"), mps.outputStream(), header, true, true, null)) {
      out.getWriter().addAlignment(rec);
    }
    mNano.check("samoutput_expected_1.sam", MainResult.run(new ExtractCli(), "--header", new File(dir, "foo.sam.gz").getPath(), "seq:500000+1").out());
    assertTrue(new File(dir, "foo.sam.gz.tbi").exists());
    assertEquals("", mps.toString());

    try (SamOutput out = SamOutput.getSamOutput(new File(dir, "foo.sam"), mps.outputStream(), header, false, true, null)) {
      out.getWriter().addAlignment(rec);
    }
    mNano.check("samoutput_expected_1.sam", FileHelper.fileToString(new File(dir, "foo.sam")));
    assertFalse(new File(dir, "foo.sam.tbi").exists());
    assertEquals("", mps.toString());

    try (SamOutput out = SamOutput.getSamOutput(new File("-"), mps.outputStream(), header, true, true, null)) {
      out.getWriter().addAlignment(rec);
    }
    mNano.check("samoutput_expected_1.sam", mps.toString());
  }
}
 
Example 10
Source File: MergeBamAlignmentTest.java    From picard with MIT License 5 votes vote down vote up
private SAMRecord makeRead(final SAMFileHeader alignedHeader, final SAMRecord unmappedRec, final HitSpec hitSpec, final int hitIndex) {
    final SAMRecord rec = new SAMRecord(alignedHeader);
    rec.setReadName(unmappedRec.getReadName());
    rec.setReadBases(unmappedRec.getReadBases());
    rec.setBaseQualities(unmappedRec.getBaseQualities());
    rec.setMappingQuality(hitSpec.mapq);
    if (!hitSpec.primary) rec.setNotPrimaryAlignmentFlag(true);
    final Cigar cigar = new Cigar();
    final int readLength = rec.getReadLength();
    if (hitSpec.filtered) {
        // Add two insertions so alignment is filtered.
        cigar.add(new CigarElement(readLength-4, CigarOperator.M));
        cigar.add(new CigarElement(1, CigarOperator.I));
        cigar.add(new CigarElement(1, CigarOperator.M));
        cigar.add(new CigarElement(1, CigarOperator.I));
        cigar.add(new CigarElement(1, CigarOperator.M));
    } else {
        cigar.add(new CigarElement(readLength, CigarOperator.M));
    }
    rec.setCigar(cigar);

    rec.setReferenceName(bigSequenceName);
    rec.setAttribute(SAMTag.HI.name(), hitIndex);
    rec.setAlignmentStart(hitIndex + 1);

    return rec;
}
 
Example 11
Source File: FastqRead.java    From cramtools with Apache License 2.0 5 votes vote down vote up
public SAMRecord toSAMRecord(SAMFileHeader header) {
	SAMRecord record = new SAMRecord(header);
	String name = null;
	if (data[nameLen - 1] == '/' && Character.isDigit(data[nameLen]))
		name = new String(data, 1, nameLen - 2);
	else
		name = new String(data, 1, nameLen - 2);
	record.setReadName(name);
	int readLen = (data.length - this.nameLen - 4 - 1) / 2;
	byte[] bases = Arrays.copyOfRange(data, nameLen + 2, nameLen + 2 + readLen);
	record.setReadBases(bases);

	byte[] scores = Arrays.copyOfRange(data, nameLen + 3 + 1 + readLen + 1, nameLen + 3 + 1 + 2 * readLen + 1);
	record.setBaseQualityString(new String(scores));

	record.setReadUnmappedFlag(true);
	switch (templateIndex) {
	case 0:
		record.setReadPairedFlag(false);
		break;
	case 1:
		record.setReadPairedFlag(true);
		record.setFirstOfPairFlag(true);
		break;
	case 2:
		record.setReadPairedFlag(true);
		record.setSecondOfPairFlag(true);
		break;

	default:
		break;
	}
	return record;
}
 
Example 12
Source File: FastWgsMetricsCollectorTest.java    From picard with MIT License 5 votes vote down vote up
@BeforeTest
public void setUp(){
    final byte[] referenceBases = ">chrM\nACCTACGTTCAATATTCTTCACCTACGTTCAATATTCTTCACCTACGTTCAATATTCTTCACCTACGTTCAATATTCTTCACCTACGTTCAATATTCTTC".getBytes();
    ref = new ReferenceSequence("chrM", 0, referenceBases);
    sequence = new SAMSequenceRecord("chrM", 100);
    record = new SAMRecord(new SAMFileHeader());
    record.setReadName("test");
    record.setBaseQualities(qualities);
    record.setReadBases(readBases);
    secondRecord = generateRecord("test1");
    thirdRecord = generateRecord("test2");
}
 
Example 13
Source File: FastWgsMetricsCollectorTest.java    From picard with MIT License 5 votes vote down vote up
private SAMRecord generateRecord(String name) {
    SAMRecord record = new SAMRecord(new SAMFileHeader());
    record.setReadName(name);
    record.setBaseQualities(highQualities);
    record.setReadBases(readBases);
    return record;
}
 
Example 14
Source File: HaplotypeBAMWriter.java    From gatk with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
/**
 * Write out a representation of this haplotype as a read
 *
 * @param haplotype a haplotype to write out, must not be null
 * @param paddedRefLoc the reference location, must not be null
 * @param isAmongBestHaplotypes true if among the best haplotypes, false if it was just one possible haplotype
 * @param callableRegion the region over which variants are being called
 */
private void writeHaplotype(final Haplotype haplotype,
                            final Locatable paddedRefLoc,
                            final boolean isAmongBestHaplotypes,
                            final Locatable callableRegion) {
    Utils.nonNull(haplotype, "haplotype cannot be null");
    Utils.nonNull(paddedRefLoc, "paddedRefLoc cannot be null");

    final SAMRecord record = new SAMRecord(output.getBAMOutputHeader());
    record.setReadBases(haplotype.getBases());
    record.setAlignmentStart(paddedRefLoc.getStart() + haplotype.getAlignmentStartHapwrtRef());
    // Use a base quality value "!" for it's display value (quality value is not meaningful)
    record.setBaseQualities(Utils.dupBytes((byte) '!', haplotype.getBases().length));
    record.setCigar(haplotype.getCigar());
    record.setMappingQuality(isAmongBestHaplotypes ? bestHaplotypeMQ : otherMQ);
    record.setReadName(output.getHaplotypeSampleTag() + uniqueNameCounter++);
    record.setAttribute(output.getHaplotypeSampleTag(), haplotype.hashCode());
    record.setReadUnmappedFlag(false);
    record.setReferenceIndex(output.getBAMOutputHeader().getSequenceIndex(paddedRefLoc.getContig()));
    record.setAttribute(SAMTag.RG.toString(), output.getHaplotypeReadGroupID());
    record.setFlags(SAMFlag.READ_REVERSE_STRAND.intValue());
    if (callableRegion != null) {
        record.setAttribute(AssemblyBasedCallerUtils.CALLABLE_REGION_TAG, callableRegion.toString());
    }

    output.add(new SAMRecordToGATKReadAdapter(record));
}
 
Example 15
Source File: DuplicateSamFilterTest.java    From rtg-tools with BSD 2-Clause "Simplified" License 4 votes vote down vote up
public void testPairedEnd() {
  final DuplicateSamFilter dsf = new DuplicateSamFilter();

  assertFalse(dsf.acceptRecord(null));

  final SAMRecord rec1f = new SAMRecord(null);

  assertFalse(dsf.acceptRecord(rec1f));

  final String firstReadName = "IronSheik";
  final String secondReadName = "Sabu";

  rec1f.setReadName(firstReadName);
  rec1f.setReadPairedFlag(true);
  rec1f.setFirstOfPairFlag(true);

  assertTrue(dsf.acceptRecord(rec1f));

  rec1f.setSecondaryAlignment(true);

  final SAMRecord rec2f = new SAMRecord(null);
  rec2f.setReadPairedFlag(true);
  rec2f.setFirstOfPairFlag(true);
  rec2f.setReadName(secondReadName);
  rec2f.setSecondaryAlignment(true);

  final SAMRecord rec3fdup = new SAMRecord(null);
  rec3fdup.setReadPairedFlag(true);
  rec3fdup.setFirstOfPairFlag(true);
  rec3fdup.setReadName(firstReadName);
  rec3fdup.setSecondaryAlignment(true);

  final SAMRecord rec1s = new SAMRecord(null);
  rec1s.setReadName(firstReadName);
  rec1s.setReadPairedFlag(true);
  rec1s.setFirstOfPairFlag(false);
  rec1s.setSecondaryAlignment(true);

  final SAMRecord rec2s = new SAMRecord(null);
  rec2s.setReadPairedFlag(true);
  rec2s.setFirstOfPairFlag(false);
  rec2s.setReadName(secondReadName);
  rec2s.setSecondaryAlignment(true);

  final SAMRecord rec3sdup = new SAMRecord(null);
  rec3sdup.setReadPairedFlag(true);
  rec3sdup.setFirstOfPairFlag(false);
  rec3sdup.setReadName(firstReadName);
  rec3sdup.setSecondaryAlignment(true);


  assertTrue(dsf.acceptRecord(rec1f));
  assertTrue(dsf.acceptRecord(rec1s));
  assertTrue(dsf.acceptRecord(rec2f));
  assertTrue(dsf.acceptRecord(rec2s));
  assertFalse(dsf.acceptRecord(rec1f));
  assertFalse(dsf.acceptRecord(rec1s));
  assertFalse(dsf.acceptRecord(rec3fdup));
  assertFalse(dsf.acceptRecord(rec3sdup));
}
 
Example 16
Source File: BamToBfqWriter.java    From picard with MIT License 4 votes vote down vote up
/**
 * Path for writing bfqs for paired-end reads
 *
 * @param iterator      the iterator witht he SAM Records to write
 * @param tagFilter     the filter for noise reads
 * @param qualityFilter the filter for PF reads
 */
private void writePairedEndBfqs(final Iterator<SAMRecord> iterator, final TagFilter tagFilter,
                                final FailsVendorReadQualityFilter qualityFilter,
                                SamRecordFilter ... otherFilters) {
    // Open the codecs for writing
    int fileIndex = 0;
    initializeNextBfqFiles(fileIndex++);

    int records = 0;

    RECORD_LOOP: while (iterator.hasNext()) {
        final SAMRecord first = iterator.next();
        if (!iterator.hasNext()) {
            throw new PicardException("Mismatched number of records in " + this.bamFile.getAbsolutePath());
        }
        final SAMRecord second = iterator.next();
        if (!second.getReadName().equals(first.getReadName()) ||
            first.getFirstOfPairFlag() == second.getFirstOfPairFlag()) {
            throw new PicardException("Unmatched read pairs in " + this.bamFile.getAbsolutePath() +
                ": " + first.getReadName() + ", " + second.getReadName() + ".");
        }

        // If *both* are noise reads, filter them out
        if (tagFilter.filterOut(first) && tagFilter.filterOut(second))  {
            continue;
        }

        // If either fails to pass filter, then exclude them as well
        if (!includeNonPfReads && (qualityFilter.filterOut(first) || qualityFilter.filterOut(second))) {
            continue;
        }

        // If either fails any of the other filters, exclude them both
        for (SamRecordFilter filter : otherFilters) {
            if (filter.filterOut(first) || filter.filterOut(second)) {
                continue RECORD_LOOP;
            }
        }

        // Otherwise, write them out
        records++;
        if (records % increment == 0) {
            first.setReadName(first.getReadName() + "/1");
            writeFastqRecord(first.getFirstOfPairFlag() ? codec1 : codec2, first);
            second.setReadName(second.getReadName() + "/2");
            writeFastqRecord(second.getFirstOfPairFlag() ? codec1 : codec2, second);
            wrote++;
            if (wrote % 1000000 == 0) {
                log.info(wrote + " records written.");
            }
            if (chunk > 0 && wrote % chunk == 0) {
                initializeNextBfqFiles(fileIndex++);
            }
        }
    }
}
 
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: TestBAMRecordView.java    From cramtools with Apache License 2.0 4 votes vote down vote up
@Test
public void test() throws IOException {
	byte[] buf = new byte[1024];
	BAMRecordView view = new BAMRecordView(buf);
	view.setRefID(0);
	view.setAlignmentStart(77);
	view.setMappingScore(44);
	view.setIndexBin(99);
	view.setFlags(555);
	view.setMateRefID(0);
	view.setMateAlStart(78);
	view.setInsertSize(133);

	view.setReadName("name1");
	view.setCigar(TextCigarCodec.decode("10M"));
	view.setBases("AAAAAAAAAA".getBytes());
	view.setQualityScores("BBBBBBBBBB".getBytes());

	int id = 'A' << 16 | 'M' << 8 | 'A';
	view.addTag(id, "Q".getBytes(), 0, 1);

	int len = view.finish();

	System.out.println(Arrays.toString(Arrays.copyOf(buf, len)));

	ByteArrayOutputStream baos = new ByteArrayOutputStream();

	SAMFileHeader header = new SAMFileHeader();
	header.addSequence(new SAMSequenceRecord("14", 14));

	ByteArrayOutputStream baos2 = new ByteArrayOutputStream();
	SAMFileWriter writer = new SAMFileWriterFactory().makeBAMWriter(header, true, baos2);
	SAMRecord record = new SAMRecord(header);
	record.setReferenceIndex(0);
	record.setAlignmentStart(1);
	record.setCigarString("10M");
	record.setFlags(555);
	record.setMappingQuality(44);
	record.setMateReferenceIndex(0);
	record.setMateAlignmentStart(0);
	record.setInferredInsertSize(133);
	record.setReadName("name1");
	record.setReadBases("AAAAAAAAAA".getBytes());
	record.setBaseQualities("BBBBBBBBBB".getBytes());
	record.setAttribute("AM", 'Q');

	System.out.println("BAMFileWriter.addAlignment():");
	writer.addAlignment(record);
	System.out.println(".");
	writer.close();

	System.out.println("------------------------------------------");
	System.out.println();
	System.out.println(new String(baos2.toByteArray()));
	System.out.println();

	SAMFileReader.setDefaultValidationStringency(ValidationStringency.SILENT);
	SAMFileReader reader2 = new SAMFileReader(new ByteArrayInputStream(baos2.toByteArray()));
	SAMRecordIterator iterator = reader2.iterator();
	while (iterator.hasNext()) {
		record = iterator.next();
		System.out.println(record.getSAMString());
	}
	System.out.println("------------------------------------------");

	BlockCompressedOutputStream bcos = new BlockCompressedOutputStream(baos, null);
	bcos.write("BAM\1".getBytes());
	bcos.write(toByteArray(header));
	CramInt.writeInt32(header.getSequenceDictionary().size(), bcos);
	for (final SAMSequenceRecord sequenceRecord : header.getSequenceDictionary().getSequences()) {
		byte[] bytes = sequenceRecord.getSequenceName().getBytes();
		CramInt.writeInt32(bytes.length + 1, bcos);
		bcos.write(sequenceRecord.getSequenceName().getBytes());
		bcos.write(0);
		CramInt.writeInt32(sequenceRecord.getSequenceLength(), bcos);
	}
	bcos.write(buf, 0, len);
	bcos.close();

	System.out.println(new String(baos.toByteArray()));

	SAMFileReader reader = new SAMFileReader(new ByteArrayInputStream(baos.toByteArray()));
	iterator = reader.iterator();
	while (iterator.hasNext()) {
		record = iterator.next();
		System.out.println(record.getSAMString());
	}
	reader.close();

}
 
Example 19
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 20
Source File: DuplicateSamFilterTest.java    From rtg-tools with BSD 2-Clause "Simplified" License 3 votes vote down vote up
public void testSingleEnd() {
  final DuplicateSamFilter dsf = new DuplicateSamFilter();

  assertFalse(dsf.acceptRecord(null));

  final SAMRecord rec1 = new SAMRecord(null);

  assertFalse(dsf.acceptRecord(rec1));

  rec1.setReadName("IronSheik");

  assertTrue(dsf.acceptRecord(rec1));

  rec1.setSecondaryAlignment(true);

  final SAMRecord rec2 = new SAMRecord(null);
  rec2.setReadName("Sabu");
  rec2.setSecondaryAlignment(true);

  final SAMRecord rec3 = new SAMRecord(null);
  rec3.setReadName("IronSheik");
  rec3.setSecondaryAlignment(true);

  assertTrue(dsf.acceptRecord(rec1));
  assertTrue(dsf.acceptRecord(rec2));
  assertFalse(dsf.acceptRecord(rec1));
  assertFalse(dsf.acceptRecord(rec3));
}