MAPQ undefined being 255 handling

Author: cmdcolin

plugins/comparative-adapters/src/PAFAdapter/PAFAdapter.test.ts

@@ -3,6 +3,60 @@ import { toArray } from 'rxjs/operators'
 
 import Adapter from './PAFAdapter.ts'
 import MyConfigSchema from './configSchema.ts'
+import { getWeightedMeans } from './util.ts'
+
+import type { PAFRecord } from './util.ts'
+
+function makeRecord(
+  qname: string,
+  tname: string,
+  mappingQual: number | undefined,
+  blockLen: number,
+): PAFRecord {
+  return {
+    qname,
+    qstart: 0,
+    qend: blockLen,
+    tname,
+    tstart: 0,
+    tend: blockLen,
+    strand: 1,
+    extra: { mappingQual, blockLen },
+  }
+}
+
+test('getWeightedMeans treats MAPQ 255 as missing, not as quality 255', () => {
+  const records = [
+    makeRecord('q1', 't1', 255, 1000),
+    makeRecord('q1', 't1', 60, 1000),
+  ]
+  getWeightedMeans(records)
+  // If 255 were used as-is the mean would be (255+60)/2 = 157.5, so min=max=157.5, meanScore=0.5.
+  // With the fix, 255 → 1, so mean = (1+60)/2 = 30.5, still min=max=30.5, meanScore=0.5.
+  // The key assertion is that the record with MAPQ 255 doesn't skew the result toward 255.
+  expect(records[0]!.extra.meanScore).toBe(0.5)
+  expect(records[1]!.extra.meanScore).toBe(0.5)
+})
+
+test('getWeightedMeans returns 0.5 when all pairs have identical quality', () => {
+  const records = [
+    makeRecord('q1', 't1', 60, 1000),
+    makeRecord('q2', 't2', 60, 500),
+  ]
+  getWeightedMeans(records)
+  expect(records[0]!.extra.meanScore).toBe(0.5)
+  expect(records[1]!.extra.meanScore).toBe(0.5)
+})
+
+test('getWeightedMeans normalizes scores correctly across different quality pairs', () => {
+  const records = [
+    makeRecord('q1', 't1', 60, 1000),
+    makeRecord('q2', 't2', 20, 1000),
+  ]
+  getWeightedMeans(records)
+  expect(records[0]!.extra.meanScore).toBe(1)
+  expect(records[1]!.extra.meanScore).toBe(0)
+})
 
 function makeAdapter() {
   return new Adapter(

plugins/comparative-adapters/src/PAFAdapter/util.ts

@@ -59,7 +59,10 @@ export function getWeightedMeans(ret: PAFRecord[]) {
   const scoreMap: Record<string, { valueSum: number; weightSum: number }> = {}
   for (const entry of ret) {
     const key = `${entry.qname}-${entry.tname}`
-    const qual = entry.extra.mappingQual ?? 1
+    // MAPQ 255 is the PAF spec sentinel for "not computed" — treat as missing,
+    // not as a high-quality signal (which would massively upweight such records).
+    const mapq = entry.extra.mappingQual
+    const qual = mapq !== undefined && mapq !== 255 ? mapq : 1
     const len = entry.extra.blockLen ?? 1
     const existing = scoreMap[key]
     if (existing) {
@@ -86,7 +89,9 @@ export function getWeightedMeans(ret: PAFRecord[]) {
   for (const entry of ret) {
     const key = `${entry.qname}-${entry.tname}`
     const score = meanScoreMap[key]!
-    entry.extra.meanScore = range > 0 ? (score - min) / range : 1
+    // When all pairs have identical quality (range === 0), use 0.5 rather than
+    // 1 so the encoding reads as "neutral" instead of "maximum quality."
+    entry.extra.meanScore = range > 0 ? (score - min) / range : 0.5
   }
 
   return ret