Skip to content

Commit

Permalink
add different secondary mapping count for seqs less than segLength
Browse files Browse the repository at this point in the history
This helps us to map collections of reads with different sizes, where
some are below the desired segment length for the entire collection. We
can still use these sequences, but we might not want to map them many
places. This makes that configurable, so by default we at least align
them to their best mapping, but we can potentially multi-map too.
  • Loading branch information
ekg committed Oct 30, 2020
1 parent 4bb309a commit 7b5ebf9
Show file tree
Hide file tree
Showing 3 changed files with 14 additions and 3 deletions.
9 changes: 6 additions & 3 deletions src/map/include/computeMap.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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 <MinVec_Type> Q;
Q.seq = &(input->seq)[0u];
Expand Down Expand Up @@ -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
Expand Down
1 change: 1 addition & 0 deletions src/map/include/map_parameters.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<std::string> refSequences; //reference sequence(s)
std::vector<std::string> querySequences; //query sequence(s)
Expand Down
7 changes: 7 additions & 0 deletions src/yeet/include/parse_args.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,7 @@ void parse_args(int argc,
args::ValueFlag<int> kmer_size(parser, "N", "kmer size <= 16 [default: 16]", {'k', "kmer"});
args::ValueFlag<std::string> 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<int> 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<int> 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<char> 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"});
Expand Down Expand Up @@ -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 {
Expand Down

0 comments on commit 7b5ebf9

Please sign in to comment.