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.

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

Make vg inject use AlignmentEmitter and support GAF #4319

Merged
merged 1 commit into from
Jul 1, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
36 changes: 27 additions & 9 deletions src/subcommand/inject_main.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
#include "../alignment.hpp"
#include "../vg.hpp"
#include "../xg.hpp"
#include "../hts_alignment_emitter.hpp"
#include <vg/io/stream.hpp>
#include <vg/io/vpkg.hpp>
#include <bdsg/overlays/overlay_helper.hpp>
Expand All @@ -23,10 +24,11 @@ using namespace vg;
using namespace vg::subcommand;

void help_inject(char** argv) {
cerr << "usage: " << argv[0] << " inject [options] input.[bam|sam|cram] >output.gam" << endl
cerr << "usage: " << argv[0] << " inject -x graph.xg [options] input.[bam|sam|cram] >output.gam" << endl
<< endl
<< "options:" << endl
<< " -x, --xg-name FILE use this graph or xg index (required, non-XG formats also accepted)" << endl
<< " -o, --output-format NAME output the alignments in NAME format (gam / gaf / json) [gam]" << endl
<< " -t, --threads N number of threads to use" << endl;
}

Expand All @@ -37,6 +39,8 @@ int main_inject(int argc, char** argv) {
}

string xg_name;
string output_format = "GAM";
std::set<std::string> output_formats = { "GAM", "GAF", "JSON" };
int threads = get_thread_count();

int c;
Expand All @@ -46,12 +50,13 @@ int main_inject(int argc, char** argv) {
{
{"help", no_argument, 0, 'h'},
{"xg-name", required_argument, 0, 'x'},
{"output-format", required_argument, 0, 'o'},
{"threads", required_argument, 0, 't'},
{0, 0, 0, 0}
};

int option_index = 0;
c = getopt_long (argc, argv, "hx:t:",
c = getopt_long (argc, argv, "hx:o:t:",
long_options, &option_index);

// Detect the end of the options.
Expand All @@ -63,6 +68,19 @@ int main_inject(int argc, char** argv) {
case 'x':
xg_name = optarg;
break;

case 'o':
{
output_format = optarg;
for (char& c : output_format) {
c = std::toupper(c);
}
if (output_formats.find(output_format) == output_formats.end()) {
std::cerr << "error: [vg inject] Invalid output format: " << optarg << std::endl;
std::exit(1);
}
}
break;

case 't':
threads = parse<int>(optarg);
Expand Down Expand Up @@ -90,14 +108,14 @@ int main_inject(int argc, char** argv) {
}
unique_ptr<PathHandleGraph> path_handle_graph = vg::io::VPKG::load_one<PathHandleGraph>(xg_name);
bdsg::PathPositionOverlayHelper overlay_helper;
PathPositionHandleGraph* xgidx = overlay_helper.apply(path_handle_graph.get());
PathPositionHandleGraph* xgidx = overlay_helper.apply(path_handle_graph.get());

vg::io::ProtobufEmitter<Alignment> buf(cout);
function<void(Alignment&)> lambda = [&buf](Alignment& aln) {
#pragma omp critical (buf)
{
buf.write(std::move(aln));
}
// We don't do HTS output formats but we do need an empty paths collection to make an alignment emitter
vector<tuple<path_handle_t, size_t, size_t>> paths;
unique_ptr<vg::io::AlignmentEmitter> alignment_emitter = get_alignment_emitter("-", output_format, paths, threads, xgidx);

function<void(Alignment&)> lambda = [&](Alignment& aln) {
alignment_emitter->emit_mapped_single({std::move(aln)});
};
if (threads > 1) {
hts_for_each_parallel(file_name, lambda, xgidx);
Expand Down
4 changes: 3 additions & 1 deletion test/t/39_vg_inject.t
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 12
plan tests 13

vg construct -r small/x.fa > j.vg
vg index -x j.xg j.vg
Expand Down Expand Up @@ -52,5 +52,7 @@ cat <(samtools view -H small/x.bam) <(printf "name\t4\t*\t0\t0\t*\t*\t0\t0\tACGT
is "$(vg inject -x x.xg unmapped.sam | vg view -aj - | grep "path" | wc -l)" 0 "vg inject does not make an alignment for an umapped read"
is "$(echo $?)" 0 "vg inject does not crash on an unmapped read"

is $(vg inject -x x.xg small/x.bam -o GAF | wc -l) \
1000 "vg inject supports GAF output"

rm j.vg j.xg x.vg x.gcsa x.gcsa.lcp x.xg unmapped.sam
Loading