diff --git a/app/scripts/nmr-cli/src/index.ts b/app/scripts/nmr-cli/src/index.ts index cc63fce..3d20466 100755 --- a/app/scripts/nmr-cli/src/index.ts +++ b/app/scripts/nmr-cli/src/index.ts @@ -25,9 +25,10 @@ Options for 'parse-spectra' command: -p, --path Directory path -s, --capture-snapshot Capture snapshot - + Options for 'predict' command: - -n, --nucleus Predicted nucleus "1H" or "13C" (required) + -ps,--peakShape Peak shape algorithm (default: "lorentzian") choices: ["gaussian", "lorentzian"] + -n, --nucleus Predicted nucleus, choices: ["1H","13C"] (required) -i, --id Input ID (default: 1) -t, --type NMR type (default: "nmr;1H;1d") -s, --shifts Chemical shifts (default: "1") @@ -35,10 +36,10 @@ Options for 'predict' command: -m, --molText MOL text (required) --from From in (ppm) --to To in (ppm) - --nbPoints Number of points (default: 128) + --nbPoints Number of points (default: 1024) --lineWidth Line width (default: 1) --frequency NMR frequency (MHz) (default: 400) - --tolerance Tolerance (default: 0.001) + --tolerance Tolerance to group peaks with close shift (default: 0.001) diff --git a/app/scripts/nmr-cli/src/prediction/generatePredictedSpectrumData.ts b/app/scripts/nmr-cli/src/prediction/generatePredictedSpectrumData.ts index ec8ca2c..78a8647 100644 --- a/app/scripts/nmr-cli/src/prediction/generatePredictedSpectrumData.ts +++ b/app/scripts/nmr-cli/src/prediction/generatePredictedSpectrumData.ts @@ -8,23 +8,21 @@ export interface ShiftsItem { spheres: number } -interface LorentzianOptions { +type PeakShape = 'gaussian' | 'lorentzian' + +interface PeakShapeOptions { x: number fwhm: number } -function lorentzian2(options: LorentzianOptions) { - const { x, fwhm } = options - return fwhm ** 2 / (4 * x ** 2 + fwhm ** 2) -} - -interface GenerateSpectrumOptions { +export interface GenerateSpectrumOptions { from?: number to?: number nbPoints?: number lineWidth?: number frequency?: number tolerance?: number + peakShape?: PeakShape } interface GroupItem { @@ -38,7 +36,37 @@ interface Data1D { re: number[] } -const getLorentzianFactor = (area = 0.9999) => { +const GAUSSIAN_EXP_FACTOR = -4 * Math.LN2 + +function lorentzian(options: PeakShapeOptions) { + const { x, fwhm } = options + return fwhm ** 2 / (4 * x ** 2 + fwhm ** 2) +} + +function gaussian(options: PeakShapeOptions) { + const { x, fwhm } = options + return Math.exp(GAUSSIAN_EXP_FACTOR * Math.pow(x / fwhm, 2)) +} + +function erfinv(x: number): number { + let a = 0.147 + if (x === 0) return 0 + let ln1MinusXSqrd = Math.log(1 - x * x) + let lnEtcBy2Plus2 = ln1MinusXSqrd / 2 + 2 / (Math.PI * a) + let firstSqrt = Math.sqrt(lnEtcBy2Plus2 ** 2 - ln1MinusXSqrd / a) + let secondSqrt = Math.sqrt(firstSqrt - lnEtcBy2Plus2) + return secondSqrt * (x > 0 ? 1 : -1) +} + +function getGaussianFactor(area = 0.9999) { + if (area >= 1) { + throw new Error('area should be (0 - 1)') + } + + return Math.sqrt(2) * erfinv(area) +} + +function getLorentzianFactor(area = 0.9999) { if (area >= 1) { throw new Error('area should be (0 - 1)') } @@ -49,6 +77,22 @@ const getLorentzianFactor = (area = 0.9999) => { ) } +function peakShapeFunction(options: PeakShapeOptions, shape: PeakShape) { + if (shape === 'lorentzian') { + return lorentzian(options) + } + + return gaussian(options) +} + +function getPeakShapeFactor(area: number, shape: PeakShape) { + if (shape === 'lorentzian') { + return getLorentzianFactor(area) + } + + return getGaussianFactor(area) +} + function groupEquivalentShifts(shifts: ShiftsItem[], tolerance = 0.001) { const groups: GroupItem[] = [] @@ -81,6 +125,7 @@ export function generatePredictedSpectrumData( frequency = 400, lineWidth = 1, tolerance = 0.001, + peakShape = 'lorentzian', } = options if (!shifts || shifts.length === 0) return [] @@ -107,24 +152,22 @@ export function generatePredictedSpectrumData( const stepSize = (to - from) / (nbPoints - 1) const groupedShifts = groupEquivalentShifts(acceptedShifts, tolerance) - console.log(groupedShifts) - - const limit = (lineWidth * getLorentzianFactor(0.99)) / frequency - console.log(limit) + const limit = (lineWidth * getPeakShapeFactor(0.99, peakShape)) / frequency for (let i = 0; i < nbPoints; i++) { const x = from + i * stepSize let intensity = 0 for (const { prediction, count } of groupedShifts) { if (Math.abs(x - prediction) <= limit) { intensity += - lorentzian2({ x: x - prediction, fwhm: lineWidth / frequency }) * - count + peakShapeFunction( + { x: x - prediction, fwhm: lineWidth / frequency }, + peakShape + ) * count } } - if (intensity > 0) { - data.x.push(x) - data.re.push(intensity) - } + + data.x.push(x) + data.re.push(intensity) } return data diff --git a/app/scripts/nmr-cli/src/prediction/parsePredictionCommand.ts b/app/scripts/nmr-cli/src/prediction/parsePredictionCommand.ts index 7200a9e..578c88c 100644 --- a/app/scripts/nmr-cli/src/prediction/parsePredictionCommand.ts +++ b/app/scripts/nmr-cli/src/prediction/parsePredictionCommand.ts @@ -1,247 +1,248 @@ import { Argv, CommandModule, Options } from 'yargs' import { - generatePredictedSpectrumData, - ShiftsItem, + generatePredictedSpectrumData, + GenerateSpectrumOptions, + ShiftsItem, } from './generatePredictedSpectrumData' import { v4 } from '@lukeed/uuid' import { CURRENT_EXPORT_VERSION } from 'nmr-load-save' import https from 'https' import axios from 'axios' -interface PredictionOptions { - from?: number - to?: number - nbPoints?: number - lineWidth?: number - frequency?: number - tolerance?: number -} - interface PredictionParameters { - molText: string - id: number - type: string - shifts: string - solvent: string - nucleus: string + molText: string + id: number + type: string + shifts: string + solvent: string + nucleus: string } -const predictionOptions: { [key in keyof PredictionOptions]: Options } = { - from: { - type: 'number', - description: 'From in (ppm)', - }, - to: { - type: 'number', - description: 'To in (ppm)', - }, - nbPoints: { - type: 'number', - description: 'Number of points', - default: 1024, - }, - lineWidth: { - type: 'number', - description: 'Line width', - default: 1, - }, - frequency: { - type: 'number', - description: 'NMR frequency (MHz)', - default: 400, - }, - tolerance: { - type: 'number', - description: 'Tolerance', - default: 0.001, - }, +const predictionOptions: { [key in keyof GenerateSpectrumOptions]: Options } = { + from: { + type: 'number', + description: 'From in (ppm)', + }, + to: { + type: 'number', + description: 'To in (ppm)', + }, + nbPoints: { + type: 'number', + description: 'Number of points', + default: 1024, + }, + lineWidth: { + type: 'number', + description: 'Line width', + default: 1, + }, + frequency: { + type: 'number', + description: 'NMR frequency (MHz)', + default: 400, + }, + tolerance: { + type: 'number', + description: 'Tolerance', + default: 0.001, + }, + peakShape: { + alias: 'ps', + type: 'string', + description: 'Peak shape algorithm', + default: 'lorentzian', + choices: ['gaussian', 'lorentzian'], + }, } as const const nmrOptions: { [key in keyof PredictionParameters]: Options } = { - id: { - alias: 'i', - type: 'number', - description: 'Input ID', - default: 1, - }, - type: { - alias: 't', - type: 'string', - description: 'NMR type', - default: 'nmr;1H;1d', - choices: ['nmr;1H;1d', 'nmr;13C;1d'], - }, - shifts: { - alias: 's', - type: 'string', - description: 'Chemical shifts', - default: '1', - }, - solvent: { - type: 'string', - description: 'NMR solvent', - default: 'Dimethylsulphoxide-D6 (DMSO-D6, C2D6SO)', - choices: [ - 'Any', - 'Chloroform-D1 (CDCl3)', - 'Dimethylsulphoxide-D6 (DMSO-D6, C2D6SO)', - 'Methanol-D4 (CD3OD)', - 'Deuteriumoxide (D2O)', - 'Acetone-D6 ((CD3)2CO)', - 'TETRACHLORO-METHANE (CCl4)', - 'Pyridin-D5 (C5D5N)', - 'Benzene-D6 (C6D6)', - 'neat', - 'Tetrahydrofuran-D8 (THF-D8, C4D4O)', - ], - }, - molText: { - alias: 'm', - type: 'string', - description: 'MOL file content', - requiresArg: true, - }, - nucleus: { - alias: 'n', - type: 'string', - description: 'Predicted nucleus', - requiresArg: true, - choices: ['1H', '13C'], - }, + id: { + alias: 'i', + type: 'number', + description: 'Input ID', + default: 1, + }, + type: { + alias: 't', + type: 'string', + description: 'NMR type', + default: 'nmr;1H;1d', + choices: ['nmr;1H;1d', 'nmr;13C;1d'], + }, + shifts: { + alias: 's', + type: 'string', + description: 'Chemical shifts', + default: '1', + }, + solvent: { + type: 'string', + description: 'NMR solvent', + default: 'Dimethylsulphoxide-D6 (DMSO-D6, C2D6SO)', + choices: [ + 'Any', + 'Chloroform-D1 (CDCl3)', + 'Dimethylsulphoxide-D6 (DMSO-D6, C2D6SO)', + 'Methanol-D4 (CD3OD)', + 'Deuteriumoxide (D2O)', + 'Acetone-D6 ((CD3)2CO)', + 'TETRACHLORO-METHANE (CCl4)', + 'Pyridin-D5 (C5D5N)', + 'Benzene-D6 (C6D6)', + 'neat', + 'Tetrahydrofuran-D8 (THF-D8, C4D4O)', + ], + }, + molText: { + alias: 'm', + type: 'string', + description: 'MOL file content', + requiresArg: true, + }, + nucleus: { + alias: 'n', + type: 'string', + description: 'Predicted nucleus', + requiresArg: true, + choices: ['1H', '13C'], + }, } as const interface PredictionResponseItem { - id: number - type: string - statistics: { - accept: number - warning: number - reject: number - missing: number - total: number - } - shifts: ShiftsItem[] + id: number + type: string + statistics: { + accept: number + warning: number + reject: number + missing: number + total: number + } + shifts: ShiftsItem[] } interface PredictionResponse { - result: PredictionResponseItem[] + result: PredictionResponseItem[] } async function predictNMR(options: PredictionArgs): Promise { - const url = process.env['NMR_PREDICTION_URL'] - - if (!url) { - throw new Error('Environment variable NMR_PREDICTION_URL is not defined.') - } - - try { - new URL(url).toString() - } catch { - throw new Error(`Invalid URL in NMR_PREDICTION_URL: "${url}"`) - } - - try { - const { - id, - type, - shifts, - solvent, - from, - to, - nbPoints = 1024, - frequency = 400, - lineWidth = 1, - tolerance = 0.001, - molText, - nucleus, - } = options - - const payload: any = { - inputs: [ - { - id, - type, - shifts, - solvent, - }, - ], - moltxt: molText.replaceAll(/\\n/g, '\n'), + const url = process.env['NMR_PREDICTION_URL'] + + if (!url) { + throw new Error('Environment variable NMR_PREDICTION_URL is not defined.') } - const httpsAgent = new https.Agent({ - rejectUnauthorized: false, - }) - - // Axios POST request with httpsAgent - const response = await axios.post(url, payload, { - headers: { - 'Content-Type': 'application/json', - }, - httpsAgent, - }) - - const responseResult: PredictionResponse = response.data - const spectra = [] - - for (const result of responseResult.result) { - const name = v4() - const data = generatePredictedSpectrumData(result.shifts, { - from, - to, - nbPoints, - lineWidth, - frequency, - tolerance, - }) - - const info = { - isFid: false, - isComplex: false, - dimension: 1, - originFrequency: frequency, - baseFrequency: frequency, - pulseSequence: '', - solvent, - isFt: true, - name, - nucleus, - } - - spectra.push({ - id: v4(), - data, - info, - }) + try { + new URL(url).toString() + } catch { + throw new Error(`Invalid URL in NMR_PREDICTION_URL: "${url}"`) } - const nmrium = { data: { spectra }, version: CURRENT_EXPORT_VERSION } - console.log(JSON.stringify(nmrium, null, 2)) - } catch (error) { - console.error( - 'Error:', - error instanceof Error ? error.message : String(error) - ) - - if (axios.isAxiosError(error) && error.response) { - console.error('Response data:', error.response.data) - } else if (error instanceof Error && error.cause) { - console.error('Network Error:', error.cause) + try { + const { + id, + type, + shifts, + solvent, + from, + to, + nbPoints = 1024, + frequency = 400, + lineWidth = 1, + tolerance = 0.001, + molText, + nucleus, + peakShape = "lorentzian", + } = options + + const payload: any = { + inputs: [ + { + id, + type, + shifts, + solvent, + }, + ], + moltxt: molText.replaceAll(/\\n/g, '\n'), + } + + const httpsAgent = new https.Agent({ + rejectUnauthorized: false, + }) + + // Axios POST request with httpsAgent + const response = await axios.post(url, payload, { + headers: { + 'Content-Type': 'application/json', + }, + httpsAgent, + }) + + const responseResult: PredictionResponse = response.data + const spectra = [] + + for (const result of responseResult.result) { + const name = v4() + const data = generatePredictedSpectrumData(result.shifts, { + from, + to, + nbPoints, + lineWidth, + frequency, + tolerance, + peakShape, + }) + + const info = { + isFid: false, + isComplex: false, + dimension: 1, + originFrequency: frequency, + baseFrequency: frequency, + pulseSequence: '', + solvent, + isFt: true, + name, + nucleus, + } + + spectra.push({ + id: v4(), + data, + info, + }) + } + + const nmrium = { data: { spectra }, version: CURRENT_EXPORT_VERSION } + console.log(JSON.stringify(nmrium, null, 2)) + } catch (error) { + console.error( + 'Error:', + error instanceof Error ? error.message : String(error) + ) + + if (axios.isAxiosError(error) && error.response) { + console.error('Response data:', error.response.data) + } else if (error instanceof Error && error.cause) { + console.error('Network Error:', error.cause) + } } - } } -type PredictionArgs = PredictionParameters & PredictionOptions +type PredictionArgs = PredictionParameters & GenerateSpectrumOptions // Define the prediction string command export const parsePredictionCommand: CommandModule<{}, PredictionArgs> = { - command: ['predict', 'p'], - describe: 'Predict NMR spectrum from mol text', - builder: (yargs: Argv<{}>): Argv => { - return yargs.options({ - ...nmrOptions, - ...predictionOptions, - }) as Argv - }, - handler: async argv => { - await predictNMR(argv) - }, + command: ['predict', 'p'], + describe: 'Predict NMR spectrum from mol text', + builder: (yargs: Argv<{}>): Argv => { + return yargs.options({ + ...nmrOptions, + ...predictionOptions, + }) as Argv + }, + handler: async argv => { + await predictNMR(argv) + }, }