Merge pull request #4012 from vgteam/surject-score

Report alignment scores in surject

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());
     }

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");

src/subcommand/convert_main.cpp

@@ -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);

test/t/15_vg_surject.t

@@ -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
@@ -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