Skip to content

Commit

Permalink
Small test harness for readfeatures
Browse files Browse the repository at this point in the history
  • Loading branch information
cmdcolin committed Nov 21, 2022
1 parent 3e79385 commit 3431a85
Show file tree
Hide file tree
Showing 4 changed files with 133 additions and 100 deletions.
103 changes: 3 additions & 100 deletions plugins/alignments/src/CramAdapter/CramSlightlyLazyFeature.ts
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ import {
} from '@jbrowse/core/util/simpleFeature'
import { CramRecord } from '@gmod/cram'
import CramAdapter from './CramAdapter'
import { readFeaturesToMismatches } from './util'

export interface Mismatch {
qual?: number
Expand Down Expand Up @@ -295,107 +296,9 @@ export default class CramSlightlyLazyFeature implements Feature {
}

_get_mismatches(): Mismatch[] {
const readFeatures = this.get('cram_read_features')
const readFeatures = this.record.readFeatures
const qual = this.qualRaw()
if (!readFeatures) {
return []
}
const start = this.get('start')
const mismatches: Mismatch[] = new Array(readFeatures.length)
let j = 0
let insLen = 0

let refPos = 0
let sublen = 0
let lastPos = this.record.alignmentStart

for (let i = 0; i < readFeatures.length; i++) {
const f = readFeatures[i]
const { code, pos, data, sub, ref } = f
sublen = refPos - lastPos
lastPos = refPos

if (sublen && insLen > 0) {
mismatches[j++] = {
start: refPos,
type: 'insertion',
base: `${insLen}`,
length: 0,
}
insLen = 0
}
refPos = f.refPos - 1 - start

if (code === 'X') {
// substitution
mismatches[j++] = {
start: refPos,
length: 1,
base: sub,
qual: qual?.[pos],
altbase: ref,
type: 'mismatch',
}
} else if (code === 'I') {
// insertion
mismatches[j++] = {
start: refPos,
type: 'insertion',
base: `${data.length}`,
length: 0,
}
} else if (code === 'N') {
// reference skip
mismatches[j++] = {
type: 'skip',
length: data,
start: refPos,
base: 'N',
}
} else if (code === 'S') {
// soft clip
const len = data.length
mismatches[j++] = {
start: refPos,
type: 'softclip',
base: `S${len}`,
cliplen: len,
length: 1,
}
} else if (code === 'P') {
// padding
} else if (code === 'H') {
// hard clip
const len = data
mismatches[j++] = {
start: refPos,
type: 'hardclip',
base: `H${len}`,
cliplen: len,
length: 1,
}
} else if (code === 'D') {
// deletion
mismatches[j++] = {
type: 'deletion',
length: data,
start: refPos,
base: '*',
}
} else if (code === 'b') {
// stretch of bases
} else if (code === 'q') {
// stretch of qual scores
} else if (code === 'B') {
// a pair of [base, qual]
} else if (code === 'i') {
// single-base insertion
// insertion
insLen++
} else if (code === 'Q') {
// single quality value
}
}
return mismatches.slice(0, j)
return readFeaturesToMismatches(readFeatures, start, qual)
}
}
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@
// Jest Snapshot v1, https://goo.gl/fbAQLP

exports[`cram read features 1`] = `[]`;
7 changes: 7 additions & 0 deletions plugins/alignments/src/CramAdapter/util.test.ts
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
import { readFeaturesToMismatches } from './util'

test('cram read features', () => {
expect(
readFeaturesToMismatches([{ code: 'i', pos: 5, refPos: 5, data: 'x' }], 0),
).toMatchSnapshot()
})
120 changes: 120 additions & 0 deletions plugins/alignments/src/CramAdapter/util.ts
Original file line number Diff line number Diff line change
@@ -0,0 +1,120 @@
import { CramRecord } from '@gmod/cram'

type ReadFeatures = CramRecord['readFeatures']

export interface Mismatch {
qual?: number
start: number
length: number
type: string
base: string | undefined
altbase?: string
seq?: string
cliplen?: number
}

export function readFeaturesToMismatches(
readFeatures: ReadFeatures,
start: number,
qual?: number[] | null,
) {
if (!readFeatures) {
return []
}
const mismatches: Mismatch[] = new Array(readFeatures.length)
let j = 0
let insLen = 0

let refPos = 0
let sublen = 0
let lastPos = start

for (let i = 0; i < readFeatures.length; i++) {
const f = readFeatures[i]
const { code, pos, data, sub, ref } = f
sublen = refPos - lastPos
lastPos = refPos

if (sublen && insLen > 0) {
mismatches[j++] = {
start: refPos,
type: 'insertion',
base: `${insLen}`,
length: 0,
}
insLen = 0
}
refPos = f.refPos - 1 - start

if (code === 'X') {
// substitution
mismatches[j++] = {
start: refPos,
length: 1,
base: sub,
qual: qual?.[pos],
altbase: ref,
type: 'mismatch',
}
} else if (code === 'I') {
// insertion
mismatches[j++] = {
start: refPos,
type: 'insertion',
base: `${data.length}`,
length: 0,
}
} else if (code === 'N') {
// reference skip
mismatches[j++] = {
type: 'skip',
length: data,
start: refPos,
base: 'N',
}
} else if (code === 'S') {
// soft clip
const len = data.length
mismatches[j++] = {
start: refPos,
type: 'softclip',
base: `S${len}`,
cliplen: len,
length: 1,
}
} else if (code === 'P') {
// padding
} else if (code === 'H') {
// hard clip
const len = data
mismatches[j++] = {
start: refPos,
type: 'hardclip',
base: `H${len}`,
cliplen: len,
length: 1,
}
} else if (code === 'D') {
// deletion
mismatches[j++] = {
type: 'deletion',
length: data,
start: refPos,
base: '*',
}
} else if (code === 'b') {
// stretch of bases
} else if (code === 'q') {
// stretch of qual scores
} else if (code === 'B') {
// a pair of [base, qual]
} else if (code === 'i') {
// single-base insertion
// insertion
insLen++
} else if (code === 'Q') {
// single quality value
}
}
return mismatches.slice(0, j)
}

0 comments on commit 3431a85

Please sign in to comment.