Skip to content

More verbose error message for nanopore counts #328

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
andrewdavidsmith merged 1 commit into from
Dec 22, 2025
Merged

More verbose error message for nanopore counts #328

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
37 changes: 30 additions & 7 deletions src/analysis/nanopore.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -313,6 +313,11 @@ struct mod_prob_buffer {
bam_parse_basemod(aln.b, m.get());
// ADS: or bam_parse_basemod2(aln.b, m, HTS_MOD_REPORT_UNCHECKED)

int n_types{};
const auto types = bam_mods_recorded(m.get(), &n_types);
if (n_types < 2 || (types[h_idx] != 'h' && types[m_idx] != 'm'))
return false;

int pos{};
int n{};
while ((n = bam_next_basemod(aln.b, m.get(), d, max_mods, &pos)) > 0) {
Expand Down Expand Up @@ -502,7 +507,7 @@ struct mod_prob_stats {

mod_prob_stats() : m{hts_base_mod_state_alloc(), &hts_base_mod_state_free} {};

[[nodiscard]] auto
auto
operator()(const bamxx::bam_rec &aln) {
static constexpr auto h_idx = 0;
static constexpr auto m_idx = 1;
Expand All @@ -514,6 +519,11 @@ struct mod_prob_stats {
bam_parse_basemod(aln.b, m.get());
// ADS: or bam_parse_basemod2(aln.b, m, HTS_MOD_REPORT_UNCHECKED)

int n_types{};
const auto types = bam_mods_recorded(m.get(), &n_types);
if (n_types < 2 || (types[h_idx] != 'h' && types[m_idx] != 'm'))
return;

int pos{};
int n{};
while ((n = bam_next_basemod(aln.b, m.get(), d, max_mods, &pos)) > 0) {
Expand Down Expand Up @@ -907,7 +917,8 @@ struct read_processor {

[[nodiscard]] static auto
valid_modification_types(const std::string &infile,
const std::uint32_t n_reads_to_check) -> bool {
const std::uint32_t n_reads_to_check)
-> std::pair<bool, std::string> {
using mstate = hts_base_mod_state;
std::unique_ptr<mstate, void (*)(mstate *)> m(hts_base_mod_state_alloc(),
&hts_base_mod_state_free);
Expand All @@ -918,6 +929,8 @@ valid_modification_types(const std::string &infile,
if (!hdr)
throw std::runtime_error("failed to read header");

std::string message;

std::uint32_t read_count{};
bool valid_types{true};
bamxx::bam_rec aln;
Expand All @@ -929,9 +942,17 @@ valid_modification_types(const std::string &infile,
// ADS: or bam_parse_basemod2(aln.b, m, HTS_MOD_REPORT_UNCHECKED)
int n_types{};
const auto types = bam_mods_recorded(m.get(), &n_types);
valid_types = n_types >= 2 && types[0] == 'h' && types[1] == 'm';
valid_types = (n_types == 1 && types[0] == 'C') ||
(n_types >= 2 && types[0] == 'h' && types[1] == 'm');
if (!valid_types) {
message = "n_types: " + std::to_string(n_types) + "\n";
for (auto i = 0; i < n_types; ++i) {
message += "type[" + std::to_string(i) +
"]=" + std::to_string(static_cast<char>(types[i])) + "\n";
}
}
}
return valid_types;
return std::make_pair(valid_types, message);
}

[[nodiscard]] static auto
Expand Down Expand Up @@ -1053,9 +1074,11 @@ main_nanocount(int argc, char *argv[]) { // NOLINT(*-avoid-c-arrays)
std::ostringstream cmd;
std::copy(argv, argv + argc, std::ostream_iterator<const char *>(cmd, " "));

if (!valid_modification_types(mapped_reads_file, n_reads_to_check)) {
std::cerr << "modification types are not valid\n"
<< "expected h=0 and m=1\n";
const auto [is_valid, message] =
valid_modification_types(mapped_reads_file, n_reads_to_check);
if (!is_valid) {
std::cerr << "modification types are not valid. Violation:\n"
<< message << '\n';
return EXIT_FAILURE;
}

Expand Down