Skip to content

Fixing A-rich vs. T-rich tags for reads that map to the reverse strand #107

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 2 commits into from
May 19, 2026
+36 −37
Merged

Fixing A-rich vs. T-rich tags for reads that map to the reverse strand #107

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ebe89bc
src/abismal.cpp: toggling the conversion T-rich vs. A-rich if the rea…
andrewdavidsmith May 19, 2026
ff88e65
data/md5sum.txt: updating test checksums
andrewdavidsmith May 19, 2026
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
20 changes: 10 additions & 10 deletions data/md5sum.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,16 +1,16 @@
e95a8739a378bc5628c76ffbe8293682 tests/reads_1.fq
5487eb32b5492ddba3a95a8461e2ff4c tests/reads_pbat_pe_1.fq
aff19aeae4184c38cf364b3a94527098 tests/reads_pbat_pe_2.fq
46e60a92266ed0f913f5efb9b7c0473c tests/reads.sam
514bf940e7d5f44e57291a75e1fc6629 tests/reads_pe_1.fq
5487eb32b5492ddba3a95a8461e2ff4c tests/reads_pbat_pe_1.fq
57fd6ff3b49f4656c9506688040b7ce6 tests/tRex1.idx
64e17fc1e424f9cb21879d38ddf745e5 tests/reads_pe_2.fq
e94e71292fccd255bad3f3694efe7a4b tests/reads_rpbat_pe_1.fq
8bfa00acd0639dc66acd5aa8ac0369d5 tests/reads_rpbat_pe_2.fq
e8a5fe0aec564a972ca7a3a36c0022f7 tests/reads_pbat_pe.sam
14e935a5b8a20c8a7422834c3327d1c7 tests/reads_pe.sam
98d6ad9d563e3318756bd5fc3fcbc263 tests/reads_rpbat_pe.sam
8126d46074213ad3674181f4ea4f8bd1 tests/reads.sam
981dd110d1675c77533a485204fb13cc tests/reads.mstats
bd54ba1e720039cb4b662faa4f0e04b3 tests/reads_rpbat_pe.mstats
9b78ffb855d8a94604973118f2e19457 tests/reads_pbat_pe.sam
a6a3a12238013d8081c3e06dd43d538c tests/reads_pe.sam
a6e138f4b89ef2dc0b90004baadddb9b tests/reads_pbat_pe.mstats
aff19aeae4184c38cf364b3a94527098 tests/reads_pbat_pe_2.fq
bd54ba1e720039cb4b662faa4f0e04b3 tests/reads_rpbat_pe.mstats
d7a9e3f6a5ed7373a827f07aaf32803a tests/reads_pe.mstats
57fd6ff3b49f4656c9506688040b7ce6 tests/tRex1.idx
e94e71292fccd255bad3f3694efe7a4b tests/reads_rpbat_pe_1.fq
e95a8739a378bc5628c76ffbe8293682 tests/reads_1.fq
edf2c97162f7910949520ed60d82f205 tests/reads_rpbat_pe.sam
53 changes: 26 additions & 27 deletions src/abismal.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -484,6 +484,7 @@ format_se(const bool allow_ambig, const se_element &res, const ChromLookup &cl,
// ADS: 'read' should not be used after a call to 'format_se'
static constexpr auto mapq_max_val = 255u;
static constexpr auto aux_len = 16u;
static constexpr auto cv = std::array<std::uint8_t, 2>{'T', 'A'};

const bool ambig = res.ambig();
const bool valid = !res.empty();
Expand All @@ -498,9 +499,11 @@ format_se(const bool allow_ambig, const se_element &res, const ChromLookup &cl,

// ADS: we might be doing format_se for an end in paried reads
std::uint16_t flag{};
auto is_a_rich = res.elem_is_a_rich();
if (res.rc()) {
flag |= BAM_FREVERSE;
revcomp_inplace(r.read);
is_a_rich = !is_a_rich;
}

if (allow_ambig && ambig)
Expand All @@ -511,21 +514,21 @@ format_se(const bool allow_ambig, const se_element &res, const ChromLookup &cl,
r.rec.b = bam_init1();
// clang-format off
int ret = bam_set1(r.rec.b,
std::size(r.name), // size_t l_qname,
std::data(r.name), // const char *qname,
flag, // uint16_t flag,
chrom_idx - 1, // int32_t tid (-1 for padding)
ref_s, // hts_pos_t pos,
mapq_max_val, // std::uint8_t mapq,
std::size(r.cig), // size_t n_cigar,
std::data(r.cig), // const uint32_t *cigar,
-1, // int32_t mtid,
-1, // hts_pos_t mpos,
0, // hts_pos_t isize,
std::size(r.read), // size_t l_seq,
std::data(r.read), // const char *seq,
nullptr, // const char *qual,
aux_len); // size_t l_aux);
std::size(r.name), // size_t l_qname,
std::data(r.name), // const char *qname,
flag, // uint16_t flag,
chrom_idx - 1, // int32_t tid (-1 for padding)
ref_s, // hts_pos_t pos,
mapq_max_val, // std::uint8_t mapq,
std::size(r.cig), // size_t n_cigar,
std::data(r.cig), // const uint32_t *cigar,
-1, // int32_t mtid,
-1, // hts_pos_t mpos,
0, // hts_pos_t isize,
std::size(r.read), // size_t l_seq,
std::data(r.read), // const char *seq,
nullptr, // const char *qual,
aux_len); // size_t l_aux);
// clang-format on
if (ret < 0)
throw std::runtime_error("failed to format bam");
Expand All @@ -534,10 +537,7 @@ format_se(const bool allow_ambig, const se_element &res, const ChromLookup &cl,
if (ret < 0)
throw std::runtime_error("bam_aux_update_int");

ret = bam_aux_append(
r.rec.b, "CV", 'A', 1,
// NOLINTNEXTLINE(*-reinterpret-cast)
reinterpret_cast<const std::uint8_t *>(res.elem_is_a_rich() ? "A" : "T"));
ret = bam_aux_append(r.rec.b, "CV", 'A', 1, std::data(cv) + is_a_rich);
if (ret < 0)
throw std::runtime_error("bam_aux_append");

Expand Down Expand Up @@ -652,10 +652,7 @@ format_pe(
read_holder &r1, read_holder &r2) -> map_type {
static constexpr auto mapq_max_val = 255u;
static constexpr auto aux_len = 16u;
static const std::array<std::uint8_t, 2> cv = {
'T',
'A',
};
static constexpr auto cv = std::array<std::uint8_t, 2>{'T', 'A'};

if (p.empty())
return map_unmapped;
Expand Down Expand Up @@ -689,15 +686,19 @@ format_pe(

flag2 |= BAM_FPAIRED | BAM_FPROPER_PAIR;

auto r1_is_a_rich = p.r1.elem_is_a_rich();
if (p.r1.rc()) { // ADS: is p.r1.rc() always !p.r2.rc()?
flag1 |= BAM_FREVERSE;
flag2 |= BAM_FMREVERSE;
revcomp_inplace(r1.read);
r1_is_a_rich = !r1_is_a_rich;
}
auto r2_is_a_rich = p.r2.elem_is_a_rich();
if (p.r2.rc()) {
flag2 |= BAM_FREVERSE;
flag1 |= BAM_FMREVERSE;
revcomp_inplace(r2.read);
r2_is_a_rich = !r2_is_a_rich;
}
if (allow_ambig && ambig) {
// ADS: mark ambig for both the same way?
Expand Down Expand Up @@ -733,8 +734,7 @@ format_pe(
if (ret < 0)
throw std::runtime_error("error adding aux field");

ret = bam_aux_append(r1.rec.b, "CV", 'A', 1,
std::data(cv) + p.r1.elem_is_a_rich());
ret = bam_aux_append(r1.rec.b, "CV", 'A', 1, std::data(cv) + r1_is_a_rich);
if (ret < 0)
throw std::runtime_error("error adding aux field");

Expand Down Expand Up @@ -764,8 +764,7 @@ format_pe(
if (ret < 0)
throw std::runtime_error("error adding aux field");

ret = bam_aux_append(r2.rec.b, "CV", 'A', 1,
std::data(cv) + p.r2.elem_is_a_rich());
ret = bam_aux_append(r2.rec.b, "CV", 'A', 1, std::data(cv) + r2_is_a_rich);
if (ret < 0)
throw std::runtime_error("error adding aux field");

Expand Down
Loading