Java Code Examples for htsjdk.samtools.SAMRecord#setReadBases()
The following examples show how to use
htsjdk.samtools.SAMRecord#setReadBases() .
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 |
|
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: PalindromeArtifactClipReadTransformerUnitTest.java From gatk with BSD 3-Clause "New" or "Revised" License | 6 votes |
|
private static GATKRead makeRead(final SAMFileHeader header, final String contig, final int readStart, final int fragmentLength, final byte[] bases, final byte[] qual, final String cigar) {
final SAMRecord read = new SAMRecord(header);
read.setReferenceName(contig);
read.setAlignmentStart(readStart);
read.setReadPairedFlag(true);
read.setReadUnmappedFlag(false);
read.setMateUnmappedFlag(false);
read.setMateReferenceName("Mate");
read.setReadNegativeStrandFlag(false);
read.setMateNegativeStrandFlag(true);
read.setReadBases(bases);
read.setBaseQualities(qual);
final int mateStart = readStart + fragmentLength - bases.length;
read.setMateAlignmentStart(mateStart);
read.setInferredInsertSize(fragmentLength);
read.setProperPairFlag(true);
read.setCigarString(cigar);
return new SAMRecordToGATKReadAdapter(read);
}
Example 3
| Source File: Read.java From cramtools with Apache License 2.0 | 6 votes |
|
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: 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: TagReadWithGeneFunctionTest.java From Drop-seq with MIT License | 5 votes |
|
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 6
| 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 7
| Source File: FastWgsMetricsCollectorTest.java From picard with MIT License | 5 votes |
|
@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 8
| Source File: BaseErrorAggregationTest.java From picard with MIT License | 5 votes |
|
@Test
public void testBaseErrorAggregation2() {
final SAMSequenceRecord samSequenceRecord = new SAMSequenceRecord("chr1", 200);
final SAMFileHeader samFileHeader = new SAMFileHeader();
samFileHeader.addSequence(samSequenceRecord);
final SAMRecordSetBuilder builder = new SAMRecordSetBuilder();
final List<SAMRecord> samRecords = builder.addPair("Read1234", 0, 1, 1,
false, false, "16M", "16M", true, false, 20);
final SAMRecord samRecord1 = samRecords.get(0);
final SAMRecord samRecord2 = samRecords.get(1);
samRecord1.setReadBases("CgTGtGGAcAAAgAAA".getBytes());
samRecord2.setReadBases("CcTGGtGAcAAAgAAA".getBytes());
final byte[] refBases = "CATGGGGAAAAAAAAA".getBytes();
BaseErrorAggregation<?> baseErrorAggregation = new BaseErrorAggregation<>(OverlappingReadsErrorCalculator::new, ReadBaseStratification.readDirectionStratifier);
final int length = getLengthAndAddBases(samSequenceRecord, samRecord1, samRecord2, refBases, baseErrorAggregation);
final ErrorMetric[] metrics = baseErrorAggregation.getMetrics();
final OverlappingErrorMetric metric1 = (OverlappingErrorMetric) metrics[0];
metric1.calculateDerivedFields();
Assert.assertEquals(metric1.COVARIATE, ReadBaseStratification.ReadDirection.POSITIVE.toString());
Assert.assertEquals(metric1.TOTAL_BASES, length);
Assert.assertEquals(metric1.NUM_BASES_WITH_OVERLAPPING_READS, length);
Assert.assertEquals(metric1.NUM_DISAGREES_WITH_REFERENCE_ONLY, 2L);
Assert.assertEquals(metric1.NUM_DISAGREES_WITH_REF_AND_MATE, 1L);
Assert.assertEquals(metric1.NUM_THREE_WAYS_DISAGREEMENT, 1L);
final OverlappingErrorMetric metric2 = (OverlappingErrorMetric) metrics[1];
metric2.calculateDerivedFields();
Assert.assertEquals(metric2.COVARIATE, ReadBaseStratification.ReadDirection.NEGATIVE.toString());
Assert.assertEquals(metric2.TOTAL_BASES, length);
Assert.assertEquals(metric2.NUM_BASES_WITH_OVERLAPPING_READS, length);
Assert.assertEquals(metric2.NUM_DISAGREES_WITH_REFERENCE_ONLY, 2L);
Assert.assertEquals(metric2.NUM_DISAGREES_WITH_REF_AND_MATE, 1L);
Assert.assertEquals(metric2.NUM_THREE_WAYS_DISAGREEMENT, 1L);
}
Example 9
| Source File: BaseErrorAggregationTest.java From picard with MIT License | 5 votes |
|
@Test
public void testBaseErrorAggregation() {
final SAMSequenceRecord samSequenceRecord = new SAMSequenceRecord("chr1", 200);
final SAMFileHeader samFileHeader = new SAMFileHeader();
samFileHeader.addSequence(samSequenceRecord);
final SAMRecordSetBuilder builder = new SAMRecordSetBuilder();
final List<SAMRecord> samRecords = builder.addPair("Read1234", 0, 1, 1,
false, false, "16M", "16M", true, false, 20);
final SAMRecord samRecord1 = samRecords.get(0);
final SAMRecord samRecord2 = samRecords.get(1);
samRecord1.setReadBases("CgTGtGGAcAAAgAAA".getBytes());
samRecord2.setReadBases("CcTGGtGAcAAAgAAA".getBytes());
final byte[] refBases = "CATGGGGAAAAAAAAA".getBytes();
BaseErrorAggregation<?> baseErrorAggregation = new BaseErrorAggregation<>(OverlappingReadsErrorCalculator::new, ReadBaseStratification.readOrdinalityStratifier);
final int length = getLengthAndAddBases(samSequenceRecord, samRecord1, samRecord2, refBases, baseErrorAggregation);
final ErrorMetric[] metrics = baseErrorAggregation.getMetrics();
final OverlappingErrorMetric metric1 = (OverlappingErrorMetric) metrics[0];
metric1.calculateDerivedFields();
Assert.assertEquals(metric1.COVARIATE, ReadBaseStratification.ReadOrdinality.FIRST.name());
Assert.assertEquals(metric1.TOTAL_BASES, length);
Assert.assertEquals(metric1.NUM_BASES_WITH_OVERLAPPING_READS, length);
Assert.assertEquals(metric1.NUM_DISAGREES_WITH_REFERENCE_ONLY, 2L);
Assert.assertEquals(metric1.NUM_DISAGREES_WITH_REF_AND_MATE, 1L);
Assert.assertEquals(metric1.NUM_THREE_WAYS_DISAGREEMENT, 1L);
final OverlappingErrorMetric metric2 = (OverlappingErrorMetric) metrics[1];
metric2.calculateDerivedFields();
Assert.assertEquals(metric2.COVARIATE, ReadBaseStratification.ReadOrdinality.SECOND.name());
Assert.assertEquals(metric2.TOTAL_BASES, length);
Assert.assertEquals(metric2.NUM_BASES_WITH_OVERLAPPING_READS, length);
Assert.assertEquals(metric2.NUM_DISAGREES_WITH_REFERENCE_ONLY, 2L);
Assert.assertEquals(metric2.NUM_DISAGREES_WITH_REF_AND_MATE, 1L);
Assert.assertEquals(metric2.NUM_THREE_WAYS_DISAGREEMENT, 1L);
}
Example 10
| Source File: TagReadWithGeneFunctionTest.java From Drop-seq with MIT License | 5 votes |
|
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 11
| Source File: TagReadWithGeneExonFunctionTest.java From Drop-seq with MIT License | 5 votes |
|
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 12
| Source File: MergeBamAlignmentTest.java From picard with MIT License | 5 votes |
|
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 13
| Source File: TagReadWithGeneExonFunctionTest.java From Drop-seq with MIT License | 5 votes |
|
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 14
| 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"));
}
Example 15
| Source File: AbstractAlignmentMergerTest.java From picard with MIT License | 4 votes |
|
@Test(dataProvider = "overlapReadData")
public void testOverlappedReadClipping(final int readLength, final int start1, final int start2, final String cigar1, final String cigar2,
final boolean strand1, final boolean strand2,
final int r1ExpectedAlignmentStart, final int r2ExpectedAlignmentStart,
final String expectedR1Cigar, final String expectedR2Cigar, final boolean hardClipOverlappingReads,
final String read1Bases, final String read2Bases, final String expectedR1Bases, final String expectedR2Bases,
final String expectedR1ClippedBases, final String expectedR2ClippedBases, final String read1Qualities,
final String read2Qualities, final String expectedR1Qualities, final String expectedR2Qualities,
final String expectedR1ClippedQualities, final String expectedR2ClippedQualities) {
final SAMRecordSetBuilder set = new SAMRecordSetBuilder();
set.setReadLength(readLength);
final List<SAMRecord> recs = set.addPair("q1", 0, start1, start2, false, false, cigar1, cigar2, strand1, strand2, 30);
final SAMRecord r1 = recs.get(0);
final SAMRecord r2 = recs.get(1);
r1.setReadBases(StringUtil.stringToBytes(read1Bases));
r2.setReadBases(StringUtil.stringToBytes(read2Bases));
r1.setBaseQualities(SAMUtils.fastqToPhred(read1Qualities));
r2.setBaseQualities(SAMUtils.fastqToPhred(read2Qualities));
AbstractAlignmentMerger.clipForOverlappingReads(r1, r2, hardClipOverlappingReads);
Assert.assertEquals(r1.getAlignmentStart(), r1ExpectedAlignmentStart);
Assert.assertEquals(r1.getCigarString(), expectedR1Cigar);
Assert.assertEquals(r2.getAlignmentStart(), r2ExpectedAlignmentStart);
Assert.assertEquals(r2.getCigarString(), expectedR2Cigar);
Assert.assertEquals(r1.getReadString(), expectedR1Bases);
Assert.assertEquals(r2.getReadString(), expectedR2Bases);
Assert.assertEquals(SAMUtils.phredToFastq(r1.getBaseQualities()), expectedR1Qualities);
Assert.assertEquals(SAMUtils.phredToFastq(r2.getBaseQualities()), expectedR2Qualities);
Assert.assertEquals(r1.getAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASES_TAG), expectedR1ClippedBases);
Assert.assertEquals(r2.getAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASES_TAG), expectedR2ClippedBases);
Assert.assertEquals(r1.getAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASE_QUALITIES_TAG), expectedR1ClippedQualities);
Assert.assertEquals(r2.getAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASE_QUALITIES_TAG), expectedR2ClippedQualities);
// Swap first and second read to ensure logic is correct for both F1R2 and F2R1
final SAMRecordSetBuilder setSwapped = new SAMRecordSetBuilder();
setSwapped.setReadLength(readLength);
final List<SAMRecord> recsSwapped = set.addPair("q1", 0, start2, start1, false, false, cigar2, cigar1, strand2, strand1, 30);
final SAMRecord r1Swapped = recsSwapped.get(0);
final SAMRecord r2Swapped = recsSwapped.get(1);
r1Swapped.setReadBases(StringUtil.stringToBytes(read2Bases));
r2Swapped.setReadBases(StringUtil.stringToBytes(read1Bases));
r1Swapped.setBaseQualities(SAMUtils.fastqToPhred(read2Qualities));
r2Swapped.setBaseQualities(SAMUtils.fastqToPhred(read1Qualities));
AbstractAlignmentMerger.clipForOverlappingReads(r1Swapped, r2Swapped, hardClipOverlappingReads);
Assert.assertEquals(r1Swapped.getAlignmentStart(), r2ExpectedAlignmentStart);
Assert.assertEquals(r1Swapped.getCigarString(), expectedR2Cigar);
Assert.assertEquals(r2Swapped.getAlignmentStart(), r1ExpectedAlignmentStart);
Assert.assertEquals(r2Swapped.getCigarString(), expectedR1Cigar);
Assert.assertEquals(r1Swapped.getReadString(), expectedR2Bases);
Assert.assertEquals(r2Swapped.getReadString(), expectedR1Bases);
Assert.assertEquals(SAMUtils.phredToFastq(r1Swapped.getBaseQualities()), expectedR2Qualities);
Assert.assertEquals(SAMUtils.phredToFastq(r2Swapped.getBaseQualities()), expectedR1Qualities);
Assert.assertEquals(r1Swapped.getAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASES_TAG), expectedR2ClippedBases);
Assert.assertEquals(r2Swapped.getAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASES_TAG), expectedR1ClippedBases);
Assert.assertEquals(r1Swapped.getAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASE_QUALITIES_TAG), expectedR2ClippedQualities);
Assert.assertEquals(r2Swapped.getAttribute(AbstractAlignmentMerger.HARD_CLIPPED_BASE_QUALITIES_TAG), expectedR1ClippedQualities);
}
Example 16
| Source File: MergeBamAlignmentTest.java From picard with MIT License | 4 votes |
|
/**
* 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 17
| Source File: BaseErrorAggregationTest.java From picard with MIT License | 4 votes |
|
@Test
public void testBaseErrorAggregation3() {
final SAMSequenceRecord samSequenceRecord = new SAMSequenceRecord("chr1", 200);
final SAMFileHeader samFileHeader = new SAMFileHeader();
samFileHeader.addSequence(samSequenceRecord);
final SAMRecordSetBuilder builder = new SAMRecordSetBuilder();
final List<SAMRecord> samRecords = builder.addPair("Read1234", 0, 1, 1,
false, false, "16M", "16M", true, false, 20);
final SAMRecord samRecord1 = samRecords.get(0);
final SAMRecord samRecord2 = samRecords.get(1);
samRecord1.setReadBases("CgTGtGGAcAAAgAAA".getBytes());
samRecord2.setReadBases("CcTGGtGAcAAAgAAA".getBytes());
final byte[] refBases = "CATGGGGAAAAAAAAA".getBytes();
BaseErrorAggregation<?> baseErrorAggregation = new BaseErrorAggregation<>(SimpleErrorCalculator::new, ReadBaseStratification.referenceBaseStratifier);
final int length = getLengthAndAddBases(samSequenceRecord, samRecord1, samRecord2, refBases, baseErrorAggregation);
final ErrorMetric[] metrics = baseErrorAggregation.getMetrics();
final BaseErrorMetric metricA = Arrays.stream(metrics).map(a -> (BaseErrorMetric) a)
.filter(a -> a.COVARIATE.equals("A")).findFirst().get();
metricA.calculateDerivedFields();
Assert.assertEquals(metricA.COVARIATE, "A");
Assert.assertEquals(metricA.TOTAL_BASES, 11);
Assert.assertEquals(metricA.ERROR_BASES, 3);
final BaseErrorMetric metricC = Arrays.stream(metrics).map(a -> (BaseErrorMetric) a)
.filter(a -> a.COVARIATE.equals("C")).findFirst().get();
metricA.calculateDerivedFields();
Assert.assertEquals(metricC.COVARIATE, "C");
Assert.assertEquals(metricC.TOTAL_BASES, 5);
Assert.assertEquals(metricC.ERROR_BASES, 1);
final BaseErrorMetric metricG = Arrays.stream(metrics).map(a -> (BaseErrorMetric) a)
.filter(a -> a.COVARIATE.equals("G")).findFirst().get();
metricA.calculateDerivedFields();
Assert.assertEquals(metricG.COVARIATE, "G");
Assert.assertEquals(metricG.TOTAL_BASES, 5);
Assert.assertEquals(metricG.ERROR_BASES, 1);
final BaseErrorMetric metricT = Arrays.stream(metrics).map(a -> (BaseErrorMetric) a)
.filter(a -> a.COVARIATE.equals("T")).findFirst().get();
metricA.calculateDerivedFields();
Assert.assertEquals(metricT.COVARIATE, "T");
Assert.assertEquals(metricT.TOTAL_BASES, 11);
Assert.assertEquals(metricT.ERROR_BASES, 3);
}
Example 18
| Source File: Utils.java From cramtools with Apache License 2.0 | 4 votes |
|
public static void changeReadLength(SAMRecord record, int newLength) {
if (newLength == record.getReadLength())
return;
if (newLength < 1 || newLength >= record.getReadLength())
throw new IllegalArgumentException("Cannot change read length to " + newLength);
List<CigarElement> newCigarElements = new ArrayList<CigarElement>();
int len = 0;
for (CigarElement ce : record.getCigar().getCigarElements()) {
switch (ce.getOperator()) {
case D:
break;
case S:
// dump = true;
// len -= ce.getLength();
// break;
case M:
case I:
case X:
len += ce.getLength();
break;
default:
throw new IllegalArgumentException("Unexpected cigar operator: " + ce.getOperator() + " in cigar "
+ record.getCigarString());
}
if (len <= newLength) {
newCigarElements.add(ce);
continue;
}
CigarElement newCe = new CigarElement(ce.getLength() - (record.getReadLength() - newLength),
ce.getOperator());
if (newCe.getLength() > 0)
newCigarElements.add(newCe);
break;
}
byte[] newBases = new byte[newLength];
System.arraycopy(record.getReadBases(), 0, newBases, 0, newLength);
record.setReadBases(newBases);
byte[] newScores = new byte[newLength];
System.arraycopy(record.getBaseQualities(), 0, newScores, 0, newLength);
record.setCigar(new Cigar(newCigarElements));
}
Example 19
| 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 20
| Source File: ClusterDataToSamConverter.java From picard with MIT License | 4 votes |
|
/**
* Creates a new SAM record from the basecall data
*/
private SAMRecord createSamRecord(final ReadData readData, final String readName, final boolean isPf, final boolean firstOfPair,
final String unmatchedBarcode, final String barcodeQuality,
final List<String> molecularIndexes, final List<String> molecularIndexQualities) {
final SAMRecord sam = new SAMRecord(null);
sam.setReadName(readName);
sam.setReadBases(readData.getBases());
sam.setBaseQualities(readData.getQualities());
// Flag values
sam.setReadPairedFlag(isPairedEnd);
sam.setReadUnmappedFlag(true);
sam.setReadFailsVendorQualityCheckFlag(!isPf);
if (isPairedEnd) {
sam.setMateUnmappedFlag(true);
sam.setFirstOfPairFlag(firstOfPair);
sam.setSecondOfPairFlag(!firstOfPair);
}
if (filters.filterOut(sam)) {
sam.setAttribute(ReservedTagConstants.XN, 1);
}
if (this.readGroupId != null) {
sam.setAttribute(SAMTag.RG.name(), readGroupId);
}
// If it's a barcoded run and it has been decided that the original BC value should be added to the record, do it
if (unmatchedBarcode != null) {
sam.setAttribute(SAMTag.BC.name(), unmatchedBarcode);
if (barcodeQuality != null ) {
sam.setAttribute(SAMTag.QT.name(), barcodeQuality);
}
}
if (!molecularIndexes.isEmpty()) {
if (!this.MOLECULAR_INDEX_TAG.isEmpty()) {
sam.setAttribute(this.MOLECULAR_INDEX_TAG, String.join(MOLECULAR_INDEX_, molecularIndexes));
}
if (!this.MOLECULAR_INDEX_QUALITY_TAG.isEmpty()) {
sam.setAttribute(this.MOLECULAR_INDEX_QUALITY_TAG, String.join(MOLECULAR_INDEX_, molecularIndexQualities));
}
if (!this.tagPerMolecularIndex.isEmpty()) {
if (tagPerMolecularIndex.size() != molecularIndexes.size()) {
throw new PicardException("Found " + molecularIndexes.size() + " molecular indexes but only " + tagPerMolecularIndex.size() + " SAM tags given.");
}
for (int i = 0; i < this.tagPerMolecularIndex.size(); i++) {
sam.setAttribute(this.tagPerMolecularIndex.get(i), molecularIndexes.get(i));
}
}
}
return sam;
}
