Skip to content
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

Jump to bottom

Report alignment scores in surject #4012

Merged
merged 2 commits into from
Jul 7, 2023
Merged

Report alignment scores in surject #4012

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
f8cf5f1
report alignment scores from surject
jeizenga Jul 6, 2023
21a972b
add integration test for scores
jeizenga Jul 6, 2023
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
6 changes: 6 additions & 0 deletions src/alignment.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -805,6 +805,12 @@ bam1_t* alignment_to_bam_internal(bam_hdr_t* header,
}
}

if (~core.flag & BAM_FUNMAP) {
// we've decided that it is aligned
int32_t score = alignment.score();
bam_aux_append(bam, "AS", 'i', sizeof(int32_t), (uint8_t*) &score);
}

if (!alignment.read_group().empty()) {
bam_aux_append(bam, "RG", 'Z', alignment.read_group().size() + 1, (uint8_t*) alignment.read_group().c_str());
}
Expand Down
5 changes: 4 additions & 1 deletion src/multipath_alignment_emitter.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -337,10 +337,13 @@ void MultipathAlignmentEmitter::create_alignment_shim(const string& name, const
mapped = true;
}
}
// hacky way to inform the conversion code that the read is mapped
if (mapped) {
// hacky way to inform the conversion code that the read is mapped
shim.mutable_path()->add_mapping();
// and we'll also communicate the alignment score
shim.set_score(optimal_alignment_score(mp_aln, true));
}

// this tag comes from surject and is used in both
if (mp_aln.has_annotation("all_scores")) {
auto anno = mp_aln.get_annotation("all_scores");
Expand Down
2 changes: 1 addition & 1 deletion src/subcommand/convert_main.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -187,7 +187,7 @@ int main_convert(int argc, char** argv) {
{
num_threads = parse<int>(optarg);
if (num_threads <= 0) {
cerr << "error:[vg mpmap] Thread count (-t) set to " << num_threads << ", must set to a positive integer." << endl;
cerr << "error:[vg convert] Thread count (-t) set to " << num_threads << ", must set to a positive integer." << endl;
exit(1);
}
omp_set_num_threads(num_threads);
Expand Down
6 changes: 4 additions & 2 deletions test/t/15_vg_surject.t
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ BASH_TAP_ROOT=../deps/bash-tap
PATH=../bin:$PATH # for vg


plan tests 46
plan tests 47

vg construct -r small/x.fa >j.vg
vg index -x j.xg j.vg
Expand Down Expand Up @@ -36,7 +36,9 @@ is $(vg surject -p x -x x.xg -t 1 -s j.gam | grep -v "@" | cut -f3 | grep x | wc

is $(vg surject -x x.xg -t 1 -s j.gam | grep -v "@" | cut -f3 | grep x | wc -l) \
100 "vg surject doesn't need to be told which path to use"


is $(vg surject -x x.xg -t 1 -s x.gam | grep AS | wc -l) 100 "vg surject reports alignment scores"

vg paths -X -x x.vg | vg view -aj - | jq '.name = "sample#0#x#0"' | vg view -JGa - > paths.gam
vg paths -X -x x.vg | vg view -aj - | jq '.name = "ref#0#x[55]"' | vg view -JGa - >> paths.gam
vg augment x.vg -i paths.gam > x.aug.vg
Expand Down