Java Code Examples for htsjdk.samtools.SAMRecord#setReferenceIndex()
The following examples show how to use
htsjdk.samtools.SAMRecord#setReferenceIndex() .
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: SamCompareUtilsTest.java From rtg-tools with BSD 2-Clause "Simplified" License | 5 votes |
|
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 2
| Source File: UmiUtilTest.java From picard with MIT License | 5 votes |
|
@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 3
| Source File: HaplotypeBAMWriter.java From gatk with BSD 3-Clause "New" or "Revised" License | 5 votes |
|
/**
* 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 4
| Source File: BAMRecordViewTest.java From cramtools with Apache License 2.0 | 5 votes |
|
@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 5
| Source File: RevertSam.java From picard with MIT License | 4 votes |
|
/**
* Takes an individual SAMRecord and applies the set of changes/reversions to it that
* have been requested by program level options.
*/
public void revertSamRecord(final SAMRecord rec) {
if (RESTORE_ORIGINAL_QUALITIES) {
final byte[] oq = rec.getOriginalBaseQualities();
if (oq != null) {
rec.setBaseQualities(oq);
rec.setOriginalBaseQualities(null);
}
}
if (REMOVE_DUPLICATE_INFORMATION) {
rec.setDuplicateReadFlag(false);
}
if (REMOVE_ALIGNMENT_INFORMATION) {
if (rec.getReadNegativeStrandFlag()) {
rec.reverseComplement(true);
rec.setReadNegativeStrandFlag(false);
}
// Remove all alignment based information about the read itself
rec.setReferenceIndex(SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX);
rec.setAlignmentStart(SAMRecord.NO_ALIGNMENT_START);
rec.setCigarString(SAMRecord.NO_ALIGNMENT_CIGAR);
rec.setMappingQuality(SAMRecord.NO_MAPPING_QUALITY);
rec.setInferredInsertSize(0);
rec.setNotPrimaryAlignmentFlag(false);
rec.setProperPairFlag(false);
rec.setReadUnmappedFlag(true);
// Then remove any mate flags and info related to alignment
rec.setMateAlignmentStart(SAMRecord.NO_ALIGNMENT_START);
rec.setMateNegativeStrandFlag(false);
rec.setMateReferenceIndex(SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX);
rec.setMateUnmappedFlag(rec.getReadPairedFlag());
if (RESTORE_HARDCLIPS) {
String hardClippedBases = rec.getStringAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASES_TAG);
String hardClippedQualities = rec.getStringAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASE_QUALITIES_TAG);
if (hardClippedBases != null && hardClippedQualities != null) {
// Record has already been reverse complemented if this was on the negative strand
rec.setReadString(rec.getReadString() + hardClippedBases);
rec.setBaseQualities(SAMUtils.fastqToPhred(SAMUtils.phredToFastq(rec.getBaseQualities()) + hardClippedQualities));
// Remove hard clipping storage tags
rec.setAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASES_TAG, null);
rec.setAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASE_QUALITIES_TAG, null);
}
}
// And then remove any tags that are calculated from the alignment
ATTRIBUTE_TO_CLEAR.forEach(tag -> rec.setAttribute(tag, null));
}
}
Example 6
| Source File: ReorderSam.java From picard with MIT License | 4 votes |
|
/**
* 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 7
| Source File: Utils.java From cramtools with Apache License 2.0 | 4 votes |
|
/**
* Copied from net.sf.picard.sam.SamPairUtil. This is a more permissive
* version of the method, which does not reset alignment start and reference
* for unmapped reads.
*
* @param rec1
* @param rec2
* @param header
*/
public static void setLooseMateInfo(final SAMRecord rec1, final SAMRecord rec2, final SAMFileHeader header) {
if (rec1.getReferenceName() != SAMRecord.NO_ALIGNMENT_REFERENCE_NAME
&& rec1.getReferenceIndex() == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX)
rec1.setReferenceIndex(header.getSequenceIndex(rec1.getReferenceName()));
if (rec2.getReferenceName() != SAMRecord.NO_ALIGNMENT_REFERENCE_NAME
&& rec2.getReferenceIndex() == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX)
rec2.setReferenceIndex(header.getSequenceIndex(rec2.getReferenceName()));
// If neither read is unmapped just set their mate info
if (!rec1.getReadUnmappedFlag() && !rec2.getReadUnmappedFlag()) {
rec1.setMateReferenceIndex(rec2.getReferenceIndex());
rec1.setMateAlignmentStart(rec2.getAlignmentStart());
rec1.setMateNegativeStrandFlag(rec2.getReadNegativeStrandFlag());
rec1.setMateUnmappedFlag(false);
rec1.setAttribute(SAMTag.MQ.name(), rec2.getMappingQuality());
rec2.setMateReferenceIndex(rec1.getReferenceIndex());
rec2.setMateAlignmentStart(rec1.getAlignmentStart());
rec2.setMateNegativeStrandFlag(rec1.getReadNegativeStrandFlag());
rec2.setMateUnmappedFlag(false);
rec2.setAttribute(SAMTag.MQ.name(), rec1.getMappingQuality());
}
// Else if they're both unmapped set that straight
else if (rec1.getReadUnmappedFlag() && rec2.getReadUnmappedFlag()) {
rec1.setMateNegativeStrandFlag(rec2.getReadNegativeStrandFlag());
rec1.setMateUnmappedFlag(true);
rec1.setAttribute(SAMTag.MQ.name(), null);
rec1.setInferredInsertSize(0);
rec2.setMateNegativeStrandFlag(rec1.getReadNegativeStrandFlag());
rec2.setMateUnmappedFlag(true);
rec2.setAttribute(SAMTag.MQ.name(), null);
rec2.setInferredInsertSize(0);
}
// And if only one is mapped copy it's coordinate information to the
// mate
else {
final SAMRecord mapped = rec1.getReadUnmappedFlag() ? rec2 : rec1;
final SAMRecord unmapped = rec1.getReadUnmappedFlag() ? rec1 : rec2;
mapped.setMateReferenceIndex(unmapped.getReferenceIndex());
mapped.setMateAlignmentStart(unmapped.getAlignmentStart());
mapped.setMateNegativeStrandFlag(unmapped.getReadNegativeStrandFlag());
mapped.setMateUnmappedFlag(true);
mapped.setInferredInsertSize(0);
unmapped.setMateReferenceIndex(mapped.getReferenceIndex());
unmapped.setMateAlignmentStart(mapped.getAlignmentStart());
unmapped.setMateNegativeStrandFlag(mapped.getReadNegativeStrandFlag());
unmapped.setMateUnmappedFlag(false);
unmapped.setInferredInsertSize(0);
}
boolean firstIsFirst = rec1.getAlignmentStart() < rec2.getAlignmentStart();
int insertSize = firstIsFirst ? SamPairUtil.computeInsertSize(rec1, rec2) : SamPairUtil.computeInsertSize(rec2,
rec1);
rec1.setInferredInsertSize(firstIsFirst ? insertSize : -insertSize);
rec2.setInferredInsertSize(firstIsFirst ? -insertSize : insertSize);
}
Example 8
| Source File: TestBAMRecordView.java From cramtools with Apache License 2.0 | 4 votes |
|
@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 9
| Source File: AlignmentsTagsTest.java From cramtools with Apache License 2.0 | 4 votes |
|
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"));
}
