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

The following examples show how to use htsjdk.samtools.SAMRecord#setMappingQuality() . 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 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 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: ReadContextCounterTest.java    From hmftools with GNU General Public License v3.0 6 votes vote down vote up
@NotNull
public static SAMRecord buildSamRecord(final int alignmentStart, @NotNull final String cigar, @NotNull final String readString,
        @NotNull final String qualities) {
    final SAMRecord record = new SAMRecord(null);
    record.setAlignmentStart(alignmentStart);
    record.setCigarString(cigar);
    record.setReadString(readString);
    record.setReadNegativeStrandFlag(false);
    record.setBaseQualityString(qualities);
    record.setMappingQuality(20);
    record.setDuplicateReadFlag(false);
    record.setReadUnmappedFlag(false);
    record.setProperPairFlag(true);
    record.setReadPairedFlag(true);
    return record;
}
 
Example 5
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 6
Source File: ReadPairTest.java    From Drop-seq with MIT License 6 votes vote down vote up
private List<SAMRecord> getPairedRead (final String name, final int contig, final int start1, final int start2) {
	List<SAMRecord> result = new ArrayList<> ();

	SAMRecordSetBuilder builder = new SAMRecordSetBuilder();
	builder.addPair(name, contig, start1, start2);

	Collection<SAMRecord> recs = builder.getRecords();

	for (SAMRecord r: recs) {
		if (r.getFirstOfPairFlag()) result.add(0, r);
		if (r.getSecondOfPairFlag()) result.add(1, r);
		r.setMappingQuality(10);
	}

	return (result);

}
 
Example 7
Source File: RefSequenceTest.java    From hmftools with GNU General Public License v3.0 5 votes vote down vote up
@NotNull
public static SAMRecord buildSamRecord(final int alignmentStart, @NotNull final String cigar, @NotNull final String readString,
        @NotNull final String qualities) {
    final SAMRecord record = new SAMRecord(null);
    record.setAlignmentStart(alignmentStart);
    record.setCigarString(cigar);
    record.setReadString(readString);
    record.setReadNegativeStrandFlag(false);
    record.setBaseQualityString(qualities);
    record.setMappingQuality(20);
    record.setDuplicateReadFlag(false);
    record.setReadUnmappedFlag(false);
    return record;
}
 
Example 8
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 9
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 10
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 11
Source File: SamUtils.java    From rtg-tools with BSD 2-Clause "Simplified" License 5 votes vote down vote up
/**
 * Convert an aligned record to being unmapped
 * @param record the record to adjust
 */
public static void convertToUnmapped(SAMRecord record) {
  record.setReadUnmappedFlag(true);
  record.setProperPairFlag(false);
  record.setReadNegativeStrandFlag(false);
  record.setMappingQuality(0);
  record.setInferredInsertSize(0);
  record.setCigarString("*");
  record.setAttribute(ATTRIBUTE_NUM_MISMATCHES, null);
}
 
Example 12
Source File: CleanSam.java    From picard with MIT License 5 votes vote down vote up
/**
 * Do the work after command line has been parsed.
 * RuntimeException may be thrown by this method, and are reported appropriately.
 *
 * @return program exit status.
 */
@Override
protected int doWork() {
    IOUtil.assertFileIsReadable(INPUT);
    IOUtil.assertFileIsWritable(OUTPUT);
    final SamReaderFactory factory = SamReaderFactory.makeDefault().referenceSequence(REFERENCE_SEQUENCE);
    if (VALIDATION_STRINGENCY == ValidationStringency.STRICT) {
        factory.validationStringency(ValidationStringency.LENIENT);
    }
    final SamReader reader = factory.open(INPUT);
    final SAMFileWriter writer = new SAMFileWriterFactory().makeSAMOrBAMWriter(reader.getFileHeader(), true, OUTPUT);
    final CloseableIterator<SAMRecord> it = reader.iterator();
    final ProgressLogger progress = new ProgressLogger(Log.getInstance(CleanSam.class));

    // If the read (or its mate) maps off the end of the alignment, clip it
    while (it.hasNext()) {
        final SAMRecord rec = it.next();

        // If the read (or its mate) maps off the end of the alignment, clip it
        AbstractAlignmentMerger.createNewCigarsIfMapsOffEndOfReference(rec);

        // check the read's mapping quality
        if (rec.getReadUnmappedFlag() && 0 != rec.getMappingQuality()) {
            rec.setMappingQuality(0);
        }

        writer.addAlignment(rec);
        progress.record(rec);
    }

    writer.close();
    it.close();
    CloserUtil.close(reader);
    return 0;
}
 
Example 13
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 14
Source File: TagReadWithGeneExonFunctionTest.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 15
Source File: TagReadWithGeneFunctionTest.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: FilterBamByTagTest.java    From Drop-seq with MIT License 4 votes vote down vote up
@Test
public void filterReadTest() {
	SAMRecord readHasAttribute = new SAMRecord(null);
	String tag = "XT";
	readHasAttribute.setAttribute(tag, "1");

	Set<String> values = new HashSet<>();
	values.add("1");

	SAMRecord readNoAttribute = new SAMRecord(null);

	FilterBamByTag t = new FilterBamByTag();
	// read has attribute, accept any value, want to retain read.
	boolean flag1 = t.filterRead(readHasAttribute, tag, null, true, null, false);
	Assert.assertFalse(flag1);

	// read has attribute, accept any value, want to filter read.
	boolean flag2 = t.filterRead(readHasAttribute, tag, null, false, null, false);
	Assert.assertTrue(flag2);

	// read has attribute, accept certain value, want to retain read.
	boolean flag3 = t.filterRead(readHasAttribute, tag, values, true, null, false);
	Assert.assertFalse(flag3);

	// read has attribute, accept certain value, want to filter read.
	boolean flag4 = t.filterRead(readHasAttribute, tag, values, false, null, false);
	Assert.assertTrue(flag4);

	// read does not have attribute, accept any value, want to retain read.
	boolean flag5 = t.filterRead(readNoAttribute, tag, null, true, null, false);
	Assert.assertTrue(flag5);

	// read does not have attribute, accept any value, want to filter read.
	boolean flag6 = t.filterRead(readNoAttribute, tag, null, false, null, false);
	Assert.assertFalse(flag6);

	// read does not have attribute, accept certain value, want to retain read.
	boolean flag7 = t.filterRead(readNoAttribute, tag, values, true, null, false);
	Assert.assertTrue(flag7);

	// read does not have attribute, accept certain value, want to filter read.
	boolean flag8 = t.filterRead(readNoAttribute, tag, values, false, null, false);
	Assert.assertFalse(flag8);
	
	// test map quality filtering
	
	readHasAttribute.setMappingQuality(10);
	boolean flag9 = t.filterRead(readHasAttribute, tag, null, true, 10, false);
	Assert.assertFalse(flag9);
	boolean flag10 = t.filterRead(readHasAttribute, tag, null, true, 20, false);
	Assert.assertTrue(flag10);
	
	// want to test partial matching.
	SAMRecord readHasGene = new SAMRecord(null);
	String geneNameTag = "gn";
	readHasGene.setAttribute(geneNameTag, "A,B");

	Set<String> geneValues = new HashSet<> (Arrays.asList("A"));
	boolean flagExact =t.filterRead(readHasGene, geneNameTag, geneValues, true, 0, false);
	boolean flagPartial =t.filterRead(readHasGene, geneNameTag, geneValues, true, 0, true);
	// filter the exact match
	Assert.assertTrue(flagExact);
	// don't filter the partial match
	Assert.assertFalse(flagPartial);
	
}
 
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: 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 19
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 20
Source File: MergeBamAlignmentTest.java    From picard with MIT License 4 votes vote down vote up
/**
 * @return a 2-element array in which the first element is the unmapped SAM, and the second the mapped SAM.
 */
private File[] createSamFilesToBeMerged(final MultipleAlignmentSpec[] specs) {
    try {
        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 SAMRecord unmappedRecord = new SAMRecord(header);

        unmappedRecord.setReadName("theRead");
        unmappedRecord.setReadString("ACGTACGTACGTACGT");
        unmappedRecord.setBaseQualityString("5555555555555555");
        unmappedRecord.setReadUnmappedFlag(true);

        final SAMFileWriter unmappedWriter = factory.makeSAMWriter(header, false, unmappedSam);
        unmappedWriter.addAlignment(unmappedRecord);
        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);
        for (final MultipleAlignmentSpec spec : specs) {
            final SAMRecord alignedRecord = new SAMRecord(header);
            alignedRecord.setReadName(unmappedRecord.getReadName());
            alignedRecord.setReadBases(unmappedRecord.getReadBases());
            alignedRecord.setBaseQualities(unmappedRecord.getBaseQualities());
            alignedRecord.setReferenceName(sequence);
            alignedRecord.setAlignmentStart(1);
            alignedRecord.setReadNegativeStrandFlag(spec.reverseStrand);
            alignedRecord.setCigarString(spec.cigar);
            alignedRecord.setMappingQuality(spec.mapQ);
            if (spec.oneOfTheBest) {
                alignedRecord.setAttribute(ONE_OF_THE_BEST_TAG, 1);
            }
            alignedWriter.addAlignment(alignedRecord);
        }
        alignedWriter.close();

        return new File[]{unmappedSam, alignedSam};
    } catch (IOException e) {
        throw new PicardException(e.getMessage(), e);
    }
}