From f8cf5f112a3d202013746f9feac966ee5dcdcfe1 Mon Sep 17 00:00:00 2001 From: Jordan Eizenga Date: Thu, 6 Jul 2023 11:48:54 -0700 Subject: [PATCH] report alignment scores from surject --- src/alignment.cpp | 6 ++++++ src/multipath_alignment_emitter.cpp | 5 ++++- src/subcommand/convert_main.cpp | 2 +- 3 files changed, 11 insertions(+), 2 deletions(-) diff --git a/src/alignment.cpp b/src/alignment.cpp index 633deea2888..665a84b03f8 100644 --- a/src/alignment.cpp +++ b/src/alignment.cpp @@ -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()); } diff --git a/src/multipath_alignment_emitter.cpp b/src/multipath_alignment_emitter.cpp index 1d7ef25a610..705910b1dff 100644 --- a/src/multipath_alignment_emitter.cpp +++ b/src/multipath_alignment_emitter.cpp @@ -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"); diff --git a/src/subcommand/convert_main.cpp b/src/subcommand/convert_main.cpp index 3d6bb0c9b7f..710bdb4458a 100644 --- a/src/subcommand/convert_main.cpp +++ b/src/subcommand/convert_main.cpp @@ -187,7 +187,7 @@ int main_convert(int argc, char** argv) { { num_threads = parse(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);