diff --git a/ChangeLog b/ChangeLog index 3caa4217..29154675 100644 --- a/ChangeLog +++ b/ChangeLog @@ -14,6 +14,10 @@ ## Interface and handling +- report a representative minimum-free-energy structure and its energy terms + for every ensemble-selected site; expose and sort by the retained `Eens` + value in native and browser output (#319) + - select and visualize the first rank-sorted CopomuS mutation automatically when browser results arrive (#317) @@ -348,6 +352,23 @@ energies, restricted partition sums, and trackers. ################################################################################ ################################################################################ +261104 Alexander Mitrofanov + * src/IntaRNA/{Interaction,OutputConstraint,OutputHandlerMfe,OutputHandlerCsv, + OutputHandlerText,OutputHandlerSeedCandidates}.{h,cpp}, src/bin/ + * Retain the ensemble energy used for model-P ranking and trace the + representative minimum-free-energy interaction only after final output + selection, including seeded and explicit-seed predictions. + + Add the Eens CSV column, make model-P CSV output include and sort by it by + default, and report both MFE and ensemble energies in text output. + * gui/{app.js,columns.js,options.js,parameter-file.js,protocol.js,results.js}, + gui/tests/, scripts/wasm/test-browser.py + * Offer Eens only for ensemble predictions, use it for default browser + sorting and pair-relative energy differences, and preserve the MFE energy + and structure columns for visualization. + * README.md, tests/ + * Document and test ensemble-site ranking with final MFE traceback (issue + #319). + 261027 Alexander Mitrofanov * gui/copomus.js, scripts/wasm/test-browser.py * Select the first rank-sorted CopomuS mutation after populating the browser diff --git a/README.md b/README.md index 8299dd0d..9da2d57d 100644 --- a/README.md +++ b/README.md @@ -666,12 +666,16 @@ for the given site `S` multiplied with the Boltzmann weights of the respective accessibility penalties `ED1` and `ED2`. Here, `R` denotes the gas constant and `T` the temperature of the system. At `T=37` degree Celsius, the product `RT` is about `RT=0.6163173043012`. -The site's partition function can be +The site's partition function can be used to compute the ensemble energy of all interactions of a given site via -`E(S) = -RT log(Z(S))` (`log` = natural logarithm), which is reported for the optimal site -(see [CSV output](#outModeCsv)). -This abstracts from individual inter-molecular base pairing and incorporates the -dynamics and flexibility of the interactions formed by two regions. +`Eens(S) = -RT log(Z(S))` (`log` = natural logarithm). IntaRNA selects and ranks +sites by `Eens`, then traces the minimum-free-energy interaction within each +finally reported site. Thus, the reported base pairs and `E` energy terms +describe a concrete representative structure, while the `Eens` CSV column and +the additional line in text output retain the site's ensemble energy. The MFE +traceback is performed only for interactions that reach final output (see +[CSV output](#outModeCsv)). This combines the dynamics and flexibility captured +by the site ensemble with a representative inter-molecular base-pair pattern. The overall partition function is given by `Zall = sum Z(S)`, i.e. summation over all possible sites, @@ -1582,7 +1586,8 @@ are - `hybridDB` : hybrid in dot-bar notation (interactin sites only) - `hybridDBfull` : hybrid in dot-bar notation (full sequence length) - `bpList` : list of hybrid base pairs, e.g. '(4,3):(5,2):(7,1)' -- `E` : overall interaction energy = E_hybrid + ED1 + ED2 +- `E` : minimum free energy of the reported interaction structure = E_hybrid + ED1 + ED2 +- `Eens` : ensemble energy used to select and rank a `--model=P` interaction site - `ED1` : ED value of seq1 - `ED2` : ED value of seq2 - `Pu1` : probability to be accessible for seq1 @@ -1628,6 +1633,10 @@ to the RNAplfold values from which the ED values are derived from! Using `--outCsvCols '*'`, all available columns are added to the output. +For `--model=P`, the default CSV columns additionally contain `Eens`, and rows +are sorted by `Eens` unless `--outCsvSort` is set explicitly. For other models, +`Eens` is not available and is printed as `NAN` when requested. + Energies are provided in unit *kcal/mol*, probabilities in the interval [0,1]. Position annotations start indexing with 1 (or the values set via `--qIdxPos0` and `--tIdxPos0`). diff --git a/gui/app.js b/gui/app.js index 2e4eee9f..05e27ef8 100644 --- a/gui/app.js +++ b/gui/app.js @@ -157,9 +157,10 @@ csvBody.append(element('h3', 'CSV output columns', 'subheading')); csvBody.append(element('p', 'seq1 is the target; seq2 is the query. Hover over a column name for details. Select columns for the table and CSV download; table visibility can be changed after computation.', 'input-help')); const selectionTools = element('div', '', 'selection-tools'); csvBody.append(selectionTools); const selectColumns = (selection) => { - csvOrder = [...selection]; for (const input of csvInputs) input.checked = selection.includes(input.value); inputsChanged(); + const available = selection.filter((key) => !csvInputs.find((input) => input.value === key)?.disabled); + csvOrder = [...available]; for (const input of csvInputs) input.checked = available.includes(input.value); inputsChanged(); }; -for (const [name, selection] of [['Select all', () => CSV_CATALOG.map((c) => c.key)], ['Task defaults', () => defaultColumns(getTask())], ['Clear selection', () => []]]) { +for (const [name, selection] of [['Select all', () => CSV_CATALOG.map((c) => c.key)], ['Task defaults', () => defaultColumns(getTask(), byId('model').value)], ['Clear selection', () => []]]) { const button = element('button', name); button.type = 'button'; button.addEventListener('click', () => selectColumns(selection())); selectionTools.append(button); } @@ -170,7 +171,7 @@ for (const group of CSV_GROUPS) { const button = element('button', name, 'text-button'); button.type = 'button'; button.setAttribute('aria-label', `${name} ${group.label} columns`); button.addEventListener('click', () => { - for (const input of csvInputs.filter((input) => group.keys.includes(input.value))) input.checked = checked; + for (const input of csvInputs.filter((input) => group.keys.includes(input.value) && !input.disabled)) input.checked = checked; inputsChanged(); }); legend.append(button); } @@ -203,7 +204,7 @@ function updateLength(name) { byId(`${name}-length`).textContent = `${parsed.records.length} record${parsed.records.length === 1 ? '' : 's'}, ${parsed.records.reduce((sum, r) => sum + r.sequence.length, 0).toLocaleString()} nucleotides`; } catch { byId(`${name}-length`).textContent = 'Check sequence format and length'; } } -function updateAvailability() { +function updateAvailability({ preserveCsv = false } = {}) { const evaluation = getTask() === 'evaluation', mutation = getTask() === 'mutation'; for (const button of byId('example-options').querySelectorAll('[data-example-task]')) { button.hidden = button.dataset.exampleTask !== getTask(); @@ -239,6 +240,25 @@ function updateAvailability() { // Upload controls follow the text area they populate, including async reads. for (const upload of form.querySelectorAll('input[data-for]')) upload.disabled = byId(upload.dataset.for).disabled; for (const input of outputInputs) { input.disabled = mutation || evaluation && EXTRA_OUTPUTS.find((entry) => entry.key === input.value).prediction; if (input.disabled) input.checked = false; } + // A site ensemble energy exists only for model-P predictions. On first + // entering that model, select and sort by Eens by default; later user edits + // remain untouched until the model becomes unavailable again. + const eens = csvInputs.find((input) => input.value === 'Eens'); + const ensemble = !evaluation && !mutation && byId('model').value === 'P'; + const becomingAvailable = eens.disabled && ensemble; + eens.disabled = !ensemble; + eens.parentElement.hidden = !ensemble; + const eensSort = [...byId('outCsvSort').options].find((option) => option.value === 'Eens'); + eensSort.disabled = !ensemble; + if (!ensemble) { + eens.checked = false; + csvOrder = csvOrder.filter((key) => key !== 'Eens'); + if (byId('outCsvSort').value === 'Eens') byId('outCsvSort').value = ''; + } else if (becomingAvailable && !preserveCsv) { + eens.checked = true; + if (!csvOrder.includes('Eens')) csvOrder.splice(csvOrder.indexOf('E')+1,0,'Eens'); + byId('outCsvSort').value = 'Eens'; + } const explicitMutations = !!byId('copomus-mutation-encoding').value.trim(); for (const input of [byId('copomus-candidate-selection'), byId('copomus-generator'), ...copomusFilterInputs]) input.disabled = !mutation || explicitMutations; for (const input of [...copomusMeasureInputs, byId('copomus-alpha'), byId('copomus-beta'), byId('copomus-mutation-encoding')]) input.disabled = !mutation; @@ -251,7 +271,7 @@ function resetParameters({ clearInputs = false } = {}) { for (const input of form.querySelectorAll('input[type=file]')) input.value = ''; byId('examples-menu').open = false; activateGroup('sequence-input'); } - csvOrder = getTask() === 'mutation' ? CSV_COLUMNS.split(',') : defaultColumns(getTask()); + csvOrder = getTask() === 'mutation' ? CSV_COLUMNS.split(',') : defaultColumns(getTask(), UI_DEFAULTS.model); byId('seed-mode').value = 'identify'; for (const field of FIELDS) { const input = byId(field.key); if (field.type === 'checkbox') input.checked = UI_DEFAULTS[field.key]; else input.value = UI_DEFAULTS[field.key]; } syncAccessibilityFlags(); @@ -306,7 +326,7 @@ function applyParameters(request, message = 'Parameters loaded. Compute to analy byId('copomus-generator').value = copomus.generator; byId('copomus-mutation-encoding').value = copomus.mutationEncoding; byId('copomus-alpha').value = copomus.alpha; byId('copomus-beta').value = copomus.beta; - updateAvailability(); activateGroup('sequence-input'); inputsChanged(message); + updateAvailability({ preserveCsv: true }); activateGroup('sequence-input'); inputsChanged(message); } for (const strand of ['query', 'target']) { const input = byId(strand), picker = byId(`${strand}-upload`); diff --git a/gui/columns.js b/gui/columns.js index 6849075a..d6359113 100644 --- a/gui/columns.js +++ b/gui/columns.js @@ -105,7 +105,13 @@ export const CSV_CATALOG = [ { "key": "E", "label": "Energy (kcal/mol)", - "description": "overall interaction energy", + "description": "minimum free energy of the reported interaction structure", + "numeric": true + }, + { + "key": "Eens", + "label": "Ensemble energy (kcal/mol)", + "description": "ensemble energy used to select and rank a model-P interaction site", "numeric": true }, { @@ -317,19 +323,19 @@ export const CSV_CATALOG = [ // Every native column appears once. Ensemble energies and their components are // grouped with probabilities so users can see the cost/meaning of those outputs. export const CSV_GROUPS = [ - { label: 'Energy', keys: ['E', 'ED1', 'ED2', 'E_init', 'E_loops', 'E_dangleL', 'E_dangleR', 'E_endL', 'E_endR', 'E_hybrid', 'E_norm', 'E_hybridNorm', 'E_add', 'seedE', 'seedED1', 'seedED2'] }, + { label: 'Energy', keys: ['E', 'Eens', 'ED1', 'ED2', 'E_init', 'E_loops', 'E_dangleL', 'E_dangleR', 'E_endL', 'E_endR', 'E_hybrid', 'E_norm', 'E_hybridNorm', 'E_add', 'seedE', 'seedED1', 'seedED2'] }, { label: 'Output', keys: ['subseqDP', 'subseqDB', 'hybridDP', 'hybridDB', 'hybridDPfull', 'hybridDBfull', 'bpList'] }, { label: 'Sequence', keys: ['id1', 'id2', 'seq1', 'seq2', 'subseq1', 'subseq2'] }, { label: 'Interaction', keys: ['start1', 'end1', 'start2', 'end2', 'seedStart1', 'seedEnd1', 'seedStart2', 'seedEnd2'] }, { label: 'Partition Function', keys: ['Pu1', 'Pu2', 'seedPu1', 'seedPu2', 'w', 'Etotal', 'Eall', 'Eall1', 'Eall2', 'Zall', 'Zall1', 'Zall2', 'EallTotal', 'P_E', 'RT'] }, ]; -export function defaultColumns(task) { +export function defaultColumns(task, model = 'X') { const context = ['id1', 'start1', 'end1', 'id2', 'start2', 'end2', 'subseqDP', 'hybridDP', 'E']; // Seed energies remain available in evaluation, where native output is NAN. // List them explicitly: seedEnd1/2 are coordinates, not energy columns. const extra = task === 'evaluation' - ? [...CSV_CATALOG.filter(({ key }) => key.startsWith('E')).map(({ key }) => key), 'seedE', 'seedED1', 'seedED2'] + ? [...CSV_CATALOG.filter(({ key }) => key.startsWith('E') && key !== 'Eens').map(({ key }) => key), 'seedE', 'seedED1', 'seedED2'] : ['ED1', 'ED2', 'seedE', 'seedStart1', 'seedEnd1', 'seedStart2', 'seedEnd2']; - return [...new Set([...context, ...extra])]; + return [...new Set([...context, ...(task === 'prediction' && model === 'P' ? ['Eens'] : []), ...extra])]; } diff --git a/gui/options.js b/gui/options.js index 9f06d224..9c6fee4a 100644 --- a/gui/options.js +++ b/gui/options.js @@ -185,4 +185,4 @@ export const EXTRA_OUTPUTS = [ ['detailed', 'Detailed interaction report', false, 'Expanded text report with interaction structures, coordinates and energy contributions.'], ['ensemble', 'Ensemble report', false, 'Text report of interaction and intramolecular ensemble energies and partition functions.'], ].map(([key, label, predictionOnly, description]) => ({ key, label, prediction: predictionOnly, description })); -export const ENSEMBLE_COLUMNS = ['Eall', 'Zall', 'EallTotal', 'P_E']; +export const ENSEMBLE_COLUMNS = ['Eens', 'Eall', 'Zall', 'EallTotal', 'P_E']; diff --git a/gui/parameter-file.js b/gui/parameter-file.js index 376050e5..91263bb2 100644 --- a/gui/parameter-file.js +++ b/gui/parameter-file.js @@ -127,6 +127,8 @@ export function importParameters(source) { const normalized = normalizeParameters(request); request.qIdxPos0 = normalized.qIdxPos0; request.tIdxPos0 = normalized.tIdxPos0; request.parameters = normalized; + if (request.task === 'prediction' && normalized.model === 'P' + && !Object.hasOwn(document.parameters, 'outCsvCols')) request.csvColumns.push('Eens'); if (request.task === 'mutation') request.copomus = normalizeCopomus(request.copomus); return request; } diff --git a/gui/protocol.js b/gui/protocol.js index e5598e5c..b8388d5c 100644 --- a/gui/protocol.js +++ b/gui/protocol.js @@ -181,8 +181,10 @@ export function prepareRequest(request) { if (!/^(D(m-?\d+(\.\d+)?)?(b-?\d+(\.\d+)?)?|Z(b\d+(\.\d+)?)?|W)$/.test(parameters[`${p}ShapeMethod`])) throw new Error(`${p}ShapeMethod: invalid SHAPE method encoding.`); if (!/^(M|C(\d+(\.\d+)?)?|S|[LO](s-?\d+(\.\d+)?i-?\d+(\.\d+)?)?)$/.test(parameters[`${p}ShapeConversion`])) throw new Error(`${p}ShapeConversion: invalid SHAPE conversion encoding.`); } - const csvColumns = request.csvColumns ?? CSV_COLUMNS.split(','); + const csvColumns = request.csvColumns + ?? `${CSV_COLUMNS}${!evaluating && parameters.model === 'P' ? ',Eens' : ''}`.split(','); if (!Array.isArray(csvColumns) || !csvColumns.length || new Set(csvColumns).size !== csvColumns.length || csvColumns.some((key) => !CSV_CATALOG.some((c) => c.key === key))) throw new Error('Select at least one available CSV column, without duplicates.'); + if (csvColumns.includes('Eens') && (evaluating || parameters.model !== 'P')) throw new Error('Eens output is available only for ensemble model P predictions.'); if (parameters.outCsvSort && !csvColumns.includes(parameters.outCsvSort)) throw new Error('The CSV sort column must be selected for output.'); const outputs = request.outputs ?? []; if (!Array.isArray(outputs) || new Set(outputs).size !== outputs.length || outputs.some((key) => !EXTRA_OUTPUTS.some((entry) => entry.key === key && (!evaluating || !entry.prediction)))) throw new Error('Choose available additional output formats for this task.'); diff --git a/gui/results.js b/gui/results.js index 0029f18f..def67347 100644 --- a/gui/results.js +++ b/gui/results.js @@ -2,7 +2,10 @@ import { CSV_CATALOG } from './columns.js'; import { CONTEXT_COLUMNS } from './protocol.js'; import { InteractionVisualization } from './visualization.js'; -export const COLUMNS = Object.freeze([...CSV_CATALOG, { key: 'deltaE', label: 'mfe Δ (kcal/mol)', numeric: true }]); +export const COLUMNS = Object.freeze([...CSV_CATALOG, { + key: 'deltaE', label: 'Energy Δ (kcal/mol)', numeric: true, + description: 'Difference from the best ranking energy for this sequence pair (Eens for model P, otherwise E).', +}]); const catalog = new Map(COLUMNS.map((column) => [column.key, column])); // Native CSV is unquoted, semicolon-separated; unavailable optional fields use @@ -28,11 +31,12 @@ export function parseResults(csv) { }); // Each sequence pair has its own optimum in a multi-FASTA screen. const best = new Map(); + const rankingEnergy = (row) => Number.isFinite(Number(row.Eens)) ? Number(row.Eens) : Number(row.E); for (const row of rows) { const pair = JSON.stringify([row.id1, row.id2]); - best.set(pair, Math.min(best.get(pair) ?? Infinity, row.E)); + best.set(pair, Math.min(best.get(pair) ?? Infinity, rankingEnergy(row))); } - return rows.map((row) => ({ ...row, deltaE: Number((row.E - best.get(JSON.stringify([row.id1, row.id2]))).toFixed(2)) })); + return rows.map((row) => ({ ...row, deltaE: Number((rankingEnergy(row) - best.get(JSON.stringify([row.id1, row.id2]))).toFixed(2)) })); } function numeric(value) { return String(value).trim() && Number.isFinite(Number(value)) ? Number(value) : null; } export function selectRows(rows, search, key, direction) { @@ -96,7 +100,7 @@ export class ResultsView { }); label.append(input, document.createTextNode(column.label)); this.picker.append(label); } - this.search.value = ''; this.key = columns.includes('E') ? 'E' : columns[0]; this.direction = 1; + this.search.value = ''; this.key = columns.includes('Eens') ? 'Eens' : columns.includes('E') ? 'E' : columns[0]; this.direction = 1; this.root.hidden = false; this.render(true); return this.rows.length; } render(selectFirst = false) { @@ -142,7 +146,7 @@ export class ResultsView { action.append(view); tr.append(action); for (const column of columns) { const cell = document.createElement('td'); - cell.textContent = ['E', 'deltaE'].includes(column.key) ? Number(row[column.key]).toFixed(2) : String(row[column.key]); + cell.textContent = ['E', 'Eens', 'deltaE'].includes(column.key) ? Number(row[column.key]).toFixed(2) : String(row[column.key]); if (column.numeric) cell.className = 'number'; tr.append(cell); } this.body.append(tr); diff --git a/gui/tests/parameter-file.test.mjs b/gui/tests/parameter-file.test.mjs index 4bbe922e..ff673d8f 100644 --- a/gui/tests/parameter-file.test.mjs +++ b/gui/tests/parameter-file.test.mjs @@ -43,6 +43,8 @@ test('invalid/native-only settings are rejected before replacing the form', () = assert.throws(() => importParameters('{invalid')); assert.throws(() => importParameters(' '.repeat(5000001)), /5 MB/); assert.deepEqual(importParameters('\ufeff'+JSON.stringify(base)).csvColumns, CSV_COLUMNS.split(',')); + assert.deepEqual(importParameters({ ...base, parameters: { ...base.parameters, model: 'P' } }).csvColumns, + [...CSV_COLUMNS.split(','), 'Eens']); }); test('maintainer examples use the same JSON import and prediction validation', async () => { diff --git a/gui/tests/protocol.test.mjs b/gui/tests/protocol.test.mjs index 14dbb8f8..504639fa 100644 --- a/gui/tests/protocol.test.mjs +++ b/gui/tests/protocol.test.mjs @@ -36,6 +36,10 @@ test('CLI defaults, ranges and strict typed controls reject invalid values', () assert.equal(sequenceOffset(-1, -3), 2); assert.throws(() => prepareRequest(request({}, { csvColumns: [] })), /at least one/); assert.throws(() => prepareRequest(request({}, { outputs: ['../other'] })), /available/); + const ensemble = prepareRequest(request({ model: 'P' })); + assert.ok(ensemble.csvColumns.includes('Eens')); + assert.ok(makeArguments(ensemble).some((arg) => arg.includes('Eens'))); + assert.throws(() => prepareRequest(request({ model: 'X' }, { csvColumns: [...CSV_COLUMNS.split(','), 'Eens'] })), /model P/); }); test('evaluation validates antiparallel supplied structures in signed coordinates', () => { @@ -135,10 +139,12 @@ test('browser batches use longest selected strands and round-trip full FASTA inp test('CSV groups cover native columns once and task defaults serialize explicitly', async () => { assert.deepEqual(CSV_GROUPS.flatMap((group) => group.keys).toSorted(), CSV_CATALOG.map((column) => column.key).toSorted()); - const energy = [...CSV_GROUPS.find((group) => group.label === 'Energy').keys, 'Etotal', 'Eall', 'Eall1', 'Eall2', 'EallTotal']; + const energy = [...CSV_GROUPS.find((group) => group.label === 'Energy').keys.filter((key) => key !== 'Eens'), 'Etotal', 'Eall', 'Eall1', 'Eall2', 'EallTotal']; for (const key of energy) assert.ok(defaultColumns('evaluation').includes(key), key); for (const key of ['seedStart1', 'seedStart2', 'seedEnd1', 'seedEnd2']) assert.ok(!defaultColumns('evaluation').includes(key), key); for (const key of ['ED1', 'ED2', 'seedE', 'seedStart1', 'seedEnd1', 'seedStart2', 'seedEnd2']) assert.ok(defaultColumns('prediction').includes(key), key); + assert.ok(defaultColumns('prediction', 'P').includes('Eens')); + assert.ok(!defaultColumns('prediction', 'X').includes('Eens')); const { exportParameters, importParameters } = await import('../parameter-file.js'); const input = request({ outNumber: 10, outNumberPerSeed: 2 }, { csvColumns: defaultColumns('prediction') }); const saved = exportParameters(input, 'test'); diff --git a/gui/tests/results.test.mjs b/gui/tests/results.test.mjs index 5ea65311..ce7299e5 100644 --- a/gui/tests/results.test.mjs +++ b/gui/tests/results.test.mjs @@ -15,6 +15,15 @@ test('native CSV orientation, numeric coordinates, energy and relative suboptima assert.equal(selectRows(rows, 'missing', 'id2', 1).length, 0); }); +test('ensemble output computes pair-relative differences and sorting from Eens', () => { + const csv = 'id1;start1;end1;id2;start2;end2;E;Eens\n' + + 'target;1;8;query;1;8;-10;-5.2\n' + + 'target;2;7;query;2;7;-12;-4.1\n'; + const rows = parseResults(csv); + assert.deepEqual(rows.map((row) => row.deltaE), [0, 1.1]); + assert.deepEqual(selectRows(rows, '', 'Eens', 1).map((row) => row.E), [-10, -12]); +}); + test('empty output and literal native ID quotes; malformed output fails clearly', () => { assert.deepEqual(parseResults(''), []); assert.deepEqual(parseResults(header), []); diff --git a/scripts/wasm/test-browser.py b/scripts/wasm/test-browser.py index e12d1066..0b583db9 100644 --- a/scripts/wasm/test-browser.py +++ b/scripts/wasm/test-browser.py @@ -223,9 +223,24 @@ def predict(): assert groups.count() == 5 for group in groups.all(): group.get_by_role('button', name=re.compile('^Select all')).click() - assert group.locator('input:not(:checked)').count() == 0 + assert group.locator('input:enabled:not(:checked)').count() == 0 group.get_by_role('button', name='Deselect all').click() assert group.locator('input:checked').count() == 0 + eens = page.locator('.csv-group input[value="Eens"]') + expect(eens).to_be_disabled() + expect(eens.locator('xpath=..')).not_to_be_visible() + open_group('interaction-options') + page.locator('#model').select_option('P') + open_group('csv-options') + expect(eens).to_be_enabled() + expect(eens.locator('xpath=..')).to_be_visible() + expect(eens).to_be_checked() + expect(page.locator('#outCsvSort')).to_have_value('Eens') + open_group('interaction-options') + page.locator('#model').select_option('X') + open_group('csv-options') + expect(eens).to_be_disabled() + expect(eens).not_to_be_checked() page.get_by_role('button', name='Task defaults', exact=True).click() for key in ['ED1', 'ED2', 'seedE', 'seedStart1', 'seedStart2']: expect(page.locator('.csv-group input[value="' + key + '"]')).to_be_checked() diff --git a/src/IntaRNA/Interaction.cpp b/src/IntaRNA/Interaction.cpp index 87f55695..cca12e83 100644 --- a/src/IntaRNA/Interaction.cpp +++ b/src/IntaRNA/Interaction.cpp @@ -171,6 +171,8 @@ operator= ( const Interaction & toCopy ) // copy energy value energy = toCopy.energy; + // copy optional ensemble ranking energy + ensembleEnergy = toCopy.ensembleEnergy; // copy seed data if (toCopy.seed != NULL) { @@ -216,6 +218,8 @@ operator= ( const InteractionRange & range ) // copy energy value energy = range.energy; + // ranges do not encode an ensemble ranking energy + ensembleEnergy = E_INF; return *this; } @@ -229,6 +233,8 @@ operator == ( const Interaction &i ) const return s1 == i.s1 && s2 == i.s2 && E_equal( energy, i.energy ) + && ((E_isINF(ensembleEnergy) && E_isINF(i.ensembleEnergy)) + || E_equal(ensembleEnergy, i.ensembleEnergy)) && basePairs == i.basePairs && (seed == i.seed || (seed != NULL && i.seed != NULL && *seed == *(i.seed))) diff --git a/src/IntaRNA/Interaction.h b/src/IntaRNA/Interaction.h index b1806d3e..d8993af7 100644 --- a/src/IntaRNA/Interaction.h +++ b/src/IntaRNA/Interaction.h @@ -231,6 +231,9 @@ class Interaction { //! energy of the interaction E_type energy; + //! optional ensemble energy used to rank an ensemble-selected interaction + E_type ensembleEnergy; + //! optional: seed information SeedSet * seed; @@ -346,11 +349,10 @@ class Interaction { /** * Checks whether or not this interaction is considered better (smaller) - * than another interaction of the same sequences + * than another interaction of the same sequences. The ensemble energy is + * used when available, otherwise the interaction energy is used. * @param i the interaction to compare to (for the same sequences!) - * @return energy < i.energy - * || (energy == i.energy && bpLeft < i.bpLeft) - * || (energy == i.energy && bpLeft == i.bpLeft && bpRight < i.bpRight) + * @return rankingEnergy < i.rankingEnergy, with boundary tie breaking */ bool operator < ( const Interaction &i ) const; @@ -366,12 +368,28 @@ class Interaction { * * @param energy the energy to compare to * @param hasLargerE the interaction that is supposed to have larger energy - * @return (energy < (hasLargerE.energy-precisionDelta)) + * @return (energy < (hasLargerE.rankingEnergy-precisionDelta)) */ static bool compareEnergy( const E_type& energy, const Interaction & hasLargerE ); + /** + * Whether this interaction was selected by an ensemble-based predictor. + * @return true if ensembleEnergy is available + */ + bool + hasEnsembleEnergy() const; + + /** + * Energy used for ranking and output filtering. Ensemble-selected + * interactions retain their ensemble energy after energy is replaced by the + * representative structure's MFE. + * @return ensembleEnergy when available, energy otherwise + */ + E_type + getRankingEnergy() const; + /** * Prints the interacting base pair to stream * @param out the ostream to write to @@ -468,6 +486,7 @@ Interaction::Interaction( const RnaSequence & s1, const RnaSequence & s2 ) , s2(&s2) , basePairs() , energy( E_INF ) + , ensembleEnergy( E_INF ) , seed( NULL ) { } @@ -481,6 +500,7 @@ Interaction::Interaction( const Interaction & toCopy ) , s2(toCopy.s2) , basePairs(toCopy.basePairs) , energy( toCopy.energy ) + , ensembleEnergy( toCopy.ensembleEnergy ) , seed( toCopy.seed == NULL ? NULL : new SeedSet() ) { // copy seed if needed @@ -499,6 +519,7 @@ Interaction::Interaction( const InteractionRange & range ) , s2(NULL) , basePairs() , energy( E_INF ) + , ensembleEnergy( E_INF ) , seed( NULL ) { // init data @@ -545,6 +566,8 @@ clear() basePairs.clear(); // clear energy energy = E_INF; + // clear optional ensemble ranking energy + ensembleEnergy = E_INF; // undo seed information INTARNA_CLEANUP(seed); @@ -557,7 +580,28 @@ bool Interaction:: compareEnergy( const E_type & energy, const Interaction & hasLargerE ) { - return energy < hasLargerE.energy && !(E_equal(energy,hasLargerE.energy)); + return energy < hasLargerE.getRankingEnergy() + && !(E_equal(energy,hasLargerE.getRankingEnergy())); +} + +//////////////////////////////////////////////////////////////////////////// + +inline +bool +Interaction:: +hasEnsembleEnergy() const +{ + return E_isNotINF(ensembleEnergy); +} + +//////////////////////////////////////////////////////////////////////////// + +inline +E_type +Interaction:: +getRankingEnergy() const +{ + return hasEnsembleEnergy() ? ensembleEnergy : energy; } //////////////////////////////////////////////////////////////////////////// @@ -572,11 +616,13 @@ operator < ( const Interaction &i ) const if (basePairs.empty() || i.basePairs.empty()) throw std::runtime_error("Interaction::operator < () : comparing empty interactions"); #endif // check if energy is larger - if (energy > i.energy && !(E_equal(energy, i.energy))) { + const E_type rankEnergy = getRankingEnergy(); + const E_type otherRankEnergy = i.getRankingEnergy(); + if (rankEnergy > otherRankEnergy && !(E_equal(rankEnergy, otherRankEnergy))) { return false; } // check if energy is equal - if (E_equal(energy, i.energy)) { + if (E_equal(rankEnergy, otherRankEnergy)) { return // i1 is smaller basePairs.begin()->first < i.basePairs.begin()->first // i1 diff --git a/src/IntaRNA/Makefile.am b/src/IntaRNA/Makefile.am index e061f0d2..b2f1bd4a 100644 --- a/src/IntaRNA/Makefile.am +++ b/src/IntaRNA/Makefile.am @@ -61,6 +61,7 @@ libIntaRNA_a_HEADERS = \ OutputHandlerEnsemble.h \ OutputHandlerHub.h \ OutputHandlerInteractionList.h \ + OutputHandlerMfe.h \ OutputHandlerSeedCandidates.h \ OutputHandlerText.h \ PredictionTracker.h \ @@ -133,6 +134,7 @@ libIntaRNA_a_SOURCES = \ OutputHandlerEnsemble.cpp \ OutputHandlerCsv.cpp \ OutputHandlerInteractionList.cpp \ + OutputHandlerMfe.cpp \ OutputHandlerSeedCandidates.cpp \ OutputHandlerText.cpp \ PredictionTrackerPairMinE.cpp \ diff --git a/src/IntaRNA/OutputConstraint.h b/src/IntaRNA/OutputConstraint.h index c48e63d0..4d29fadf 100644 --- a/src/IntaRNA/OutputConstraint.h +++ b/src/IntaRNA/OutputConstraint.h @@ -38,10 +38,10 @@ class OutputConstraint //! defines whether and where overlapping interaction sites are allowed for reporting const ReportOverlap reportOverlap; - //! upper bound (exclusive) for the energy of a reported interaction (E(interaction) < maxE) + //! upper bound (exclusive) for the ranking energy of a reported interaction const E_type maxE; - //! the maximal energy difference to the mfe of a reported interaction + //! maximal ranking-energy difference to the best reported interaction const E_type deltaE; //! whether or not only the best or all putative seeds are to be reported @@ -74,8 +74,8 @@ class OutputConstraint * reported to the output handler * @param reportOverlap defines whether and where overlapping interaction * sites are allowed for reporting - * @param maxE upper bound (exclusive) for the energy of a reported interaction (E(interaction) < maxE) - * @param deltaE maximal energy difference of a reported interaction to mfe + * @param maxE upper bound (exclusive) for a reported interaction's ranking energy + * @param deltaE maximal ranking-energy difference to the best interaction * @param bestSeedOnly whether or not only the best seed is to be reported * @param noLP whether or not lonely (non-stacked) inter-molecular bps are allowed * @param noGUend whether or not inter-molecular UG base pairs are allowed at interaction ends diff --git a/src/IntaRNA/OutputHandlerCsv.cpp b/src/IntaRNA/OutputHandlerCsv.cpp index 34d04a7e..f438c142 100644 --- a/src/IntaRNA/OutputHandlerCsv.cpp +++ b/src/IntaRNA/OutputHandlerCsv.cpp @@ -22,7 +22,7 @@ const std::string OutputHandlerCsv::notAvailable = "NAN"; const OutputHandlerCsv::ColTypeList OutputHandlerCsv::colTypeNumericSort( OutputHandlerCsv::string2list( "start1,end1,start2,end2" - ",E,ED1,ED2,Pu1,Pu2,E_init,E_loops,E_dangleL,E_dangleR,E_endL,E_endR,E_hybrid,E_norm,E_hybridNorm,E_add" + ",E,Eens,ED1,ED2,Pu1,Pu2,E_init,E_loops,E_dangleL,E_dangleR,E_endL,E_endR,E_hybrid,E_norm,E_hybridNorm,E_add" ",w" ",seedStart1,seedEnd1,seedStart2,seedEnd2,seedE,seedED1,seedED2,seedPu1,seedPu2" ",Eall,Eall1,Eall2" @@ -204,6 +204,11 @@ add( const Interaction & i ) outTmp < +#include +#include +#include +#include + +namespace IntaRNA { +namespace { + +const OutputConstraint & requireConstraint(OutputHandler * output) +{ + if (output == NULL) throw std::invalid_argument("OutputHandlerMfe: output is NULL"); + return output->getOutputConstraint(); +} + +/** Expose fixed-boundary traceback without running a complete prediction. */ +class FixedBoundaryMfe : public PredictorMfe2d { +public: + using PredictorMfe2d::PredictorMfe2d; + bool trace(Interaction & interaction); +}; + +bool +FixedBoundaryMfe::trace(Interaction & interaction) +{ + interaction.sort(); + const auto first = interaction.basePairs.front(); + const auto last = interaction.basePairs.back(); + energy.setOffset1(first.first); + energy.setOffset2(energy.getAccessibility2().getReversedIndex(first.second)); + const size_t n = last.first-first.first+1; + const size_t m = first.second-last.second+1; + hybridE_pq.resize(n,m); + + // A singleton has no interior recurrence but still has a complete energy. + if (n == 1 && m == 1) { + interaction.basePairs.resize(1); + interaction.energy = energy.getE(0,0,0,0,energy.getE_init()); + return E_isNotINF(interaction.energy); + } + + fillHybridE(n-1,m-1,0,0,false); + if (E_isINF(hybridE_pq(0,0))) return false; + PredictorMfe2d::traceBack(interaction); + interaction.energy = energy.getE(0,n-1,0,m-1,hybridE_pq(0,0)); + return E_isNotINF(interaction.energy); +} + +/** Fixed-boundary MFE traceback constrained to contain an admitted seed. */ +class FixedBoundarySeededMfe : public PredictorMfe2dSeed { +public: + FixedBoundarySeededMfe(const InteractionEnergy & energy, OutputHandler & output, + SeedHandler * seedHandler); + bool trace(Interaction & interaction); +}; + +FixedBoundarySeededMfe::FixedBoundarySeededMfe(const InteractionEnergy & energy, + OutputHandler & output, SeedHandler * seedHandler) + : PredictorMfe2dSeed(energy,output,NULL,seedHandler) +{} + +bool +FixedBoundarySeededMfe::trace(Interaction & interaction) +{ + interaction.sort(); + const auto first = interaction.basePairs.front(); + const auto last = interaction.basePairs.back(); + const size_t offset2 = energy.getAccessibility2().getReversedIndex(first.second); + energy.setOffset1(first.first); + energy.setOffset2(offset2); + seedHandler.setOffset1(first.first); + seedHandler.setOffset2(offset2); + const size_t n = last.first-first.first+1; + const size_t m = first.second-last.second+1; + hybridE_pq.resize(n,m); + hybridE_pq_seed.resize(n,m); + if (seedHandler.fillSeed(0,n-1,0,m-1) == 0) return false; + + // Explicit singleton seeds need no interior recurrence. + if (n == 1 && m == 1) { + if (!seedHandler.isSeedBound(0,0)) return false; + const E_type hybrid = seedHandler.getSeedE(0,0)+energy.getE_init(); + interaction.basePairs.resize(1); + interaction.energy = energy.getE(0,0,0,0,hybrid); + interaction.setSeedRange(interaction.basePairs.front(),interaction.basePairs.front(),interaction.energy); + return E_isNotINF(interaction.energy); + } + + fillHybridE(n-1,m-1,0,0,false); + if (E_isINF(hybridE_pq_seed(0,0))) return false; + interaction.energy = energy.getE(0,n-1,0,m-1,hybridE_pq_seed(0,0)); + PredictorMfe2dSeed::traceBack(interaction); + return E_isNotINF(interaction.energy); +} + +// Explicit seeds can contain singletons. The standard seeded MFE recurrence +// represents a singleton inside a larger interaction through an adjacent pair. +class ExplicitMfeSeedConstraint : public SeedConstraint { +public: + explicit ExplicitMfeSeedConstraint(const SeedConstraint & original); +}; + +ExplicitMfeSeedConstraint::ExplicitMfeSeedConstraint(const SeedConstraint & original) + : SeedConstraint(original) +{ + if (!getExplicitSeeds().empty()) bp = 2; +} + +} // namespace + +OutputHandlerMfe::OutputHandlerMfe(OutputHandler * output, + const InteractionEnergy & energy, SeedHandler * seedHandler) + : OutputHandler(requireConstraint(output)) + , output(output) + , energy(energy) + , seedHandler(seedHandler) +{} + +OutputHandlerMfe::~OutputHandlerMfe() = default; + +void +OutputHandlerMfe::traceMfe(Interaction & interaction) +{ + const OutputConstraint & constraint = getOutputConstraint(); + OutputHandlerInteractionList unused(constraint,1); + if (!seedHandler) { + FixedBoundaryMfe trace(energy,unused,NULL); + if (!trace.trace(interaction)) + throw std::runtime_error("ensemble output: no MFE structure in selected site"); + return; + } + + const auto & explicitSeeds = seedHandler->getConstraint().getExplicitSeeds(); + if (explicitSeeds.empty()) { + // The predictor owns its seed argument, so interpose a non-owning adapter. + FixedBoundarySeededMfe trace(energy,unused, + new SeedHandlerIdxOffset(seedHandler.get(),false)); + if (!trace.trace(interaction)) + throw std::runtime_error("ensemble output: no MFE structure in selected seeded site"); + return; + } + + // Trace every explicit alternative that occurs within the selected site and + // retain the lowest-energy representative. Mixed seed lengths and repeated + // starts therefore remain distinguishable. + ExplicitMfeSeedConstraint seedConstraint(seedHandler->getConstraint()); + std::vector variants; + const auto first = interaction.basePairs.front(); + const auto last = interaction.basePairs.back(); + if (first == last) { + FixedBoundarySeededMfe trace(energy,unused, + new SeedHandlerExplicit(energy,seedConstraint)); + if (!trace.trace(interaction)) + throw std::runtime_error("ensemble output: no MFE structure in selected explicit-seed site"); + return; + } + const size_t lo2 = energy.getAccessibility2().getReversedIndex(first.second); + const size_t hi2 = energy.getAccessibility2().getReversedIndex(last.second); + std::istringstream encodings(explicitSeeds); + for (std::string encoding; std::getline(encodings,encoding,','); ) { + seedConstraint.getExplicitSeeds() = encoding; + SeedHandlerExplicit seed(energy,seedConstraint); + size_t a = RnaSequence::lastPos, b = RnaSequence::lastPos; + if (!seed.updateToNextSeed(a,b,first.first,last.first,lo2,hi2)) continue; + if (seed.getSeedLength1(a,b) > 1) { + variants.push_back(encoding); + continue; + } + // Expand a singleton into each feasible adjacent two-pair seed. With + // noLP the neighbor must stack; otherwise any allowed interior loop works. + for (int direction : {-1,1}) { + const size_t max1 = constraint.noLP ? 1 : energy.getMaxInternalLoopSize1()+1; + const size_t max2 = constraint.noLP ? 1 : energy.getMaxInternalLoopSize2()+1; + for (size_t d1=1; d1<=max1; ++d1) for (size_t d2=1; d2<=max2; ++d2) { + if ((direction<0 && (d1>a-first.first || d2>b-lo2)) + || (direction>0 && (d1>last.first-a || d2>hi2-b))) continue; + const size_t c = direction<0 ? a-d1 : a+d1; + const size_t d = direction<0 ? b-d2 : b+d2; + if (!energy.areComplementary(c,d)) continue; + variants.push_back(std::to_string(interaction.s1->getInOutIndex(std::min(a,c))) + +"|"+std::string(d1-1,'.')+"|&" + +std::to_string(interaction.s2->getInOutIndex( + energy.getAccessibility2().getReversedIndex(std::max(b,d)))) + +"|"+std::string(d2-1,'.')+"|"); + } + } + } + + E_type best = E_INF; + Interaction representative(interaction); + for (const auto & encoding : variants) { + seedConstraint.getExplicitSeeds() = encoding; + FixedBoundarySeededMfe trace(energy,unused, + new SeedHandlerExplicit(energy,seedConstraint)); + Interaction candidate(interaction); + if (trace.trace(candidate) && candidate.energy < best) { + best = candidate.energy; + representative = candidate; + } + } + if (E_isINF(best)) + throw std::runtime_error("ensemble output: no MFE structure in selected explicit-seed site"); + interaction = representative; +} + +void +OutputHandlerMfe::add(const Interaction & interaction) +{ + if (interaction.isEmpty() || interaction.hasEnsembleEnergy()) { + output->add(interaction); + return; + } + Interaction representative(interaction); + representative.ensembleEnergy = interaction.energy; + // Ensemble predictors can already have traced seed pairs for output. Keep + // only the selected site boundaries so the fixed-boundary MFE traceback + // starts from the same canonical input for every predictor variant. + representative.sort(); + const Interaction::BasePair first = representative.basePairs.front(); + const Interaction::BasePair last = representative.basePairs.back(); + representative.basePairs.clear(); + representative.basePairs.push_back(first); + if (last != first) representative.basePairs.push_back(last); + INTARNA_CLEANUP(representative.seed); + traceMfe(representative); + output->add(representative); +} + +size_t +OutputHandlerMfe::reported() const +{ + return output->reported(); +} + +void +OutputHandlerMfe::incrementZ(const Z_type subZ) +{ + output->incrementZ(subZ); +} + +Z_type +OutputHandlerMfe::getZ() const +{ + return output->getZ(); +} + +void +OutputHandlerMfe::setExactPartition(bool exact) +{ + OutputHandler::setExactPartition(exact); + output->setExactPartition(exact); +} + +} // namespace IntaRNA diff --git a/src/IntaRNA/OutputHandlerMfe.h b/src/IntaRNA/OutputHandlerMfe.h new file mode 100644 index 00000000..ec4b8b6b --- /dev/null +++ b/src/IntaRNA/OutputHandlerMfe.h @@ -0,0 +1,67 @@ +#ifndef INTARNA_OUTPUTHANDLERMFE_H_ +#define INTARNA_OUTPUTHANDLERMFE_H_ + +#include "IntaRNA/InteractionEnergy.h" +#include "IntaRNA/OutputHandler.h" +#include "IntaRNA/SeedHandler.h" + +#include + +namespace IntaRNA { + +/** + * Final-output adapter for interactions selected by ensemble energy. + * + * The predictor-provided energy is retained as Interaction::ensembleEnergy. + * For each interaction that reaches this adapter, the minimum-free-energy + * structure with the same boundaries and seed requirements is traced and its + * energy is stored in Interaction::energy before forwarding. This keeps the + * additional MFE computation proportional to the final reported result count. + */ +class OutputHandlerMfe : public OutputHandler { +public: + /** + * Construction takes ownership of both supplied handlers. + * @param output final output handler receiving representative structures + * @param energy interaction energy model used for the ensemble prediction + * @param seedHandler optional seed handler; NULL denotes seed-free prediction + */ + OutputHandlerMfe(OutputHandler * output, const InteractionEnergy & energy, + SeedHandler * seedHandler = NULL); + + //! Destruction of the owned output and seed handlers. + virtual ~OutputHandlerMfe() override; + + /** + * Trace and forward one representative MFE interaction. + * @param interaction ensemble-selected interaction boundaries + */ + void add(const Interaction & interaction) override; + + //! @return number of interactions forwarded by the wrapped output handler + size_t reported() const override; + + //! Forward partition totals to the wrapped output handler. + void incrementZ(Z_type subZ) override; + + //! @return partition total of the wrapped output handler + Z_type getZ() const override; + + //! Forward exact-partition presentation policy. + void setExactPartition(bool exact) override; + +private: + //! final format-specific output handler + std::unique_ptr output; + //! unshifted energy model shared with the prediction + const InteractionEnergy & energy; + //! optional seed constraint used for fixed-boundary traceback + std::unique_ptr seedHandler; + + /** Fill interaction with the representative structure and MFE. */ + void traceMfe(Interaction & interaction); +}; + +} // namespace IntaRNA + +#endif /* INTARNA_OUTPUTHANDLERMFE_H_ */ diff --git a/src/IntaRNA/OutputHandlerSeedCandidates.cpp b/src/IntaRNA/OutputHandlerSeedCandidates.cpp index ba63a817..894c2f9a 100644 --- a/src/IntaRNA/OutputHandlerSeedCandidates.cpp +++ b/src/IntaRNA/OutputHandlerSeedCandidates.cpp @@ -31,7 +31,7 @@ void OutputHandlerSeedCandidates::add(const Interaction &) void OutputHandlerSeedCandidates::addSeed(const Interaction::BasePair & seed, const Interaction & interaction) { - if (interaction.isEmpty() || E_isINF(interaction.energy) || !outConstraint.reportMax) return; + if (interaction.isEmpty() || E_isINF(interaction.getRankingEnergy()) || !outConstraint.reportMax) return; // Window workers share this collector. Forward exceptions outside the block. std::exception_ptr failure; #if INTARNA_MULITHREADING @@ -74,10 +74,11 @@ void OutputHandlerSeedCandidates::report(OutputHandler & destination) const // Apply final restrictions once to the selected union, without quota refills. const auto & c = outConstraint; - const auto cutoff = std::int64_t(ranked.front().energy) + c.deltaE; + const auto cutoff = std::int64_t(ranked.front().getRankingEnergy()) + c.deltaE; std::vector reported; for (const auto & interaction : ranked) { - if (reported.size() >= c.reportMax || interaction.energy >= c.maxE || interaction.energy > cutoff) break; + if (reported.size() >= c.reportMax || interaction.getRankingEnergy() >= c.maxE + || interaction.getRankingEnergy() > cutoff) break; bool overlaps = false; for (const auto * previous : reported) { const bool target = interaction.basePairs.front().first <= previous->basePairs.back().first diff --git a/src/IntaRNA/OutputHandlerText.cpp b/src/IntaRNA/OutputHandlerText.cpp index fe1154f5..514e5e32 100644 --- a/src/IntaRNA/OutputHandlerText.cpp +++ b/src/IntaRNA/OutputHandlerText.cpp @@ -338,6 +338,9 @@ add( const Interaction & i ) <<"\n" <<"interaction energy = "< #include #include #include @@ -77,6 +78,7 @@ extern "C" { #include "IntaRNA/OutputHandlerCsv.h" #include "IntaRNA/OutputHandlerEnsemble.h" +#include "IntaRNA/OutputHandlerMfe.h" #include "IntaRNA/OutputHandlerText.h" @@ -89,6 +91,10 @@ const std::string CommandLineParsing::outCsvCols_default = "id1,start1,end1,id2, //////////////////////////////////////////////////////////////////////////// +const std::string CommandLineParsing::outCsvCols_ensemble_default = outCsvCols_default + ",Eens"; + +//////////////////////////////////////////////////////////////////////////// + const std::string CommandLineParsing::outCsvLstSep = ":"; //////////////////////////////////////////////////////////////////////////// @@ -1009,7 +1015,8 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) , value(&(outMaxE.val)) ->default_value(outMaxE.def) ->notifier(boost::bind(&CommandLineParsing::validate_numberArgument,this,outMaxE,_1)) - , std::string("only interactions with E < maxE are reported" + , std::string("only interactions with ranking energy < maxE are reported" + " (Eens for model P, E otherwise)" " (arg in range ["+toString(outMaxE.min)+","+toString(outMaxE.max)+"])").c_str()) (outMinPu.name.c_str() , value(&(outMinPu.val)) @@ -1021,7 +1028,8 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) , value(&(outDeltaE.val)) ->default_value(outDeltaE.def) ->notifier(boost::bind(&CommandLineParsing::validate_numberArgument,this,outDeltaE,_1)) - , std::string("suboptimal output : only interactions with E <= (minE+deltaE) are reported;" + , std::string("suboptimal output : only interactions with ranking energy <= (minimum+deltaE) are reported;" + " ranking energy is Eens for model P and E otherwise;" " minE is per sequence pair, or per region combination with --outPerRegion" " (arg in range ["+toString(outDeltaE.min)+","+toString(outDeltaE.max)+"])").c_str()) ("outBestSeedOnly" @@ -1046,12 +1054,13 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) , std::string("output : comma separated list of CSV column IDs to print if outMode=C." " Using '*' or an empty argument ('') prints all possible columns from the following available ID list: " + OutputHandlerCsv::list2string(OutputHandlerCsv::string2list(""),", ")+"." - + "\nDefault = '"+outCsvCols+"'." + + "\nDefault = '"+outCsvCols+"'; model P additionally reports Eens." ).c_str()) ("outCsvSort" , value(&(outCsvSort)) ->notifier(boost::bind(&CommandLineParsing::validate_outCsvSort,this,_1)) - , std::string("output : column ID from [outCsvCols] to be used for CSV row sorting if outMode=C.").c_str()) + , std::string("output : column ID from [outCsvCols] to be used for CSV row sorting if outMode=C;" + " defaults to Eens for model P.").c_str()) ("outPerRegion" , value(&outPerRegion) ->default_value(outPerRegion) @@ -1247,6 +1256,16 @@ parse(int argc, char** argv) // if parsing was successful, continue with final parsing if (parsingCode == ReturnCode::KEEP_GOING) { try { + // Ensemble prediction reports a representative MFE as E and keeps the + // site-selection objective in Eens. Apply these defaults only when the + // user did not provide an explicit CSV selection or sort override. + if (rri.empty() && model.val == 'P' && outMode.val == 'C') { + if (vm.at("outCsvCols").defaulted()) outCsvCols = outCsvCols_ensemble_default; + const auto columns = OutputHandlerCsv::string2list(outCsvCols); + if (vm.count("outCsvSort") == 0 + && std::find(columns.begin(),columns.end(),OutputHandlerCsv::Eens) != columns.end()) + outCsvSort = "Eens"; + } // parsing escape literals @@ -2846,27 +2865,44 @@ OutputHandler* CommandLineParsing:: getOutputHandler( const InteractionEnergy & energy ) const { + std::unique_ptr handler; switch (outMode.val) { case 'N' : - return new OutputHandlerText( getOutputConstraint(energy), outStreamHandler->getOutStream(), energy, 10, false ); + handler.reset(new OutputHandlerText( getOutputConstraint(energy), outStreamHandler->getOutStream(), energy, 10, false )); + break; case 'D' : - return new OutputHandlerText( getOutputConstraint(energy), outStreamHandler->getOutStream(), energy, 10, true ); + handler.reset(new OutputHandlerText( getOutputConstraint(energy), outStreamHandler->getOutStream(), energy, 10, true )); + break; case 'E' : // ensure that Zall is computed outNeedsZall = true; // no interaction details needed outNeedsBPs = false; - return new OutputHandlerEnsemble( getOutputConstraint(energy), outStreamHandler->getOutStream(), energy ); + handler.reset(new OutputHandlerEnsemble( getOutputConstraint(energy), outStreamHandler->getOutStream(), energy )); + break; case 'C' : // ensure that Zall is computed if needed outNeedsZall = outNeedsZall || OutputHandlerCsv::needsZall(OutputHandlerCsv::string2list( outCsvCols )); // check whether interaction details are needed outNeedsBPs = OutputHandlerCsv::needBPs(OutputHandlerCsv::string2list( outCsvCols ));; // create output handler - return new OutputHandlerCsv( getOutputConstraint(energy), outStreamHandler->getOutStream(), energy, OutputHandlerCsv::string2list( outCsvCols ), outSep, false, outCsvLstSep ); + handler.reset(new OutputHandlerCsv( getOutputConstraint(energy), outStreamHandler->getOutStream(), energy, OutputHandlerCsv::string2list( outCsvCols ), outSep, false, outCsvLstSep )); + break; default : INTARNA_NOT_IMPLEMENTED("Output mode "+toString(outMode.val)+" not implemented yet"); } + // Ensemble-only mode has no per-site structure. All other model-P formats + // receive one representative MFE traceback for each final reported site. + if (rri.empty() && model.val == 'P' && outMode.val != 'E') { + std::unique_ptr representativeSeed( + noSeedRequired ? NULL : getSeedHandler(energy)); + auto decorated = std::make_unique( + handler.get(),energy,representativeSeed.get()); + handler.release(); + representativeSeed.release(); + return decorated.release(); + } + return handler.release(); } //////////////////////////////////////////////////////////////////////////// diff --git a/src/bin/CommandLineParsing.h b/src/bin/CommandLineParsing.h index 266b14e6..989723ac 100644 --- a/src/bin/CommandLineParsing.h +++ b/src/bin/CommandLineParsing.h @@ -778,7 +778,7 @@ class CommandLineParsing { NumberParameter outNumberPerSeed; //! whether or not reported interactions can to be overlapping CharParameter outOverlap; - //! deltaE to mfe allowed to report an interaction + //! delta to the best ranking energy allowed for reported interactions NumberParameter outDeltaE; //! max E allowed to report an interaction NumberParameter outMaxE; @@ -800,6 +800,8 @@ class CommandLineParsing { std::string outCsvSort; //! the CSV column selection static const std::string outCsvCols_default; + //! default CSV columns for ensemble-based prediction + static const std::string outCsvCols_ensemble_default; //! whether or not best interaction output should be provided independently //! for all region combinations or only the best for each query-target //! combination diff --git a/src/bin/IntaRNA.cpp b/src/bin/IntaRNA.cpp index 8ae31d62..f91101ea 100644 --- a/src/bin/IntaRNA.cpp +++ b/src/bin/IntaRNA.cpp @@ -371,11 +371,11 @@ int main(int argc, char **argv){ // Independent per-region output retains each region's local window. const E_type maxMergedE = parameters.reportBestPerRegion() || bestInteractions.empty() ? E_INF - : (*bestInteractions.begin())->energy + bestInteractions.getOutputConstraint().deltaE; + : (*bestInteractions.begin())->getRankingEnergy() + bestInteractions.getOutputConstraint().deltaE; // The collector is sorted by energy, so all later entries are worse. for( const Interaction * inter : bestInteractions) { if (output->getOutputConstraint().reportMax==0) break; - if (inter->energy > maxMergedE) { break; } + if (inter->getRankingEnergy() > maxMergedE) { break; } output->add(*inter); } } diff --git a/tests/Makefile.am b/tests/Makefile.am index 72e6485f..ee20eec4 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -79,6 +79,7 @@ runApiTests_SOURCES = \ OutputStreamHandlerSortedCsv_test.cpp \ OutputHandlerHub_test.cpp \ OutputHandlerInteractionList_test.cpp \ + OutputHandlerMfe_test.cpp \ OutputHandlerSeedCandidates_test.cpp \ SeedHandlerExplicit_test.cpp \ SeedHandlerNoBulge_test.cpp \ diff --git a/tests/OutputHandlerMfe_test.cpp b/tests/OutputHandlerMfe_test.cpp new file mode 100644 index 00000000..22465419 --- /dev/null +++ b/tests/OutputHandlerMfe_test.cpp @@ -0,0 +1,76 @@ +#include "catch.hpp" + +#undef NDEBUG + +#include "IntaRNA/AccessibilityDisabled.h" +#include "IntaRNA/InteractionEnergyBasePair.h" +#include "IntaRNA/OutputHandlerCsv.h" +#include "IntaRNA/OutputHandlerInteractionList.h" +#include "IntaRNA/OutputHandlerMfe.h" +#include "IntaRNA/SeedHandlerExplicit.h" + +#include + +using namespace IntaRNA; + +TEST_CASE("OutputHandlerMfe", "[OutputHandlerMfe]") { + + // setup easylogging++ stuff if not already done + #include "testEasyLoggingSetup.icc" + + RnaSequence target("target","GGGG"); + RnaSequence query("query","CCCC"); + AccessibilityDisabled targetAcc(target,target.size(),NULL); + AccessibilityDisabled queryAcc(query,query.size(),NULL); + ReverseAccessibility queryAccReversed(queryAcc); + InteractionEnergyBasePair energy(targetAcc,queryAccReversed,16,16,false, + 1,Ekcal_2_E(-1.0),1); + OutputConstraint constraint(10,OutputConstraint::OVERLAP_BOTH, + Ekcal_2_E(100.0),E_INF,false,false,false,false,true); + + // Keep a raw observer while the adapter owns the final collector. + OutputHandlerInteractionList * stored = + new OutputHandlerInteractionList(constraint,10); + OutputHandlerMfe output(stored,energy); + Interaction site(target,query); + // Include predictor-traced interior pairs and seed metadata. The adapter must + // normalize these back to boundaries before its representative traceback. + for (size_t i=0; ibegin() != stored->end()); + const Interaction & representative = **stored->begin(); + REQUIRE(representative.basePairs.size() == 4); + REQUIRE(representative.seed == NULL); + REQUIRE(E_equal(representative.energy,Ekcal_2_E(-4.0))); + REQUIRE(E_equal(representative.ensembleEnergy,site.energy)); + REQUIRE(E_equal(representative.getRankingEnergy(),site.energy)); + + std::ostringstream csv; + OutputHandlerCsv csvOutput(constraint,csv,energy, + OutputHandlerCsv::string2list("E,Eens"),";",true); + csvOutput.add(representative); + REQUIRE(csv.str() == "E;Eens\n-4;-5.25\n"); + + SeedConstraint explicitConstraint(2,0,0,0,0, + AccessibilityDisabled::ED_UPPER_BOUND,0,IndexRangeList(""), + IndexRangeList(""),"2|&2|",false,false,false); + OutputHandlerInteractionList * explicitStored = + new OutputHandlerInteractionList(constraint,1); + OutputHandlerMfe explicitOutput(explicitStored,energy, + new SeedHandlerExplicit(energy,explicitConstraint)); + Interaction singleton(target,query); + singleton.basePairs.push_back(Interaction::BasePair(1,1)); + singleton.energy = Ekcal_2_E(-1.5); + explicitOutput.add(singleton); + REQUIRE(explicitStored->begin() != explicitStored->end()); + const Interaction & explicitRepresentative = **explicitStored->begin(); + REQUIRE(explicitRepresentative.basePairs.size() == 1); + REQUIRE(explicitRepresentative.seed != NULL); + REQUIRE(E_equal(explicitRepresentative.energy,Ekcal_2_E(-1.0))); + REQUIRE(E_equal(explicitRepresentative.ensembleEnergy,singleton.energy)); +}