1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
/*
* Copyright (C) 2015-2020 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>.
*/
#pragma once
#include <cstddef>
#include <cstdint>
#include <cstring>
#include <limits>
#include <string_view>
namespace IsoSpec{
// We will work with C H N O S Se tuples */
extern const int aa_isotope_numbers[6];
extern const double aa_elem_masses[19];
extern const double aa_elem_nominal_masses[19];
extern const double aa_elem_probabilities[19];
extern const int aa_symbol_to_elem_counts[256*6];
//! aa_symbol_to_elem_counts again, each row packed into one integer: count ii
//! in bits [10*ii, 10*ii+10). Summing packed rows sums all six counts at once,
//! as long as no field passes 1023 -- see parse_fasta().
struct PackedElemCounts { std::uint64_t row[256]; };
extern const PackedElemCounts aa_symbol_to_packed_elem_counts;
constexpr int AA_PACKED_FIELD_BITS = 10;
//! Residues summed in packed form before unpacking; fasta.cpp checks that the
//! largest count in the table times this still fits a field.
constexpr int AA_PACKED_FLUSH_EVERY = 48;
//! Ceiling the sequence parsers clamp each atom count to, instead of letting a
//! sum overflow `int` (undefined behaviour) -- reachable from a sequence of
//! ~10^8 residues, or from a few modifications with huge override-table
//! deltas. Any count this large is far past the per-element limit Iso enforces
//! (verify_atom_cnt), so a clamped molecule is still rejected there; the
//! headroom up to INT_MAX absorbs the terminal water and modification deltas
//! later added to a clamped count.
constexpr int AA_ATOM_COUNT_CEILING = (std::numeric_limits<int>::max)() / 2;
//! Count elemental composition of an unmodificed sequence of amino acids, resulting in CHNOSSe counts.
/*!
WARNING!!! This function does not add the terminating H and OH groups, resulting in a residue backbone skeleton formula.
If you don't know what that means, you should probably be using parse_fasta_full function.
*/
inline void parse_fasta(std::string_view fasta, int atomCounts[6])
{
// One load and one add per residue, using the packed table, instead of
// six of each: the six-count row add vectorised under clang but not under
// gcc, which made this loop -- most of the cost of reading a protein --
// markedly slower there. The length is taken up front so that the inner
// loop is a counted one, with no per-residue test for the terminator.
// Summed in 64 bits and clamped once at the end (AA_ATOM_COUNT_CEILING):
// 21 hydrogens per residue overflow an int after ~10^8 residues.
// Exactly the bytes of the view are read: a NUL in it is one more
// character with no amino acid, ignored like 'X'.
const unsigned char* p = reinterpret_cast<const unsigned char*>(fasta.data());
const size_t len = fasta.size();
std::uint64_t totals[6] = {0, 0, 0, 0, 0, 0};
for(size_t block = 0; block < len; block += AA_PACKED_FLUSH_EVERY)
{
const size_t block_end = len - block < static_cast<size_t>(AA_PACKED_FLUSH_EVERY) ? len : block + AA_PACKED_FLUSH_EVERY;
std::uint64_t packed = 0;
for(size_t idx = block; idx < block_end; ++idx)
packed += aa_symbol_to_packed_elem_counts.row[p[idx]];
for(int ii = 0; ii < 6; ++ii)
totals[ii] += (packed >> (AA_PACKED_FIELD_BITS * ii)) & ((1u << AA_PACKED_FIELD_BITS) - 1);
}
for(int ii = 0; ii < 6; ++ii)
atomCounts[ii] = totals[ii] < static_cast<std::uint64_t>(AA_ATOM_COUNT_CEILING)
? static_cast<int>(totals[ii]) : AA_ATOM_COUNT_CEILING;
}
//! As above, for a C string: the sequence ends at the first NUL.
inline void parse_fasta(const char* fasta, int atomCounts[6])
{
parse_fasta(std::string_view(fasta), atomCounts);
}
//! Turn an input FASTA aminoacid sequence into atom counts, in CHNOSSe order.
/*!
Unlike parse_fasta, this function includes the H and OH groups at the N- and C- termini of the skeleton, resulting in a formula of a full (inert) molecule.
*/
inline void parse_fasta_full(std::string_view fasta, int atomCounts[6])
{
parse_fasta(fasta, atomCounts);
// Add terminal water (H2O) for either precursor or fragment.
const int H_INDEX = 1; // Indexing: 0=C, 1=H, 2=N, 3=O, 4=S, 5=Se
const int O_INDEX = 3;
atomCounts[H_INDEX] += 2;
atomCounts[O_INDEX] += 1;
}
inline void parse_fasta_full(const char* fasta, int atomCounts[6])
{
parse_fasta_full(std::string_view(fasta), atomCounts);
}
} // namespace IsoSpec