Skip to content

Give vg surject opinions about paired-end input #5021

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.

Sign up for GitHub

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

Merged
adamnovak merged 4 commits into from
Aug 26, 2026
+112 −37
Merged

Give vg surject opinions about paired-end input #5021

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
fb5ea18
Give vg surject opinions about paired-end input
adamnovak Aug 25, 2026
87a5bf1
Revise message text
adamnovak Aug 25, 2026
c098382
Remove no-longer-generated annotation from test cases
adamnovak Aug 25, 2026
606ceac
make helptext allowable length
faithokamoto Aug 25, 2026
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
5 changes: 3 additions & 2 deletions src/minimizer_mapper.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2803,14 +2803,15 @@ pair<vector<Alignment>, vector<Alignment>> MinimizerMapper::map_paired(Alignment
}
}

//Annotate top pair with its fragment distance, properly-paired-ness, fragment length distrubution, and secondary scores
// Annotate top pair with its fragment distance, properly-paired-ness, fragment length distrubution, and secondary scores
bool properly_paired = distances.front() == std::numeric_limits<int64_t>::max() ? false :
(std::abs(distances.front()-fragment_length_distr.mean()) <= 6.0*fragment_length_distr.std_dev()) ;
string distribution = "-I " + to_string(fragment_length_distr.mean()) + " -D " + to_string(fragment_length_distr.std_dev());
for (auto r : {0, 1}) {
set_annotation(mappings[r].front(), "fragment_length", distance_to_annotation(distances.front()));
set_annotation(mappings[r].front(), "proper_pair", properly_paired);
set_annotation(mappings[r].front(),"fragment_length_distribution", distribution);
// GAM actually has a fragment_length_distribution field.
mappings[r].front().set_fragment_length_distribution(distribution);
set_annotation(mappings[r].front(),"secondary_scores", scores);
}
}
Expand Down
112 changes: 88 additions & 24 deletions src/subcommand/surject_main.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,15 +46,19 @@ void help_surject(char** argv) {
<< " -F, --into-paths FILE surject into path names listed in" << endl
<< " HTSlib sequence dictionary or path list FILE" << endl
<< " -n, --into-ref NAME surject into this reference assembly" << endl
<< " -i, --interleaved GAM is interleaved paired-ended, so pair reads" << endl
<< " when outputting HTS formats" << endl
<< " -M, --multimap include secondary alignments to all" << endl
<< " overlapping paths instead of just primary" << endl
<< " -G, --gaf-input input file is GAF instead of GAM" << endl
<< " -m, --gamp-input input file is GAMP instead of GAM" << endl
<< " -c, --cram-output write CRAM instead of GAM to stdout" << endl
<< " -b, --bam-output write BAM instead of GAM to stdout" << endl
<< " -s, --sam-output write SAM instead of GAM to stdout" << endl
<< " -i, --interleaved input is interleaved paired-ended, so pair reads" << endl
<< " when outputting HTS formats (SAM/BAM/CRAM)" << endl
<< " -U, --force-unpaired surject reads as unpaired even if they appear paired" << endl
<< " -f, --max-frag-len N reads with fragment lengths greater than N won't be" << endl
<< " marked properly paired in HTS formats." << endl
<< " 0 for unlimited. (default: unlimited)" << endl
<< " -u, --supplementary divide into supplementary alignments as necessary" << endl
<< " -B, --left-align attempt to left-align indels" << endl
<< " -l, --subpath-local let the multipath mapping surjection produce local" << endl
Expand All @@ -78,8 +82,6 @@ void help_surject(char** argv) {
<< " of the 10x-scaled scoring parameters" << endl
<< " -N, --sample NAME set this sample name for all reads" << endl
<< " -R, --read-group NAME set this read group for all reads" << endl
<< " -f, --max-frag-len N reads with fragment lengths greater than N won't be" << endl
<< " marked properly paired in SAM/BAM/CRAM" << endl
<< " -L, --list-all-paths annotate SAM records with a list of all attempted" << endl
<< " re-alignments to paths in SS tag" << endl
<< " -H, --graph-aln annotate SAM records with cs-style difference string" << endl
Expand Down Expand Up @@ -136,10 +138,48 @@ static void ensure_alignment_is_for_graph(const Logger& logger, const MultipathA
}
}

/// Returns true if the given alignment seems like it was aligned paired-end,
/// and false otherwise. Should not detect paired-end reads that were aligned
/// single-ended.
static bool smells_paired(const Alignment& aln) {
// The vg way to find this out is the fragment_prev/fragment_next fields.
// vg GAF has tags for them
return aln.has_fragment_prev() ||
aln.has_fragment_next() ||
// If an alignment just came in from BAM somehow, this would be set.
aln.read_paired() ||
// vg GAF will record a tag that carries this.
// Maybe some other GAF writer supports this but not next/prev?
has_annotation(aln, "proper_pair") ||
// Supplementary alignments might have this instead.
has_annotation(aln, "mate_info");
}

/// Returns true if the given multipath alignment seems like it was aligned
/// paired-end, and false otherwise. Should not detect paired-end reads that
/// were aligned single-ended.
static bool smells_paired(const MultipathAlignment& mpaln) {
return has_annotation(mpaln, "fragment_length_distribution") ||
has_annotation(mpaln, "proper_pair") ||
has_annotation(mpaln, "mate_info");
}

/// Helper for producing an error message about a read that smells paired when
/// running in single-ended mode.
static void smells_paired_error(const Logger& logger, const string& name) {
#pragma omp critical (cerr)
logger.error() << "Read " << name << " seems to have been aligned as a "
<< "paired-end read. Provide either the -i/--interleaved "
<< "option (and, ideally, a value for -f/--max-frag-len) "
<< "to surject as paired-end, or -U/--force-unpaired "
<< "to surject as single-ended anyway." << endl;
}

/// This is a helper for the cases where we find two alignments
/// that are adjacent in the input but not paired.
static void adjacent_but_not_paired_error(const Logger& logger, const string& name1, const string& name2) {


#pragma omp critical (cerr)
logger.error() << "alignments " << name1 << " and "
<< name2 << " are adjacent but not paired" << endl;
}
Expand All @@ -163,9 +203,10 @@ int main_surject(int argc, char** argv) {
string input_format = "GAM";
bool spliced = false;
bool interleaved = false;
bool force_unpaired = false;
std::optional<int32_t> max_frag_len;
string sample_name;
string read_group;
int32_t max_frag_len = 0;
int compress_level = 9;
int64_t min_splice_length = 20;
size_t watchdog_timeout = 10;
Expand Down Expand Up @@ -204,13 +245,15 @@ int main_surject(int argc, char** argv) {
{"subpath-local", no_argument, 0, 'l'},
{"max-tail-len", required_argument, 0, 'T'},
{"max-graph-scale", required_argument, 0, 'g'},
{"interleaved", no_argument, 0, 'i'},
{"multimap", no_argument, 0, 'M'},
{"gaf-input", no_argument, 0, 'G'},
{"gamp-input", no_argument, 0, 'm'},
{"cram-output", no_argument, 0, 'c'},
{"bam-output", no_argument, 0, 'b'},
{"sam-output", no_argument, 0, 's'},
{"interleaved", no_argument, 0, 'i'},
{"force-unpaired", no_argument, 0, 'U'},
{"max-frag-len", required_argument, 0, 'f'},
{"supplementary", no_argument, 0, 'u'},
{"off-ref-position", no_argument, 0, OPT_OFF_REF_POS},
{"left-align", no_argument, 0, 'B'},
Expand All @@ -224,7 +267,6 @@ int main_surject(int argc, char** argv) {
{"extra-gap-cost", required_argument, 0, 'E'},
{"sample", required_argument, 0, 'N'},
{"read-group", required_argument, 0, 'R'},
{"max-frag-len", required_argument, 0, 'f'},
{"list-all-paths", no_argument, 0, 'L'},
{"graph-aln", no_argument, 0, 'H'},
{"compression", required_argument, 0, 'C'},
Expand All @@ -235,7 +277,7 @@ int main_surject(int argc, char** argv) {
};

int option_index = 0;
c = getopt_long (argc, argv, "h?x:p:F:n:lT:g:iGmcbsuBN:R:f:C:t:D:SPI:a:AE:LHMVw:r",
c = getopt_long (argc, argv, "h?x:p:F:n:lT:g:GmcbsiUf:uBN:R:C:t:D:SPI:a:AE:LHMVw:r",
long_options, &option_index);

// Detect the end of the options.
Expand Down Expand Up @@ -273,10 +315,6 @@ int main_surject(int argc, char** argv) {
max_graph_scale.reset(new double(parse<double>(optarg)));
break;

case 'i':
interleaved = true;
break;

case 'M':
multimap = true;
break;
Expand All @@ -302,6 +340,18 @@ int main_surject(int argc, char** argv) {
output_format = "SAM";
break;

case 'i':
interleaved = true;
break;

case 'U':
force_unpaired = true;
break;

case 'f':
max_frag_len = parse<int32_t>(optarg);
break;

case 'u':
report_supplementary = true;
break;
Expand Down Expand Up @@ -350,10 +400,6 @@ int main_surject(int argc, char** argv) {
read_group = optarg;
break;

case 'f':
max_frag_len = parse<int32_t>(optarg);
break;

case 'C':
compress_level = parse<int>(optarg);
break;
Expand Down Expand Up @@ -407,6 +453,21 @@ int main_surject(int argc, char** argv) {
prune_anchors = (read_length == "long");
}

// Validate configuration
if (!interleaved && max_frag_len.has_value()) {
logger.error() << "-f/--max-frag-len can only be used with paired-end reads, "
<< "but -i/--interleaved was not provided." << endl;
}
if (interleaved && !max_frag_len.has_value()) {
// TODO: Once we get fragment distribution learning, this will be less of a problem.
logger.warn() << "Running in paired-end mode without -f/--max-frag-len. "
<< "Reads will be assumed to be properly paired at any distance!" << endl;
}
if (interleaved && force_unpaired) {
logger.error() << "Cannot use -U/--force-unpaired and -i/--interleaved at the same time. "
<< "To surject interleaved reads as unpaired, use just -U/--force-unpaired." << endl;
}

string file_name = get_input_file_name(optind, argc, argv);

if (have_input_file(optind, argc, argv)) {
Expand Down Expand Up @@ -562,7 +623,6 @@ int main_surject(int argc, char** argv) {
!src2.has_fragment_prev() ||
src2.fragment_prev().name() != src1.name()) {

#pragma omp critical (cerr)
adjacent_but_not_paired_error(logger, src1.name(), src2.name());

}
Expand All @@ -572,13 +632,11 @@ int main_surject(int argc, char** argv) {
!src1.has_fragment_prev() ||
src1.fragment_prev().name() != src2.name()) {

#pragma omp critical (cerr)
adjacent_but_not_paired_error(logger, src1.name(), src2.name());

}
} else {
// Alignments aren't paired up at all
#pragma omp critical (cerr)
adjacent_but_not_paired_error(logger, src1.name(), src2.name());
}

Expand Down Expand Up @@ -616,7 +674,7 @@ int main_surject(int argc, char** argv) {
auto it = strand_idx2.find(make_pair(pos.name(), !pos.is_reverse()));
if (!is_supplementary(surjected1[i]) && it != strand_idx2.end()) {
// the alignments are paired on this strand
alignment_emitter->emit_pair(std::move(surjected1[i]), std::move(surjected2[it->second]), max_frag_len);
alignment_emitter->emit_pair(std::move(surjected1[i]), std::move(surjected2[it->second]), max_frag_len.value_or(0));
}
else {
// supplementary or unpaired
Expand Down Expand Up @@ -695,6 +753,10 @@ int main_surject(int argc, char** argv) {
ensure_alignment_is_for_graph(logger, src, *xgidx);
}

if (!force_unpaired && smells_paired(src)) {
smells_paired_error(logger, src.name());
}

// Preprocess read to set metadata before surjection
set_metadata(src);

Expand Down Expand Up @@ -750,13 +812,11 @@ int main_surject(int argc, char** argv) {
// TODO: Integrate into for_each_interleaved_pair_parallel when running on Alignments.
if (src1.paired_read_name() != src2.name() || src2.paired_read_name() != src1.name()) {

#pragma omp critical (cerr)
adjacent_but_not_paired_error(logger, src1.name(), src2.name());

}
else if (src1.paired_read_name().empty() || src2.paired_read_name().empty()) {
// Alignments aren't paired up at all
#pragma omp critical (cerr)
adjacent_but_not_paired_error(logger, src1.name(), src2.name());
}

Expand Down Expand Up @@ -861,7 +921,7 @@ int main_surject(int argc, char** argv) {
}

// write to output
vector<int64_t> tlen_limits(surjected.size(), max_frag_len);
vector<int64_t> tlen_limits(surjected.size(), max_frag_len.value_or(0));
mp_alignment_emitter.emit_pairs(src1.name(), src2.name(), std::move(surjected), &positions, &tlen_limits);
mp_alignment_emitter.emit_paired_independent(src1.name(), src2.name(),
std::move(surjected_unpaired1), std::move(surjected_unpaired2),
Expand Down Expand Up @@ -891,6 +951,10 @@ int main_surject(int argc, char** argv) {
ensure_alignment_is_for_graph(logger, src, *xgidx);
}

if (!force_unpaired && smells_paired(src)) {
smells_paired_error(logger, src.name());
}

multipath_alignment_t mp_src;
from_proto_multipath_alignment(src, mp_src);

Expand Down
2 changes: 1 addition & 1 deletion src/surjector.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ using namespace std;

Surjector::Surjector(const PathPositionHandleGraph* graph) : graph(graph), choose_band_padding(algorithms::pad_band_constant(1)) {
if (!graph) {
cerr << "error:[Surjector] Failed to provide an graph to the Surjector" << endl;
cerr << "error:[Surjector] Failed to provide a graph to the Surjector" << endl;
exit(1);
}

Expand Down
26 changes: 18 additions & 8 deletions test/t/15_vg_surject.t
View file
Open in desktop
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 78
plan tests 81

vg construct -r small/x.fa >j.vg
vg index -x j.xg j.vg
Expand Down Expand Up @@ -192,11 +192,11 @@ is "$(vg surject -x j.vg -b --graph-aln r.gam | samtools view | grep 'GR:Z:' | w

rm -f h.vg h.gcsa r.gam r.sam x.sub.fa j.vg j.gcsa j.gcsa.lcp j.sub.vg j.sub.gcsa j.sub.gcsa.lcp r.sub.gam r.sub.sam r.sub.sam path_info.tsv

vg surject -s -x surject/perpendicular.vg surject/perpendicular.gam > perpendicular.sam
vg surject -U -s -x surject/perpendicular.vg surject/perpendicular.gam > perpendicular.sam
is "$?" 0 "vg surject does not crash when surjecting a read that grazes the reference with a deletion"
is "$(cat perpendicular.sam | grep -v "^@" | cut -f2)" "4" "vg surject leaves a read that grazes the reference with a deletion unmapped"

vg surject -s --prune-low-cplx -x surject/perpendicular.vg surject/perpendicular.gam > perpendicular.sam
vg surject -U -s --prune-low-cplx -x surject/perpendicular.vg surject/perpendicular.gam > perpendicular.sam
is "$?" 0 "vg surject does not crash when surjecting a read that grazes the reference with a deletion and pruning low complexity anchors"
is "$(cat perpendicular.sam | grep -v "^@" | cut -f2)" "4" "vg surject leaves a read that grazes the reference with a deletion unmapped when pruning low complexity anchors"

Expand All @@ -209,12 +209,22 @@ vg sim -x x.xg -n 20 -l 40 -p 60 -v 10 -a --random-seed 123 > x.gam
vg mpmap -x x.xg -g x.gcsa -n dna --suppress-mismapping -B -G x.gam -i -F GAM -I 60 -D 10 -t 1 > mapped.gam
vg mpmap -x x.xg -g x.gcsa -n dna --suppress-mismapping -B -G x.gam -i -F GAMP -I 60 -D 10 -t 1 > mapped.gamp

is "$(vg surject -x x.xg -s -t 1 mapped.gam | grep -v '@' | wc -l)" 40 "GAM surject can return only primaries"
is "$(vg surject -x x.xg -M -s -t 1 mapped.gam | grep -v '@' | wc -l)" 80 "GAM surject can return multimappings"
vg surject -x x.xg -s -t 1 mapped.gam >/dev/null 2>/dev/null
is "$?" "1" "GAM surject produces an error when given paired reads when not in interleaved mode"

vg surject -x x.xg -m -s -t 1 mapped.gamp >/dev/null 2>/dev/null
is "$?" "1" "GAMP surject produces an error when given paired reads when not in interleaved mode"

vg convert x.xg -G mapped.gam -t 1 >mapped.gaf
vg surject -x x.xg -G -s -t 1 mapped.gaf >/dev/null 2>/dev/null
is "$?" 1 "GAF surject produces an error when given paired reads when not in interleaved mode"

is "$(vg surject -x x.xg -U -s -t 1 mapped.gam | grep -v '@' | wc -l)" 40 "GAM surject can return only primaries"
is "$(vg surject -x x.xg -M -U -s -t 1 mapped.gam | grep -v '@' | wc -l)" 80 "GAM surject can return multimappings"
is "$(vg surject -x x.xg -M -i -s -t 1 mapped.gam | grep -v '@' | wc -l)" 80 "GAM surject can return paired multimappings"
is "$(vg surject -x x.xg -s -m -t 1 mapped.gamp | grep -v '@' | wc -l)" 40 "GAMP surject can return only primaries"
is "$(vg surject -x x.xg -M -m -s -t 1 mapped.gamp | grep -v '@' | wc -l)" 80 "GAMP surject can return multimappings"
is "$(vg surject -x x.xg -M -m -s -i -t 1 mapped.gamp | grep -v '@' | wc -l)" 80 "GAMP surject can return multimappings"
is "$(vg surject -x x.xg -U -s -m -t 1 mapped.gamp | grep -v '@' | wc -l)" 40 "GAMP surject can return only primaries"
is "$(vg surject -x x.xg -M -U -m -s -t 1 mapped.gamp | grep -v '@' | wc -l)" 80 "GAMP surject can return multimappings"
is "$(vg surject -x x.xg -M -i -m -s -i -t 1 mapped.gamp | grep -v '@' | wc -l)" 80 "GAMP surject can return paired multimappings"

vg construct -r tiny/tiny.fa > tiny.vg
vg surject -x tiny.vg -s -t 1 mapped.gam >/dev/null 2>err.txt
Expand Down
Loading
Loading