Skip to content

Commit

Permalink
Merge pull request #4028 from vgteam/require-dist
Browse files Browse the repository at this point in the history 
Require a distance index in vg minimizer
  • Loading branch information
@jltsiren
jltsiren committed Jul 24, 2023
2 parents 2ceff87 + cbc6bbe commit fcb40ed
Show file tree
Hide file tree
Showing 3 changed files with 37 additions and 23 deletions.
25 changes: 18 additions & 7 deletions src/subcommand/minimizer_main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,13 +57,14 @@ int get_default_threads() {
size_t estimate_hash_table_size(const gbwtgraph::GBZ& gbz, bool progress);

void help_minimizer(char** argv) {
std::cerr << "usage: " << argv[0] << " minimizer [options] graph" << std::endl;
std::cerr << "usage: " << argv[0] << " minimizer [options] -d graph.dist -o graph.min graph" << std::endl;
std::cerr << std::endl;
std::cerr << "Builds a (w, k)-minimizer index or a (k, s)-syncmer index of the threads in the GBWT" << std::endl;
std::cerr << "index. The graph can be any HandleGraph, which will be transformed into a GBWTGraph." << std::endl;
std::cerr << "The transformation can be avoided by providing a GBWTGraph or a GBZ graph." << std::endl;
std::cerr << std::endl;
std::cerr << "Required options:" << std::endl;
std::cerr << " -d, --distance-index X annotate the hits with positions in this distance index" << std::endl;
std::cerr << " -o, --output-name X store the index to file X" << std::endl;
std::cerr << std::endl;
std::cerr << "Minimizer options:" << std::endl;
Expand All @@ -81,13 +82,13 @@ void help_minimizer(char** argv) {
std::cerr << " --hash-table N use 2^N-cell hash tables for kmer counting (default: guess)" << std::endl;
std::cerr << std::endl;
std::cerr << "Other options:" << std::endl;
std::cerr << " -d, --distance-index X annotate the hits with positions in this distance index" << std::endl;
std::cerr << " -l, --load-index X load the index from file X and insert the new kmers into it" << std::endl;
std::cerr << " (overrides minimizer / weighted minimizer options)" << std::endl;
std::cerr << " -g, --gbwt-name X use the GBWT index in file X (required with a non-GBZ graph)" << std::endl;
std::cerr << " -p, --progress show progress information" << std::endl;
std::cerr << " -t, --threads N use N threads for index construction (default " << get_default_threads() << ")" << std::endl;
std::cerr << " (using more than " << DEFAULT_MAX_THREADS << " threads rarely helps)" << std::endl;
std::cerr << " --no-dist build the index without distance index annotations (not recommended)" << std::endl;
std::cerr << std::endl;
}

Expand All @@ -105,19 +106,22 @@ int main_minimizer(int argc, char** argv) {
size_t threshold = DEFAULT_THRESHOLD, iterations = DEFAULT_ITERATIONS, hash_table_size = 0;
bool progress = false;
int threads = get_default_threads();
bool require_distance_index = true;

constexpr int OPT_THRESHOLD = 1001;
constexpr int OPT_ITERATIONS = 1002;
constexpr int OPT_FAST_COUNTING = 1003;
constexpr int OPT_SAVE_MEMORY = 1004;
constexpr int OPT_HASH_TABLE = 1005;
constexpr int OPT_NO_DIST = 1100;

int c;
optind = 2; // force optind past command positional argument
while (true) {
static struct option long_options[] =
{
{ "gbwt-name", required_argument, 0, 'g' },
{ "distance-index", required_argument, 0, 'd' },
{ "output-name", required_argument, 0, 'o' },
{ "index-name", required_argument, 0, 'i' }, // deprecated
{ "kmer-length", required_argument, 0, 'k' },
Expand All @@ -131,23 +135,26 @@ int main_minimizer(int argc, char** argv) {
{ "fast-counting", no_argument, 0, OPT_FAST_COUNTING },
{ "save-memory", no_argument, 0, OPT_SAVE_MEMORY },
{ "hash-table", required_argument, 0, OPT_HASH_TABLE },
{ "distance-index", required_argument, 0, 'd' },
{ "load-index", required_argument, 0, 'l' },
{ "gbwt-graph", no_argument, 0, 'G' }, // deprecated
{ "progress", no_argument, 0, 'p' },
{ "threads", required_argument, 0, 't' },
{ "no-dist", no_argument, 0, OPT_NO_DIST },
{ 0, 0, 0, 0 }
};

int option_index = 0;
c = getopt_long(argc, argv, "g:o:i:k:w:bcs:Wd:l:Gpt:h", long_options, &option_index);
c = getopt_long(argc, argv, "g:d:o:i:k:w:bcs:Wl:Gpt:h", long_options, &option_index);
if (c == -1) { break; } // End of options.

switch (c)
{
case 'g':
gbwt_name = optarg;
break;
case 'd':
distance_name = optarg;
break;
case 'o':
output_name = optarg;
break;
Expand Down Expand Up @@ -199,9 +206,6 @@ int main_minimizer(int argc, char** argv) {
}
break;

case 'd':
distance_name = optarg;
break;
case 'l':
load_index = optarg;
break;
Expand All @@ -216,6 +220,9 @@ int main_minimizer(int argc, char** argv) {
threads = std::min(threads, omp_get_max_threads());
threads = std::max(threads, 1);
break;
case OPT_NO_DIST:
require_distance_index = false;
break;

case 'h':
case '?':
Expand All @@ -234,6 +241,10 @@ int main_minimizer(int argc, char** argv) {
return 1;
}
graph_name = argv[optind];
if (require_distance_index && distance_name.empty()) {
std::cerr << "[vg minimizer] error: one of options --distance-index and --no-dist is required" << std::endl;
return 1;
}
if (!load_index.empty() || use_syncmers) {
weighted = false;
}
Expand Down
22 changes: 13 additions & 9 deletions test/t/46_vg_minimizer.t
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ BASH_TAP_ROOT=../deps/bash-tap

PATH=../bin:$PATH # for vg

plan tests 15
plan tests 16


# Indexing a single graph
Expand All @@ -15,33 +15,37 @@ vg gbwt -x x.vg -o x.gbwt -v small/xy2.vcf.gz
vg index -j x.dist x.xg

# Default construction
vg minimizer -o x.mi -g x.gbwt x.xg
vg minimizer --no-dist -o x.mi -g x.gbwt x.xg
is $? 0 "default parameters"

# Construction fails without a distance index
vg minimizer -o x.mi -g x.gbwt x.xg 2> /dev/null
is $? 1 "distance index or --no-dist is required"

# Single-threaded for deterministic results
vg minimizer -t 1 -o x.mi -g x.gbwt x.xg
vg minimizer --no-dist -t 1 -o x.mi -g x.gbwt x.xg
is $? 0 "single-threaded construction"
is $(md5sum x.mi | cut -f 1 -d\ ) 0d75343d78d1e7d9e9fbc3d7d2386ce2 "construction is deterministic"

# Indexing syncmers
vg minimizer -t 1 -o x.mi -c -g x.gbwt x.xg
vg minimizer --no-dist -t 1 -o x.mi -c -g x.gbwt x.xg
is $? 0 "syncmer index"
is $(md5sum x.mi | cut -f 1 -d\ ) 74d836b46799590835c7d61e283df4f0 "construction is deterministic"

# Minimizer parameters
vg minimizer -t 1 -k 7 -w 3 -o x.mi -g x.gbwt -p x.xg
vg minimizer --no-dist -t 1 -k 7 -w 3 -o x.mi -g x.gbwt x.xg
is $? 0 "minimizer parameters"
is $(md5sum x.mi | cut -f 1 -d\ ) 7260e5ea22f063a8dff1f9dd60f92288 "setting -k -w works correctly"

# Construction from GBWTGraph
vg gbwt -x x.xg -g x.gg x.gbwt
vg minimizer -t 1 -g x.gbwt -o x.mi x.gg
vg minimizer --no-dist -t 1 -g x.gbwt -o x.mi x.gg
is $? 0 "construction from GBWTGraph"
is $(md5sum x.mi | cut -f 1 -d\ ) 0d75343d78d1e7d9e9fbc3d7d2386ce2 "construction is deterministic"

# Construction from GBZ
vg gbwt -x x.xg -g x.gbz --gbz-format x.gbwt
vg minimizer -t 1 -o x.mi x.gbz
vg minimizer --no-dist -t 1 -o x.mi x.gbz
is $? 0 "construction from GBZ"
is $(md5sum x.mi | cut -f 1 -d\ ) 0d75343d78d1e7d9e9fbc3d7d2386ce2 "construction is deterministic"

Expand All @@ -62,9 +66,9 @@ vg index -x x.xg -G x.gbwt -v small/xy2.vcf.gz x.vg
vg index -x y.xg -G y.gbwt -v small/xy2.vcf.gz y.vg

# Appending to the index
vg minimizer -t 1 -o x.mi -g x.gbwt x.xg
vg minimizer --no-dist -t 1 -o x.mi -g x.gbwt x.xg
is $? 0 "multiple graphs: first"
vg minimizer -t 1 -l x.mi -o xy.mi -g y.gbwt y.xg
vg minimizer --no-dist -t 1 -l x.mi -o xy.mi -g y.gbwt y.xg
is $? 0 "multiple graphs: second"
is $(md5sum xy.mi | cut -f 1 -d\ ) 1ca39921b15cc3e7d27919a3ec7f47fa "construction is deterministic"

Expand Down
13 changes: 6 additions & 7 deletions test/t/50_vg_giraffe.t
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,9 +13,8 @@ vg gbwt -o x-haps.gbwt -x x.vg -v small/x.vcf.gz
vg gbwt -o x-paths.gbwt -x x.vg --index-paths
vg gbwt -o x-merged.gbwt -m x-haps.gbwt x-paths.gbwt
vg gbwt -o x.gbwt --augment-gbwt -x x.vg x-merged.gbwt
vg snarls --include-trivial x.vg > x.snarls
vg index -j x.dist x.vg
vg minimizer -k 29 -w 11 -g x.gbwt -o x.min x.xg
vg minimizer -k 29 -w 11 -d x.dist -g x.gbwt -o x.min x.xg

vg giraffe -x x.xg -H x.gbwt -m x.min -d x.dist -f reads/small.middle.ref.fq > mapped1.gam
is "${?}" "0" "a read can be mapped with xg + gbwt + min + dist specified without crashing"
Expand Down Expand Up @@ -50,7 +49,7 @@ vg giraffe -Z x.giraffe.gbz -f reads/small.middle.ref.fq --full-l-bonus 0 > mapp
is "$(vg view -aj mapped-nobonus.gam | jq '.score')" "63" "Mapping without a full length bonus produces the correct score"
rm -f mapped-nobonus.gam

vg minimizer -k 29 -b -s 18 -g x.gbwt -o x.sync x.xg
vg minimizer -k 29 -b -s 18 -d x.dist -g x.gbwt -o x.sync x.xg

vg giraffe -x x.xg -H x.gbwt -m x.sync -d x.dist -f reads/small.middle.ref.fq > mapped.sync.gam
is "${?}" "0" "a read can be mapped with syncmer indexes without crashing"
Expand All @@ -59,7 +58,7 @@ chmod 400 x.dist
vg giraffe -x x.xg -H x.gbwt -m x.min -d x.dist -f reads/small.middle.ref.fq >/dev/null
is "${?}" "0" "a read can be mapped when the distance index is not writable"

rm -f x.vg x.xg x.gbwt x.snarls x.min x.sync x.dist x.gg
rm -f x.vg x.xg x.gbwt x.min x.sync x.dist x.gg
rm -f x.giraffe.gbz


Expand Down Expand Up @@ -114,10 +113,10 @@ is "$(cat surjected.sam | grep -v '^@' | cut -f 1 | sort | uniq | wc -l)" "2" "s
is "$(cat surjected.sam | grep -v '^@' | cut -f 7)" "$(printf '*\n*')" "surjection of unpaired reads to SAM produces absent partner contigs"
is "$(cat surjected.sam | grep -v '^@' | sort -k4 | cut -f 2)" "$(printf '0\n16')" "surjection of unpaired reads to SAM produces correct flags"

rm -f x.vg x.gbwt x.xg x.snarls x.min x.dist x.gg x.fa x.fa.fai x.vcf.gz x.vcf.gz.tbi single.gam paired.gam surjected.sam
rm -f x.vg x.gbwt x.xg x.min x.dist x.gg x.fa x.fa.fai x.vcf.gz x.vcf.gz.tbi single.gam paired.gam surjected.sam
rm -f x.giraffe.gbz

rm -f xy.vg xy.gbwt xy.xg xy.snarls xy.min xy.dist xy.gg xy.fa xy.fa.fai xy.vcf.gz xy.vcf.gz.tbi
rm -f xy.vg xy.gbwt xy.xg xy.min xy.dist xy.gg xy.fa xy.fa.fai xy.vcf.gz xy.vcf.gz.tbi
cp small/xy.fa .
cp small/xy.vcf.gz .
cp small/xy.vcf.gz.tbi .
Expand Down Expand Up @@ -171,7 +170,7 @@ is "$(cat xy.json | grep "correct-minimizer-coverage" | wc -l)" "2000" "paired r
is "$(cat xy.json | jq '.annotation["correct-minimizer-coverage"] | select(. > 0)' | wc -l)" "2000" "paired reads all have nonzero correct minimizer coverage"

rm -f x.vg xy.sam xy.json
rm -f xy.vg xy.gbwt xy.xg xy.snarls xy.min xy.dist xy.gg xy.fa xy.fa.fai xy.vcf.gz xy.vcf.gz.tbi
rm -f xy.vg xy.gbwt xy.xg xy.min xy.dist xy.gg xy.fa xy.fa.fai xy.vcf.gz xy.vcf.gz.tbi

vg giraffe -Z xy.giraffe.gbz -G x.gam -o BAM >xy.bam
is $? "0" "vg giraffe can map to BAM using a GBZ alone"
Expand Down

1 comment on commit fcb40ed

@adamnovak
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

vg CI tests complete for merge to master. View the full report here.

15 tests passed, 0 tests failed and 0 tests skipped in 36822 seconds

Please sign in to comment.