isospec 0.1.0+2.6.0

/*
 *   Copyright (C) 2015-2026 Mateusz Łącki and Michał Startek.
 *
 *   This file is part of IsoSpec.
 *
 *   IsoSpec is free software: you can redistribute it and/or modify
 *   it under the terms of the Simplified ("2-clause") BSD licence.
 *
 *   IsoSpec is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 *
 *   You should have received a copy of the Simplified BSD Licence
 *   along with IsoSpec.  If not, see <https://opensource.org/licenses/BSD-2-Clause>.
 */

#include "fasta_mods.h"

#include <algorithm>
#include <cctype>
#include <cstdint>
#include <cstring>
#include <limits>
#include <stdexcept>
#include <string>
#include <string_view>
#include <vector>

#include "composition.h"
#include "element_lookup.h"
#include "element_tables.h"
#include "fasta.h"

namespace IsoSpec {

namespace {

//! Indices into a dense, direct-indexed accumulator: one slot per row of
//! element_tables.h's flat arrays (292), addressed by an element's *first*
//! isotope row -- every element this feature ever touches (the fixed 6 of
//! plain amino acids, plus the 27 Unimod's shipped compositions use) fits
//! this without any map. Sized at compile time, zero-initialized per call.
constexpr int kAccumulatorSize = ISOSPEC_NUMBER_OF_ISOTOPIC_ENTRIES;

//! First-index lookups for the 6 elements a plain amino-acid sequence can
//! ever contribute (fasta.h's own CHNOSSe order), resolved once via
//! element_lookup.h rather than duplicating fasta.h's separate, narrower
//! aa_elem_masses/aa_isotope_numbers tables -- this path folds the base
//! sequence into the same general per-element-row accumulator a
//! modification's composition delta uses, so both add into one array.
struct AminoAcidElementIndex {
    int first_index[6];
    AminoAcidElementIndex() {
        static const char* const kSymbols[6] = {"C", "H", "N", "O", "S", "Se"};
        for (int i = 0; i < 6; i++)
            first_index[i] = find_element_table_first_index(kSymbols[i], strlen(kSymbols[i]));
    }
};

const AminoAcidElementIndex& amino_acid_element_index() {
    static const AminoAcidElementIndex instance;
    return instance;
}

void add_residue(std::int64_t accumulator[kAccumulatorSize], char residue) {
    int counts[6] = {0, 0, 0, 0, 0, 0};
    parse_fasta(std::string_view(&residue, 1), counts);  // zero for anything that isn't a recognized amino-acid letter

    const AminoAcidElementIndex& idx = amino_acid_element_index();
    for (int i = 0; i < 6; i++)
        if (counts[i] != 0)
            accumulator[idx.first_index[i]] += counts[i];
}

//! `p` points at '[', `end` one past the sequence (which need not be
//! NUL-terminated: nothing here reads past `end`). Parses "UNIMOD:<id>]",
//! applies that id's composition delta to `accumulator`, and advances `p` past
//! the closing ']'. Throws std::invalid_argument on anything not exactly this
//! shape, or an id not present in `mods` -- unknown and deliberately-excluded
//! ids are the same error, no special-casing.
void consume_unimod_bracket(const char*& p, const char* const end, const UnimodTable& mods, std::int64_t accumulator[kAccumulatorSize]) {
    constexpr std::string_view kPrefix = "UNIMOD:";
    const auto rest = [p, end] { return std::string(p, end); };  // built only for an error message

    const char* q = p + 1;  // skip '['
    if (static_cast<size_t>(end - q) < kPrefix.size() || std::string_view(q, kPrefix.size()) != kPrefix)
        throw std::invalid_argument("Malformed modification bracket (expected '[UNIMOD:<id>]'): " + rest());
    q += kPrefix.size();

    const char* digits_start = q;
    while (q < end && isdigit(static_cast<unsigned char>(*q)))
        q++;
    if (q == digits_start)
        throw std::invalid_argument("Malformed modification bracket (missing numeric id): " + rest());
    if (q == end || *q != ']')
        throw std::invalid_argument("Malformed modification bracket (missing closing ']'): " + rest());

    uint64_t id;
    try {
        id = std::stoull(std::string(digits_start, static_cast<size_t>(q - digits_start)));
    } catch (const std::exception&) {
        throw std::invalid_argument("Malformed modification bracket (id out of range): " + rest());
    }

    // UnimodTable::lookup takes an unsigned int, so an id past UINT_MAX would
    // wrap on the way in and silently resolve to a completely different,
    // existing modification ("[UNIMOD:4294967300]" -> id 4 -> Carbamidomethyl)
    // instead of being rejected. An id that large is unknown to any table by
    // construction (Unimod's own id space is four digits); say so, rather
    // than quietly applying the wrong composition delta.
    const UnimodEntry* entry = (id > std::numeric_limits<unsigned int>::max())
                                   ? nullptr
                                   : mods.lookup(static_cast<unsigned int>(id));
    if (entry == nullptr)
        throw std::invalid_argument("Unknown or unsupported UNIMOD id in sequence: UNIMOD:" + std::to_string(id));

    for (size_t i = 0; i < entry->element_first_index.size(); i++)
        accumulator[entry->element_first_index[i]] += entry->element_delta_count[i];

    p = q + 1;  // past ']'
}

//! The accumulator is 64-bit because an int one overflows: an override
//! table's deltas are arbitrary ints ("H2000000000"), and a sequence may
//! apply them repeatedly. Clamped on the way out, like parse_fasta's counts
//! (AA_ATOM_COUNT_CEILING): too large either way, still rejected downstream.
void accumulator_to_composition_into(const std::int64_t accumulator[kAccumulatorSize], ElementComposition& out) {
    out.clear();
    for (int i = 0; i < kAccumulatorSize; i++) {
        if (accumulator[i] != 0) {
            out.element_first_index.push_back(i);
            out.count.push_back(static_cast<int>(std::clamp<std::int64_t>(
                accumulator[i], -AA_ATOM_COUNT_CEILING, AA_ATOM_COUNT_CEILING)));
        }
    }
}

}  // namespace

void parse_fasta_with_mods_into(const char* sequence, ElementComposition& out, const UnimodTable& mods) {
    parse_fasta_with_mods_into(std::string_view(sequence), out, mods);
}

void parse_fasta_with_mods_into(std::string_view sequence, ElementComposition& out, const UnimodTable& mods) {
    out.clear();

    if (sequence.find('[') == std::string_view::npos) {
        // '[' is the only character this parser ever treats specially (it's
        // never a valid FASTA/amino-acid character on its own), so its
        // absence means the sequence cannot contain a [UNIMOD:<id>]
        // reference -- skip the general per-element accumulator below
        // entirely and go through the plain, fixed-CHNOSSe parse_fasta in
        // one pass over the whole string, matching Iso::FromFASTA's cost
        // for this (the common, unmodified) case instead of paying for
        // per-residue re-parsing plus general-element bookkeeping nothing
        // here actually needs.
        int counts[6] = {0, 0, 0, 0, 0, 0};
        parse_fasta(sequence, counts);
        const AminoAcidElementIndex& idx = amino_acid_element_index();
        for (int i = 0; i < 6; i++) {
            if (counts[i] != 0) {
                out.element_first_index.push_back(idx.first_index[i]);
                out.count.push_back(counts[i]);
            }
        }
        return;
    }

    std::int64_t accumulator[kAccumulatorSize] = {0};

    const char* p = sequence.data();
    const char* const end = p + sequence.size();
    while (p < end) {
        if (*p == '[') {
            // A [UNIMOD:<id>] bracket has the same effect wherever it
            // appears (N-terminal "[id]-SEQUENCE", internal "X[id]",
            // C-terminal "SEQUENCE-[id]") -- the '-' either side of a
            // terminal bracket is just an ordinary character, handled by
            // the branch below like any other. No position-tracking needed.
            consume_unimod_bracket(p, end, mods, accumulator);
        } else {
            add_residue(accumulator, *p);
            p++;
        }
    }

    accumulator_to_composition_into(accumulator, out);
}

void parse_fasta_with_mods_full_into(const char* sequence, ElementComposition& out, const UnimodTable& mods) {
    parse_fasta_with_mods_full_into(std::string_view(sequence), out, mods);
}

void parse_fasta_with_mods_full_into(std::string_view sequence, ElementComposition& out, const UnimodTable& mods) {
    parse_fasta_with_mods_into(sequence, out, mods);

    const AminoAcidElementIndex& idx = amino_acid_element_index();
    const int h_index = idx.first_index[1];
    const int o_index = idx.first_index[3];

    bool has_h = false, has_o = false;
    for (size_t i = 0; i < out.element_first_index.size(); i++) {
        if (out.element_first_index[i] == h_index) {
            out.count[i] += 2;
            has_h = true;
        } else if (out.element_first_index[i] == o_index) {
            out.count[i] += 1;
            has_o = true;
        }
    }
    if (!has_h) {
        out.element_first_index.push_back(h_index);
        out.count.push_back(2);
    }
    if (!has_o) {
        out.element_first_index.push_back(o_index);
        out.count.push_back(1);
    }
}

ElementComposition parse_fasta_with_mods(std::string_view sequence, const UnimodTable& mods) {
    ElementComposition result;
    parse_fasta_with_mods_into(sequence, result, mods);
    return result;
}

ElementComposition parse_fasta_with_mods(const char* sequence, const UnimodTable& mods) {
    return parse_fasta_with_mods(std::string_view(sequence), mods);
}

ElementComposition parse_fasta_with_mods_full(std::string_view sequence, const UnimodTable& mods) {
    ElementComposition result;
    parse_fasta_with_mods_full_into(sequence, result, mods);
    return result;
}

ElementComposition parse_fasta_with_mods_full(const char* sequence, const UnimodTable& mods) {
    return parse_fasta_with_mods_full(std::string_view(sequence), mods);
}

Iso Iso::FromFASTAWithMods(const char* sequence, const UnimodTable& mods, bool use_nominal_masses, bool add_water) {
    return FromFASTAWithMods(std::string_view(sequence), mods, use_nominal_masses, add_water);
}

Iso Iso::FromFASTAWithMods(std::string_view sequence, const UnimodTable& mods, bool use_nominal_masses, bool add_water) {
    if (sequence.find('[') == std::string_view::npos)
        // Nothing for `mods` to resolve -- go straight through the plain
        // path instead of building a composition via the general
        // dimNumber-vector Iso constructor build_iso_from_composition uses
        // (heap-allocated mass/probability vectors, per-element isotope-
        // count lookups) when Iso::FromFASTA's own fixed static tables
        // already do this cheaper for the unmodified case.
        return Iso::FromFASTA(sequence, use_nominal_masses, add_water);

    ElementComposition composition =
        add_water ? parse_fasta_with_mods_full(sequence, mods) : parse_fasta_with_mods(sequence, mods);
    return build_iso_from_composition(composition, use_nominal_masses);
}

Iso Iso::FromFASTAWithMods(const char* sequence, bool use_nominal_masses, bool add_water, const char* unimod_db_path) {
    return FromFASTAWithMods(std::string_view(sequence), use_nominal_masses, add_water, unimod_db_path);
}

Iso Iso::FromFASTAWithMods(std::string_view sequence, bool use_nominal_masses, bool add_water, const char* unimod_db_path) {
    if (sequence.find('[') == std::string_view::npos)
        // As above, and also skips resolving/loading unimod_db_path's table
        // (embedded or override) entirely when there's nothing to look up.
        return Iso::FromFASTA(sequence, use_nominal_masses, add_water);

    const UnimodTable& mods = (unimod_db_path == nullptr || unimod_db_path[0] == '\0')
                                   ? embedded_unimod_table()
                                   : unimod_table_for_path(unimod_db_path);
    return FromFASTAWithMods(sequence, mods, use_nominal_masses, add_water);
}

}  // namespace IsoSpec