Skip to content

Commit

Permalink
Remove Alignment::mapq field
Browse files Browse the repository at this point in the history
Keep mapping quality separate from the alignment itself appears to make more
sense:
- Mapping quality for an alignment is not a property of an alignment per se,
  but is derived from its relationship to other alignments.
- get_alignment returns an Alignment, but leaves mapq uninitialized, that
  is, we must remember to fill it in afterwards, which can easily be
  forgotten.
- Sam::add_pair is passed two mapq values and also two Alignment instances,
  which also contain mapq values. It is unclear for the caller which of the
  values is chosen.
- Adding a mapq parameter to Sam::add makes it mirror Sam::add_pair.
  • Loading branch information
marcelm committed Aug 23, 2023
1 parent 6735747 commit 011d089
Show file tree
Hide file tree
Showing 4 changed files with 8 additions and 14 deletions.
10 changes: 2 additions & 8 deletions src/aln.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -127,8 +127,7 @@ static inline void align_SE(
uint8_t mapq = 60.0 * (best_score - second_best_score) / best_score;

if (max_secondary == 0) {
best_alignment.mapq = mapq;
sam.add(best_alignment, record, read.rc, true, details);
sam.add(best_alignment, record, read.rc, mapq, true, details);
return;
}
// Sort alignments by score, highest first
Expand All @@ -145,12 +144,7 @@ static inline void align_SE(
break;
}
bool is_primary = i == 0;
if (is_primary) {
alignment.mapq = mapq;
} else {
alignment.mapq = 255;
}
sam.add(alignment, record, read.rc, is_primary, details);
sam.add(alignment, record, read.rc, mapq, is_primary, details);
}
}

Expand Down
4 changes: 3 additions & 1 deletion src/sam.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -120,6 +120,7 @@ void Sam::add(
const Alignment& alignment,
const KSeq& record,
const std::string& sequence_rc,
uint8_t mapq,
bool is_primary,
const Details& details
) {
Expand All @@ -131,8 +132,9 @@ void Sam::add(
}
if (!is_primary) {
flags |= SECONDARY;
mapq = 255;
}
add_record(record.name, flags, references.names[alignment.ref_id], alignment.ref_start, alignment.mapq, alignment.cigar, "*", -1, 0, record.seq, sequence_rc, record.qual, alignment.edit_distance, alignment.score, details);
add_record(record.name, flags, references.names[alignment.ref_id], alignment.ref_start, mapq, alignment.cigar, "*", -1, 0, record.seq, sequence_rc, record.qual, alignment.edit_distance, alignment.score, details);
}

// Add one individual record
Expand Down
3 changes: 1 addition & 2 deletions src/sam.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,6 @@ struct Alignment {
int edit_distance;
int global_ed;
int score;
uint8_t mapq;
int length;
bool is_rc;
bool is_unaligned{false};
Expand Down Expand Up @@ -79,7 +78,7 @@ class Sam {
}

/* Add an alignment */
void add(const Alignment& alignment, const klibpp::KSeq& record, const std::string& sequence_rc, bool is_primary, const Details& details);
void add(const Alignment& alignment, const klibpp::KSeq& record, const std::string& sequence_rc, uint8_t mapq, bool is_primary, const Details& details);
void add_pair(const Alignment& alignment1, const Alignment& alignment2, const klibpp::KSeq& record1, const klibpp::KSeq& record2, const std::string& read1_rc, const std::string& read2_rc, uint8_t mapq1, uint8_t mapq2, bool is_proper, bool is_primary, const std::array<Details, 2>& details);
void add_unmapped(const klibpp::KSeq& record, uint16_t flags = UNMAP);
void add_unmapped_pair(const klibpp::KSeq& r1, const klibpp::KSeq& r2);
Expand Down
5 changes: 2 additions & 3 deletions tests/test_sam.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,6 @@ TEST_CASE("Sam::add") {
aln.ref_start = 2;
aln.edit_distance = 3;
aln.score = 9;
aln.mapq = 55;
aln.cigar = Cigar("2S2=1X3=3S");

std::string read_rc = reverse_complement(record.seq);
Expand All @@ -61,15 +60,15 @@ TEST_CASE("Sam::add") {
SUBCASE("Cigar =/X") {
std::string sam_string;
Sam sam(sam_string, references);
sam.add(aln, record, read_rc, is_primary, details);
sam.add(aln, record, read_rc, 55, is_primary, details);
CHECK(sam_string ==
"readname\t16\tcontig1\t3\t55\t2S2=1X3=3S\t*\t0\t0\tACGTT\tBB#>\tNM:i:3\tAS:i:9\n"
);
}
SUBCASE("Cigar M") {
std::string sam_string;
Sam sam(sam_string, references, CigarOps::M);
sam.add(aln, record, read_rc, is_primary, details);
sam.add(aln, record, read_rc, 55, is_primary, details);
CHECK(sam_string ==
"readname\t16\tcontig1\t3\t55\t2S6M3S\t*\t0\t0\tACGTT\tBB#>\tNM:i:3\tAS:i:9\n"
);
Expand Down

0 comments on commit 011d089

Please sign in to comment.