diff --git a/src/map/include/computeMap.hpp b/src/map/include/computeMap.hpp index e91fc447..47ac8b64 100644 --- a/src/map/include/computeMap.hpp +++ b/src/map/include/computeMap.hpp @@ -152,7 +152,7 @@ namespace skch if (param.filterMode == filter::ONETOONE) qmetadata.push_back( ContigInfo{seq_name, len} ); //Is the read too short? - if(len < param.windowSize || len < param.kmerSize || len < param.segLength) + if(len < param.windowSize || len < param.kmerSize) { //#ifdef DEBUG std::cerr << "WARNING, skch::Map::mapQuery, read is not long enough for mapping" << std::endl; @@ -215,7 +215,7 @@ namespace skch //save query sequence name output->qseqName = input->seqName; - if(! param.split) + if(! param.split || input->len <= param.segLength) { QueryMetaData Q; Q.seq = &(input->seq)[0u]; @@ -294,7 +294,10 @@ namespace skch //filter mappings best over query sequence axis if(param.filterMode == filter::MAP || param.filterMode == filter::ONETOONE) { - skch::Filter::query::filterMappings(output->readMappings, param.secondaryToKeep); + skch::Filter::query::filterMappings(output->readMappings, + (input->len < param.segLength ? + param.shortSecondaryToKeep + : param.secondaryToKeep)); } //Make sure mapping boundary don't exceed sequence lengths diff --git a/src/map/include/map_parameters.hpp b/src/map/include/map_parameters.hpp index 034cc5b7..51b5c7b0 100644 --- a/src/map/include/map_parameters.hpp +++ b/src/map/include/map_parameters.hpp @@ -26,6 +26,7 @@ struct Parameters float percentageIdentity; //user defined threshold for good similarity int filterMode; //filtering mode in mashmap int secondaryToKeep; //how many secondary alignments we keep + int shortSecondaryToKeep; //how many secondary alignments we keep for reads < segLength int threads; //execution thread count std::vector refSequences; //reference sequence(s) std::vector querySequences; //query sequence(s) diff --git a/src/yeet/include/parse_args.hpp b/src/yeet/include/parse_args.hpp index 5d695f85..540dbf99 100644 --- a/src/yeet/include/parse_args.hpp +++ b/src/yeet/include/parse_args.hpp @@ -40,6 +40,7 @@ void parse_args(int argc, args::ValueFlag kmer_size(parser, "N", "kmer size <= 16 [default: 16]", {'k', "kmer"}); args::ValueFlag map_filter_mode(parser, "MODE", "filter mode for map step, either 'map', 'one-to-one', or 'none' [default: map]", {'f', "map-filter-mode"}); args::ValueFlag map_secondaries(parser, "N", "number of secondary mappings to retain in 'map' filter mode (total number of mappings is this + 1) [default: 0]", {'n', "n-secondary"}); + args::ValueFlag map_short_secondaries(parser, "N", "number of secondary mappings to retain for sequences shorter than segment length [default: 0]", {'S', "n-short-secondary"}); args::Flag skip_self(parser, "", "skip self mappings when the query and target name is the same (for all-vs-all mode)", {'X', "skip-self"}); args::ValueFlag skip_prefix(parser, "C", "skip mappings when the query and target have the same prefix before the given character C", {'Y', "skip-prefix"}); args::Flag approx_mapping(parser, "approx-map", "skip base-level alignment, producing an approximate mapping in PAF", {'m',"approx-map"}); @@ -218,6 +219,12 @@ void parse_args(int argc, map_parameters.secondaryToKeep = 0; } + if (map_short_secondaries) { + map_parameters.shortSecondaryToKeep = args::get(map_short_secondaries); + } else { + map_parameters.shortSecondaryToKeep = 0; + } + if (skip_self) { map_parameters.skip_self = true; } else {