public void ensure(final int start, final int end) { if (end - start > bases.length) throw new RuntimeException("Window is too small for start " + start + " end " + end); if (position < start) moveForwardTo(start); }
private ReferenceTracks(final int sequenceId, final String sequenceName, final byte[] reference, final int windowSize) { this.sequenceId = sequenceId; this.sequenceName = sequenceName; this.reference = reference; bases = new byte[Math.min(windowSize, reference.length)]; coverage = new short[Math.min(windowSize, reference.length)]; mismatches = new short[Math.min(windowSize, reference.length)]; position = 1; reset(); }
referenceTracks.ensureRange(alStart, alSpan); for (int i = 0; i < cigarElement.getLength(); i++) { final boolean match = samRecord.getReadBases()[pos + i] == referenceTracks .baseAt(refPos + i); mask[pos + i] = (baseCategory.type == BaseCategoryType.MATCH && match) || (baseCategory.type == BaseCategoryType.MISMATCH && !match); case X: for (int i = 0; i < cigarElement.getLength(); i++) mask[pos + i - 1] = referenceTracks.coverageAt(refPos + i) < baseCategory.param; break; default: case PILEUP: for (int i = 0; i < qualityScores.length; i++) if (referenceTracks.mismatchesAt(alStart + i) > baseCategory.param) mask[i] = true; break;
private static void updateTracks(final List<SAMRecord> samRecords, final ReferenceTracks tracks) { for (final SAMRecord samRecord : samRecords) { if (samRecord.getAlignmentStart() != SAMRecord.NO_ALIGNMENT_START && samRecord.getReferenceIndex() == tracks.getSequenceId()) { int refPos = samRecord.getAlignmentStart(); int readPos = 0; if (cigarElement.getOperator().consumesReferenceBases()) { for (int elementIndex = 0; elementIndex < cigarElement.getLength(); elementIndex++) tracks.addCoverage(refPos + elementIndex, 1); for (int pos = readPos; pos < cigarElement.getLength(); pos++) { final byte readBase = samRecord.getReadBases()[readPos + pos]; final byte refBase = tracks.baseAt(refPos + pos); if (readBase != refBase) tracks.addMismatches(refPos + pos, 1);
private static void updateTracks(final List<SAMRecord> samRecords, final ReferenceTracks tracks) { for (final SAMRecord samRecord : samRecords) { if (samRecord.getAlignmentStart() != SAMRecord.NO_ALIGNMENT_START) { int refPos = samRecord.getAlignmentStart(); int readPos = 0; for (final CigarElement cigarElement : samRecord.getCigar().getCigarElements()) { if (cigarElement.getOperator().consumesReferenceBases()) { for (int elementIndex = 0; elementIndex < cigarElement.getLength(); elementIndex++) tracks.addCoverage(refPos + elementIndex, 1); } switch (cigarElement.getOperator()) { case M: case X: case EQ: for (int pos = readPos; pos < cigarElement.getLength(); pos++) { final byte readBase = samRecord.getReadBases()[readPos + pos]; final byte refBase = tracks.baseAt(refPos + pos); if (readBase != refBase) tracks.addMismatches(refPos + pos, 1); } break; default: break; } readPos += cigarElement.getOperator().consumesReadBases() ? cigarElement.getLength() : 0; refPos += cigarElement.getOperator().consumesReferenceBases() ? cigarElement.getLength() : 0; } } } }
@Ignore("Broken test.") @Test public void test4() { String line2 = "98738 1187 20 18 29 99M1S = 1000253 432 AGCGGGGATATATAAAGGCTCAAAATTACTTTTTATATGGACAACTCTCTACTGCTTTGAGATGACTGATACTCATATTGATGGAGCTTTATCAAGAAAT !\"#$%&'()*+-./0'''''''''''#'#'#'''''''#''''#'''''''''##''''#'#''#'''''#'''''''''##''''#''##''''''''?"; String seqName = "20"; List<String> lines = Arrays.asList(new String[] { line2 }); byte[] ref = "CTGGTCTTAGTTCCGCAAGTGGGTATATATAAAGGCTCAAAATCAATCTTTATATTGACATCTCTCTACTTATTTGTGTTGTCTGATGCTCATATTGTAGGAGATTCCTCAAGAAAGG" .getBytes(); ReferenceTracks tracks = new ReferenceTracks(0, seqName, ref); QualityScorePreservation p = new QualityScorePreservation("R40X10-N40-U40"); for (int i = 0; i < ref.length; i++) tracks.addCoverage(i + 1, 66); for (String line : lines) { SAMRecord record = buildSAMRecord(seqName, line); Sam2CramRecordFactory f = new Sam2CramRecordFactory(ref, record.getHeader(), CramVersions.CRAM_v3); CramCompressionRecord cramRecord = f.createCramRecord(record); p.addQualityScores(record, cramRecord, tracks); if (!cramRecord.isForcePreserveQualityScores()) { CramNormalizer.restoreQualityScores((byte) 30, Collections.singletonList(cramRecord)); } StringBuffer sb = new StringBuffer(); sb.append(record.getBaseQualityString()); sb.append("\n"); sb.append(SAMUtils.phredToFastq(cramRecord.qualityScores)); assertArrayEquals(sb.toString(), record.getBaseQualities(), cramRecord.qualityScores); } }
tracks = new ReferenceTracks(refSeqIndex, refSeqName, refs); tracks.ensureRange(start, stop - start + 1); updateTracks(samRecords, tracks);
if (bases == null) throw new RuntimeException("Bases not found for id " + refId); tracks = new ReferenceTracks(refId, sequence.getSequenceName(), context.referenceSource.getReferenceBases(sequence, true)); context.tracks.put(refId, tracks);
tracks.ensureRange(span.getStart(), span.getSpan()); updateTracks(samRecords, tracks);
referenceTracks.ensureRange(alStart, alSpan); for (int i = 0; i < cigarElement.getLength(); i++) { final boolean match = samRecord.getReadBases()[pos + i] == referenceTracks .baseAt(refPos + i); mask[pos + i] = (baseCategory.type == BaseCategoryType.MATCH && match) || (baseCategory.type == BaseCategoryType.MISMATCH && !match); case X: for (int i = 0; i < cigarElement.getLength(); i++) mask[pos + i - 1] = referenceTracks.coverageAt(refPos + i) < baseCategory.param; break; default: case PILEUP: for (int i = 0; i < qualityScores.length; i++) if (referenceTracks.mismatchesAt(alStart + i) > baseCategory.param) mask[i] = true; break;
private static void updateTracks(final List<SAMRecord> samRecords, final ReferenceTracks tracks) { for (final SAMRecord samRecord : samRecords) { if (samRecord.getAlignmentStart() != SAMRecord.NO_ALIGNMENT_START) { int refPos = samRecord.getAlignmentStart(); int readPos = 0; for (final CigarElement cigarElement : samRecord.getCigar().getCigarElements()) { if (cigarElement.getOperator().consumesReferenceBases()) { for (int elementIndex = 0; elementIndex < cigarElement.getLength(); elementIndex++) tracks.addCoverage(refPos + elementIndex, 1); } switch (cigarElement.getOperator()) { case M: case X: case EQ: for (int pos = readPos; pos < cigarElement.getLength(); pos++) { final byte readBase = samRecord.getReadBases()[readPos + pos]; final byte refBase = tracks.baseAt(refPos + pos); if (readBase != refBase) tracks.addMismatches(refPos + pos, 1); } break; default: break; } readPos += cigarElement.getOperator().consumesReadBases() ? cigarElement.getLength() : 0; refPos += cigarElement.getOperator().consumesReferenceBases() ? cigarElement.getLength() : 0; } } } }
tracks = new ReferenceTracks(refSeqIndex, refSeqName, refs); tracks.ensureRange(start, stop - start + 1); updateTracks(samRecords, tracks);
private byte[] compressScores (SAMRecord record, byte[] ref, QualityScorePreservation p) { ReferenceTracks tracks = new ReferenceTracks(0, record.getReferenceName(), ref); Sam2CramRecordFactory f = new Sam2CramRecordFactory(ref, record.getHeader(), CramVersions.CRAM_v3); CramCompressionRecord cramRecord = f.createCramRecord(record); p.addQualityScores(record, cramRecord, tracks); if (!cramRecord.isForcePreserveQualityScores()) { CramNormalizer.restoreQualityScores((byte) 30, Collections.singletonList(cramRecord)); } return cramRecord.qualityScores; } }
t.ensureRange(alStart, alSpan); if (ce.getOperator().consumesReadBases()) { for (int i = 0; i < ce.getLength(); i++) { final boolean match = s.getReadBases()[pos + i] == t.baseAt(refPos + i); mask[pos + i] = (c.type == BaseCategoryType.MATCH && match) || (c.type == BaseCategoryType.MISMATCH && !match); case X: for (int i = 0; i < ce.getLength(); i++) mask[pos + i - 1] = t.coverageAt(refPos + i) < c.param; break; default: case PILEUP: for (int i = 0; i < qs.length; i++) if (t.mismatchesAt(alStart + i) > c.param) mask[i] = true; break;
@Ignore("Broken test.") @Test public void test3() { String line1 = "98573 1107 20 1 60 100M = 999587 -415 CTGGTCTTAGTTCCGCAAGTGGGTATATATAAAGGCTCAAAATCAATCTTTATATTGACATCTCTCTACTTATTTGTGTTGTCTGATGCTCATATTGTAG ::A<<=D@BBC;C9=7DEEBHDEHHACEEBEEEDEE=EFFHEEFFFEHEF@HFBCEFEHFEHEHFEHDHHHFHHHEHHHHDFHHHHHGHHHHHHHHHHHH"; String line2 = "98738 1187 20 18 29 99M1S = 1000253 432 AGCGGGGATATATAAAGGCTCAAAATTACTTTTTATATGGACAACTCTCTACTGCTTTGAGATGACTGATACTCATATTGATGGAGCTTTATCAAGAAAT !\"#$%&'()*+-./0'''''''''''#'#'#'''''''#''''#'''''''''##''''#'#''#'''''#'''''''''##''''#''##''''''''?"; String seqName = "20"; List<String> lines = Arrays.asList(new String[] { line2, line1 }); byte[] ref = "CTGGTCTTAGTTCCGCAAGTGGGTATATATAAAGGCTCAAAATCAATCTTTATATTGACATCTCTCTACTTATTTGTGTTGTCTGATGCTCATATTGTAGGAGATTCCTCAAGAAAGG" .getBytes(); ReferenceTracks tracks = new ReferenceTracks(0, seqName, ref); QualityScorePreservation p = new QualityScorePreservation("R8-N40-M40-D40"); for (String line : lines) { SAMRecord record = buildSAMRecord(seqName, line); Sam2CramRecordFactory f = new Sam2CramRecordFactory(ref, record.getHeader(), CramVersions.CRAM_v3); CramCompressionRecord cramRecord = f.createCramRecord(record); p.addQualityScores(record, cramRecord, tracks); if (!cramRecord.isForcePreserveQualityScores()) { CramNormalizer.restoreQualityScores((byte) 30, Collections.singletonList(cramRecord)); } StringBuffer sb = new StringBuffer(); sb.append(record.getBaseQualityString()); sb.append("\n"); sb.append(SAMUtils.phredToFastq(cramRecord.qualityScores)); assertArrayEquals(sb.toString(), record.getBaseQualities(), cramRecord.qualityScores); } }
public void ensureRange(final int start, final int length) { if (start < position) throw new RuntimeException("Cannot move the window backwards: " + start); if (start > position || start + length > position + bases.length) moveForwardTo(start); }
private ReferenceTracks(final int sequenceId, final String sequenceName, final byte[] reference, final int windowSize) { this.sequenceId = sequenceId; this.sequenceName = sequenceName; this.reference = reference; bases = new byte[Math.min(windowSize, reference.length)]; coverage = new short[Math.min(windowSize, reference.length)]; mismatches = new short[Math.min(windowSize, reference.length)]; position = 1; reset(); }
public void ensure(final int start, final int end) { if (end - start > bases.length) throw new RuntimeException("Window is too small for start " + start + " end " + end); if (position < start) moveForwardTo(start); }
private ReferenceTracks(final int sequenceId, final String sequenceName, final byte[] reference, final int windowSize) { this.sequenceId = sequenceId; this.sequenceName = sequenceName; this.reference = reference; bases = new byte[Math.min(windowSize, reference.length)]; coverage = new short[Math.min(windowSize, reference.length)]; mismatches = new short[Math.min(windowSize, reference.length)]; position = 1; reset(); }
public void ensureRange(final int start, final int length) { if (start < position) throw new RuntimeException("Cannot move the window backwards: " + start); if (start > position || start + length > position + bases.length) moveForwardTo(start); }