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

The following examples show how to use htsjdk.samtools.SAMRecord#setAlignmentStart() . 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: SmartSamWriterTest.java    From rtg-tools with BSD 2-Clause "Simplified" License 6 votes vote down vote up
private static SAMRecord createRecord(int position, String ref) {
  final SAMFileHeader header = new SAMFileHeader();
  header.getSequenceDictionary().addSequence(new SAMSequenceRecord("A", 100000));
  header.getSequenceDictionary().addSequence(new SAMSequenceRecord("B", 100000));
  final SAMRecord rec = new SAMRecord(header);
  rec.setReferenceName(ref);
  rec.setReadName("read");
  rec.setReadString("ACGT");
  rec.setBaseQualityString("####");
  if (ref.equals(SAMRecord.NO_ALIGNMENT_REFERENCE_NAME)) {
    rec.setReadUnmappedFlag(true);
  } else {
    rec.setCigarString("4=");
    rec.setAlignmentStart(position);
  }
  return rec;
}
 
Example 3
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 4
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 5
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 6
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 7
Source File: UmiUtilTest.java    From picard with MIT License 5 votes vote down vote up
@Test(dataProvider = "topStrandDataProvider")
public void testIsTopStrand(final int referenceIndex, final int alignmentStart, final int mateReferenceIndex, final int mateAlignmentStart,
                            final boolean firstOfPairFlag, final boolean negativeStrandFlag, final boolean mateNegativeStrandFlag,
                            final boolean mapped, final boolean mateMapped,
                            final UmiUtil.ReadStrand strand) {

    final int readLength = 15;
    final int contigLength = 500;
    final SAMFileHeader header = new SAMFileHeader();
    final SAMSequenceDictionary sequenceDictionary = new SAMSequenceDictionary();

    sequenceDictionary.addSequence(new SAMSequenceRecord("chr1", contigLength));
    sequenceDictionary.addSequence(new SAMSequenceRecord("chr2", contigLength));

    header.setSequenceDictionary(sequenceDictionary);

    final SAMRecord rec = new SAMRecord(header);

    rec.setReadUnmappedFlag(!mapped);
    rec.setMateUnmappedFlag(!mateMapped);
    rec.setReadPairedFlag(true);

    rec.setCigarString(readLength + "M");
    rec.setAttribute("MC", readLength + "M");

    rec.setReferenceIndex(referenceIndex);
    rec.setAlignmentStart(alignmentStart);

    rec.setMateReferenceIndex(mateReferenceIndex);
    rec.setMateAlignmentStart(mateAlignmentStart);

    rec.setFirstOfPairFlag(firstOfPairFlag);
    rec.setReadNegativeStrandFlag(negativeStrandFlag);
    rec.setMateNegativeStrandFlag(mateNegativeStrandFlag);

    Assert.assertEquals(UmiUtil.getStrand(rec), strand);
}
 
Example 8
Source File: OrientationTest.java    From rtg-tools with BSD 2-Clause "Simplified" License 5 votes vote down vote up
public static SAMRecord makeTestRecord(SAMFileHeader sfh, int flags, String sequence, String mateSequence, int pos, int matePos, String readGroup) {
  final SAMRecord sr = new SAMRecord(sfh);
  sr.setFlags(flags);
  sr.setReferenceName(sequence);
  sr.setMateReferenceName(mateSequence);
  sr.setAlignmentStart(pos);
  sr.setMateAlignmentStart(matePos);
  sr.setAttribute("RG", readGroup);
  if (pos < matePos) {
    sr.setInferredInsertSize(matePos - pos + 5);
  } else {
    sr.setInferredInsertSize(matePos - pos  - 5); //-1 * (pos - matePos + 5));
  }
  return sr;
}
 
Example 9
Source File: SamCompareUtilsTest.java    From rtg-tools with BSD 2-Clause "Simplified" License 5 votes vote down vote up
public void test() {
  final SAMFileHeader header = new SAMFileHeader();
  header.setSequenceDictionary(new SAMSequenceDictionary(Arrays.asList(new SAMSequenceRecord("raga", 100), new SAMSequenceRecord("yaga", 100), new SAMSequenceRecord("zaga", 100))));
  final SAMRecord rec1 = new SAMRecord(header);
  rec1.setReferenceIndex(1);
  final SAMRecord rec2 = new SAMRecord(header);
  rec2.setReferenceIndex(2);
  assertEquals(-1, SamCompareUtils.compareSamRecords(rec1, rec2));
  assertEquals(1, SamCompareUtils.compareSamRecords(rec2, rec1));
  rec1.setReferenceIndex(2);
  rec1.setAlignmentStart(50);
  rec2.setAlignmentStart(25);
  assertEquals(1, SamCompareUtils.compareSamRecords(rec1, rec2));
  assertEquals(-1, SamCompareUtils.compareSamRecords(rec2, rec1));
  rec1.setReadPairedFlag(true);
  rec2.setReadPairedFlag(true);
  rec1.setProperPairFlag(true);
  rec2.setProperPairFlag(false);
  rec1.setAlignmentStart(25);
  assertEquals(-1, SamCompareUtils.compareSamRecords(rec1, rec2));
  assertEquals(1, SamCompareUtils.compareSamRecords(rec2, rec1));
  rec2.setProperPairFlag(true);
  rec1.setReadUnmappedFlag(true);
  assertEquals(1, SamCompareUtils.compareSamRecords(rec1, rec2));
  assertEquals(-1, SamCompareUtils.compareSamRecords(rec2, rec1));
  rec2.setReadUnmappedFlag(true);
  assertEquals(0, SamCompareUtils.compareSamRecords(rec1, rec2));
  assertEquals(0, SamCompareUtils.compareSamRecords(rec2, rec1));
  rec1.setReferenceIndex(-1);
  assertEquals(1, SamCompareUtils.compareSamRecords(rec1, rec2));
  assertEquals(-1, SamCompareUtils.compareSamRecords(rec2, rec1));
  rec2.setReferenceIndex(-1);
  assertEquals(0, SamCompareUtils.compareSamRecords(rec2, rec1));
}
 
Example 10
Source File: AnnotationUtilsTest.java    From Drop-seq with MIT License 5 votes vote down vote up
@Test
public void testGetConsistentExons () {
	// need to build a new read that falls into 2 genes to check if it's inconsistent.
	AnnotationUtils u = AnnotationUtils.getInstance();
	OverlapDetector<Gene> geneOverlapDetector = getGeneOverlapDetector(BAM, GTF);
	Map<String, Gene> map = getGeneMap(geneOverlapDetector);
	// expected UTR and coding.
	SAMRecord rec = getReadByName("000000000-AMY9M:1:2104:20987:12124", BAM);
	// change the cigar so the read maps to two genes exons.
	// rec.setCigar(cigar);
	rec.setAlignmentStart(6269430);
	Cigar c = new Cigar ();
	c.add(new CigarElement(20, CigarOperator.M));
	c.add(new CigarElement(9741377, CigarOperator.N));
	c.add(new CigarElement(30, CigarOperator.M));
	rec.setCigar(c);
	Set<Gene> genes = new HashSet<>();
	genes.add(map.get("RPL22"));
	genes.add(map.get("PLEKHM2"));

	// This read spans 2 genes, so it's not consistent under the strict test.
	Set<Gene> actual = u.getConsistentExons(rec, genes, false);
	Set<Gene> expected = new HashSet<>();
	Assert.assertEquals(actual, expected);

	// under the less strict option, we allow the read to span 2 different genes.
	actual = u.getConsistentExons(rec, genes, true);
	expected = genes;
	Assert.assertEquals(actual, expected);

}
 
Example 11
Source File: SAMRecordSerializerUnitTest.java    From gatk with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
@Test
public void testChangingContigsOnHeaderlessSAMRecord(){
    final SparkConf conf = new SparkConf().set("spark.kryo.registrator",
            "org.broadinstitute.hellbender.engine.spark.SAMRecordSerializerUnitTest$TestGATKRegistrator");
    final SAMRecord read = ((SAMRecordToGATKReadAdapter)ArtificialReadUtils.createHeaderlessSamBackedRead("read1", "1", 100, 50)).getEncapsulatedSamRecord();

    final SAMRecord roundTrippedRead = SparkTestUtils.roundTripInKryo(read, SAMRecord.class, conf);
    Assert.assertEquals(roundTrippedRead, read, "\nActual read: " + roundTrippedRead.getSAMString() + "\nExpected read: " + read.getSAMString());

    read.setReferenceName("2");
    read.setAlignmentStart(1);

    final SAMRecord roundTrippedRead2 = SparkTestUtils.roundTripInKryo(read, SAMRecord.class, conf);
    Assert.assertEquals(roundTrippedRead2, read, "\nActual read: " + roundTrippedRead2.getSAMString() + "\nExpected read: " + read.getSAMString());
}
 
Example 12
Source File: ReadContextFactoryTest.java    From hmftools with GNU General Public License v3.0 5 votes vote down vote up
@NotNull
static SAMRecord buildSamRecord(@NotNull final String cigar, @NotNull final String readString, @NotNull final String qualities) {
    final SAMRecord record = new SAMRecord(null);
    record.setAlignmentStart(1000);
    record.setCigarString(cigar);
    record.setReadString(readString);
    record.setReadNegativeStrandFlag(false);
    record.setBaseQualityString(qualities);
    record.setMappingQuality(20);
    record.setDuplicateReadFlag(false);
    record.setReadUnmappedFlag(false);
    return record;
}
 
Example 13
Source File: AbstractAlignmentMergerTest.java    From picard with MIT License 5 votes vote down vote up
@Test(dataProvider = "readPositionIgnoringSoftClips")
public void testGetReadPositionIgnoringSoftClips(final String cigarString, final int startPosition, final int queryPosition, final int expectedReadPosititon) {
    final SAMFileHeader newHeader = SAMRecordSetBuilder.makeDefaultHeader(SAMFileHeader.SortOrder.queryname, 100000,false);
    final SAMRecord rec = new SAMRecord(newHeader);

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

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

    Assert.assertEquals(readPosition, expectedReadPosititon);
}
 
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: TagReadWithGeneExonFunctionTest.java    From Drop-seq with MIT License 5 votes vote down vote up
private SAMRecord getFakeSplitRecord (final File testBAMFile, final int start1, final int end1, final int start2, final int end2, final boolean negativeStrand) {
	SamReader inputSam = SamReaderFactory.makeDefault().open(testBAMFile);
	SAMRecord r = inputSam.iterator().next();
	r.setAlignmentStart(start1);
	int length = (end1 -start1) +1;
	int gap = ((start2-1)-(end1+1)) +1;
	int length2 = (end2 -start2) +1;
	r.setCigarString(length+"M"+gap+"N"+length2+"M");
	r.setMappingQuality(255);
	r.setReadNegativeStrandFlag(negativeStrand);
	r.setReadBases(Arrays.copyOf(r.getReadBases(), length+length2));
	r.setBaseQualities(Arrays.copyOf(r.getBaseQualities(), length+length2));
	return (r);
}
 
Example 16
Source File: TagReadWithGeneFunctionTest.java    From Drop-seq with MIT License 5 votes vote down vote up
private SAMRecord getFakeRecord (final File testBAMFile, final int start, final int end, final boolean negativeStrand) {
	SamReader inputSam = SamReaderFactory.makeDefault().open(testBAMFile);
	SAMRecord r = inputSam.iterator().next();
	r.setAlignmentStart(start);
	int length = (end -start) +1;
	r.setCigarString(length+"M");
	r.setMappingQuality(255);
	r.setReadNegativeStrandFlag(negativeStrand);
	r.setReadBases(Arrays.copyOf(r.getReadBases(), length));
	r.setBaseQualities(Arrays.copyOf(r.getBaseQualities(), length));
	return (r);
}
 
Example 17
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 18
Source File: AlignmentsTagsTest.java    From cramtools with Apache License 2.0 4 votes vote down vote up
private void doTest(byte[] ref, int alignmentStart, byte[] readBases) throws IOException,
		CloneNotSupportedException {
	SAMSequenceRecord sequenceRecord = new SAMSequenceRecord("1", ref.length);
	SAMSequenceDictionary sequenceDictionary = new SAMSequenceDictionary();
	sequenceDictionary.addSequence(sequenceRecord);

	SAMFileHeader header = new SAMFileHeader();
	header.setSequenceDictionary(sequenceDictionary);
	SAMRecord samRecord = new SAMRecord(header);
	samRecord.setReadUnmappedFlag(false);
	samRecord.setAlignmentStart(alignmentStart);
	samRecord.setReferenceIndex(0);
	samRecord.setReadBases(readBases);
	samRecord.setCigarString(samRecord.getReadLength() + "M");

	ReferenceSource referenceSource = new ReferenceSource() {
		@Override
		public synchronized ReferenceRegion getRegion(SAMSequenceRecord record, int start_1based,
				int endInclusive_1based) throws IOException {
			int zbInclusiveStart = start_1based - 1;
			int zbExlcusiveEnd = endInclusive_1based;
			return new ReferenceRegion(Arrays.copyOfRange(ref, zbInclusiveStart, zbExlcusiveEnd),
					sequenceRecord.getSequenceIndex(), sequenceRecord.getSequenceName(), start_1based);
		}
	};

	AlignmentsTags.calculateMdAndNmTags(samRecord, referenceSource, sequenceDictionary, true, true);

	SAMRecord checkRecord = (SAMRecord) samRecord.clone();
	SequenceUtil.calculateMdAndNmTags(checkRecord, ref, true, true);
	// System.out.printf("TEST: ref %s, start %d, read bases %s\n", new
	// String(ref), alignmentStart, new String(
	// readBases));
	// System.out
	// .println(referenceSource.getRegion(sequenceRecord, alignmentStart,
	// alignmentStart + readBases.length));
	// System.out.printf("NM:  %s x %s\n", samRecord.getAttribute("NM"),
	// checkRecord.getAttribute("NM"));
	// System.out.printf("MD: %s x %s\n", samRecord.getAttribute("MD"),
	// checkRecord.getAttribute("MD"));

	Assert.assertEquals(checkRecord.getAttribute("NM"), samRecord.getAttribute("NM"));
	Assert.assertEquals(checkRecord.getAttribute("MD"), samRecord.getAttribute("MD"));
}
 
Example 19
Source File: ReorderSam.java    From picard with MIT License 4 votes vote down vote up
/**
 * Helper function that writes reads from iterator it into writer out, updating each SAMRecord along the way
 * according to the newOrder mapping from dictionary index -> index.  Name is used for printing only.
 */
private void writeReads(final SAMFileWriter out,
                        final SAMRecordIterator it,
                        final Map<Integer, Integer> newOrder,
                        final String name) {
    long counter = 0;
    log.info("  Processing " + name);

    while (it.hasNext()) {
        counter++;
        final SAMRecord read = it.next();
        final int oldRefIndex = read.getReferenceIndex();
        final int oldMateIndex = read.getMateReferenceIndex();
        final int newRefIndex = newOrderIndex(read, oldRefIndex, newOrder);

        read.setHeader(out.getFileHeader());
        read.setReferenceIndex(newRefIndex);

        // read becoming unmapped
        if (oldRefIndex != NO_ALIGNMENT_REFERENCE_INDEX &&
                newRefIndex == NO_ALIGNMENT_REFERENCE_INDEX) {
            read.setAlignmentStart(NO_ALIGNMENT_START);
            read.setReadUnmappedFlag(true);
            read.setCigarString(SAMRecord.NO_ALIGNMENT_CIGAR);
            read.setMappingQuality(SAMRecord.NO_MAPPING_QUALITY);
        }

        final int newMateIndex = newOrderIndex(read, oldMateIndex, newOrder);
        if (oldMateIndex != NO_ALIGNMENT_REFERENCE_INDEX &&
                newMateIndex == NO_ALIGNMENT_REFERENCE_INDEX) { // mate becoming unmapped
            read.setMateAlignmentStart(NO_ALIGNMENT_START);
            read.setMateUnmappedFlag(true);
            read.setAttribute(SAMTag.MC.name(), null);      // Set the Mate Cigar String to null
        }
        read.setMateReferenceIndex(newMateIndex);

        out.addAlignment(read);
    }

    it.close();
    log.info("Wrote " + counter + " reads");
}
 
Example 20
Source File: WriteBAMFn.java    From dataflow-java with Apache License 2.0 4 votes vote down vote up
@ProcessElement
public void processElement(DoFn<Read, String>.ProcessContext c, BoundedWindow window)
    throws Exception {

  this.window = window;

  if (headerInfo == null) {
    headerInfo = c.sideInput(headerView);
  }
  final Read read = c.element();

  if (readCount == 0) {

    shardContig = KeyReadsFn.shardKeyForRead(read, 1);
    sequenceIndex = headerInfo.header.getSequenceIndex(shardContig.referenceName);
    final boolean isFirstShard = headerInfo.shardHasFirstRead(shardContig);
    final String outputFileName = options.getOutput();
    shardName = outputFileName + "-" + String.format("%012d", sequenceIndex) + "-"
        + shardContig.referenceName
        + ":" + String.format("%012d", shardContig.start);
    LOG.info("Writing shard file " + shardName);
    final OutputStream outputStream =
        Channels.newOutputStream(
            new GcsUtil.GcsUtilFactory().create(options)
              .create(GcsPath.fromUri(shardName),
                  BAMIO.BAM_INDEX_FILE_MIME_TYPE));
    ts = new TruncatedOutputStream(
        outputStream, BlockCompressedStreamConstants.EMPTY_GZIP_BLOCK.length);
    bw = new BAMBlockWriter(ts, null /*file*/);
    bw.setSortOrder(headerInfo.header.getSortOrder(), true);
    bw.setHeader(headerInfo.header);
    if (isFirstShard) {
      LOG.info("First shard - writing header to " + shardName);
      bw.writeHeader(headerInfo.header);
    }
  }
  SAMRecord samRecord = ReadUtils.makeSAMRecord(read, headerInfo.header);
  if (prevRead != null && prevRead.getAlignmentStart() > samRecord.getAlignmentStart()) {
    LOG.info("Out of order read " + prevRead.getAlignmentStart() + " " +
        samRecord.getAlignmentStart() + " during writing of shard " + shardName +
        " after processing " + readCount + " reads, min seen alignment is " +
        minAlignment + " and max is " + maxAlignment + ", this read is " +
        (samRecord.getReadUnmappedFlag() ? "unmapped" : "mapped") + " and its mate is " +
        (samRecord.getMateUnmappedFlag() ? "unmapped" : "mapped"));
    Metrics.counter(WriteBAMFn.class, "Out of order reads").inc();
    readCount++;
    hadOutOfOrder = true;
    return;
  }
  minAlignment = Math.min(minAlignment, samRecord.getAlignmentStart());
  maxAlignment = Math.max(maxAlignment, samRecord.getAlignmentStart());
  prevRead = samRecord;
  if (samRecord.getReadUnmappedFlag()) {
    if (!samRecord.getMateUnmappedFlag()) {
      samRecord.setReferenceName(samRecord.getMateReferenceName());
      samRecord.setAlignmentStart(samRecord.getMateAlignmentStart());
    }
    unmappedReadCount++;
  }
  bw.addAlignment(samRecord);
  readCount++;
}