Skip to content

Commit

Permalink
Merge pull request #4301 from vgteam/deconstruct
Browse files Browse the repository at this point in the history 
Deconstruct option to cluster similar alleles together
  • Loading branch information
@glennhickey
glennhickey authored May 29, 2024
2 parents db574a5 + f491fdb commit 2fea419
Show file tree
Hide file tree
Showing 8 changed files with 370 additions and 85 deletions.
117 changes: 69 additions & 48 deletions src/deconstructor.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#include "deconstructor.hpp"
#include "traversal_finder.hpp"
#include <gbwtgraph/gbwtgraph.h>
#include "traversal_clusters.hpp"

//#define debug

Expand All @@ -25,7 +26,7 @@ vector<int> Deconstructor::get_alleles(vcflib::Variant& v,
const vector<Traversal>& travs,
const vector<pair<step_handle_t, step_handle_t>>& trav_steps,
int ref_path_idx,
const vector<bool>& use_trav,
const vector<vector<int>>& trav_clusters,
char prev_char, bool use_start) const {

assert(ref_path_idx >=0 && ref_path_idx < travs.size());
Expand Down Expand Up @@ -56,10 +57,11 @@ vector<int> Deconstructor::get_alleles(vcflib::Variant& v,
bool substitution = true;

// set the other alleles (they can end up as 0 alleles too if their strings match the reference)
for (int i = 0; i < travs.size(); ++i) {
if (i != ref_path_idx) {
if (use_trav[i]) {
string allele = trav_to_string(travs[i]);
// note that we have one (unique) allele per cluster, so we take advantage of that here
for (const vector<int>& cluster : trav_clusters) {
string allele = trav_to_string(travs[cluster.front()]);
for (const int& i : cluster) {
if (i != ref_path_idx) {
auto ai_it = allele_idx.find(allele);
if (ai_it == allele_idx.end()) {
// make a new allele for this string
Expand Down Expand Up @@ -248,7 +250,7 @@ vector<int> Deconstructor::get_alleles(vcflib::Variant& v,
for (auto& c : ref_contexts) {
auto& ref_context = c.second;
auto& ref_pos = c.first;
double j = context_jaccard(ref_context, path_context);
double j = jaccard_coefficient(ref_context, path_context);
if (j > best_jaccard) {
best_jaccard = j;
best_pos = ref_pos;
Expand Down Expand Up @@ -286,13 +288,18 @@ vector<int> Deconstructor::get_alleles(vcflib::Variant& v,
}

void Deconstructor::get_genotypes(vcflib::Variant& v, const vector<string>& names,
const vector<int>& trav_to_allele) const {
const vector<int>& trav_to_allele,
const vector<pair<double, int64_t>>& trav_to_cluster_info) const {
assert(names.size() == trav_to_allele.size());
// set up our variant fields
v.format.push_back("GT");
if (show_path_info && path_to_sample_phase) {
v.format.push_back("PI");
}
if (this->cluster_threshold < 1.0) {
v.format.push_back("TS");
v.format.push_back("TL");
}

// get a list of traversals for every vcf sample
// (this will be 1:1 unless we're using the path_to_sample name map)
Expand Down Expand Up @@ -342,6 +349,18 @@ void Deconstructor::get_genotypes(vcflib::Variant& v, const vector<string>& name
}
genotype += (chosen_travs[i] != -1 && (!conflict || keep_conflicted_genotypes))
? std::to_string(trav_to_allele[chosen_travs[i]]) : ".";
if (this->cluster_threshold < 1.0) {
if (*genotype.rbegin() == '.') {
v.samples[sample_name]["TS"].push_back(".");
v.samples[sample_name]["TL"].push_back(".");
} else {
ostringstream ss;
ss.precision(3);
ss << trav_to_cluster_info[chosen_travs[i]].first;
v.samples[sample_name]["TS"].push_back(ss.str());
v.samples[sample_name]["TL"].push_back(std::to_string(trav_to_cluster_info[chosen_travs[i]].second));
}
}
}
v.samples[sample_name]["GT"] = {genotype};
if (show_path_info && path_to_sample_phase) {
Expand Down Expand Up @@ -487,38 +506,6 @@ pair<vector<int>, bool> Deconstructor::choose_traversals(const string& sample_na
return make_pair(most_frequent_travs, conflict);
}


// todo refactor if we need to reuse elsewhere in vg
// implemented inline for development
// assumes sorted input
double Deconstructor::context_jaccard(
const vector<nid_t>& target,
const vector<nid_t>& query) const {
size_t node_isec = 0;
std::set_intersection(target.begin(), target.end(),
query.begin(), query.end(),
count_back_inserter<nid_t>(node_isec));
size_t node_union = 0;
std::set_union(target.begin(), target.end(),
query.begin(), query.end(),
count_back_inserter<nid_t>(node_union));
return (double)node_isec / (double)node_union;
}

double Deconstructor::context_jaccard(
const dac_vector<>& target,
const vector<nid_t>& query) const {
size_t node_isec = 0;
std::set_intersection(target.begin(), target.end(),
query.begin(), query.end(),
count_back_inserter<nid_t>(node_isec));
size_t node_union = 0;
std::set_union(target.begin(), target.end(),
query.begin(), query.end(),
count_back_inserter<nid_t>(node_union));
return (double)node_isec / (double)node_union;
}

vector<nid_t> Deconstructor::get_context(
step_handle_t start_step,
step_handle_t end_step) const {
Expand Down Expand Up @@ -639,15 +626,17 @@ bool Deconstructor::deconstruct_site(const handle_t& snarl_start, const handle_t

// pick out the traversal corresponding to an embedded reference path, breaking ties consistently
string ref_trav_name;
for (int i = 0; i < trav_steps.size(); ++i) {
for (int i = 0; i < travs.size(); ++i) {
const string& path_trav_name = trav_path_names[i];
#ifdef debug
#pragma omp critical (cerr)
{
cerr << "Traversal " << i << ": name=" << path_trav_name << ", size=" << travs[i].size()
<< ", start=" << graph->get_position_of_step(trav_steps[i].first)
<< ", end=" << graph->get_position_of_step(trav_steps[i].second) << endl
<< " trav=" << traversal_to_string(travs[i]) << endl;
cerr << "Traversal " << i << ": name=" << path_trav_name << ", size=" << travs[i].size();
if (i < trav_steps.size()) {
cerr << ", start=" << graph->get_position_of_step(trav_steps[i].first)
<< ", end=" << graph->get_position_of_step(trav_steps[i].second) << endl;
}
cerr << " trav=" << traversal_to_string(graph, travs[i]) << endl;
}
#endif
if (ref_paths.count(path_trav_name) &&
Expand Down Expand Up @@ -738,7 +727,7 @@ bool Deconstructor::deconstruct_site(const handle_t& snarl_start, const handle_t
vector<pair<double, int>> ref_mappings;
for (uint64_t j = 0; j < ref_travs.size(); ++j) {
ref_mappings.push_back(make_pair(
context_jaccard(
jaccard_coefficient(
ref_contexts[j],
context),
ref_travs[j]));
Expand Down Expand Up @@ -828,13 +817,36 @@ bool Deconstructor::deconstruct_site(const handle_t& snarl_start, const handle_t
}
}

// Sort the traversals for clustering
vector<int> sorted_travs = get_traversal_order(graph, travs, trav_path_names, ref_travs, use_trav);

// jaccard clustering (using handles for now) on traversals
vector<pair<double, int64_t>> trav_cluster_info;
vector<vector<int>> trav_clusters = cluster_traversals(graph, travs, sorted_travs, cluster_threshold,
trav_cluster_info);

#ifdef debug
cerr << "cluster priority";
for (const auto& t: sorted_travs) {
cerr << " " << t;
}
cerr << endl;
for (const auto& tc : trav_clusters) {
cerr << "traversal cluster:";
for (const auto& t: tc) {
cerr << t << "(" << trav_cluster_info[t].first << "," << trav_cluster_info[t].second << ") ";
}
cerr << endl;
}
#endif

vector<int> trav_to_allele = get_alleles(v, travs, trav_steps,
ref_trav_idx,
use_trav,
trav_clusters,
prev_char, use_start);

// Fill in the genotypes
get_genotypes(v, trav_path_names, trav_to_allele);
get_genotypes(v, trav_path_names, trav_to_allele, trav_cluster_info);

// we only bother printing out sites with at least 1 non-reference allele
if (!std::all_of(trav_to_allele.begin(), trav_to_allele.end(), [](int i) { return (i == 0 || i == -1); })) {
Expand Down Expand Up @@ -913,7 +925,7 @@ string Deconstructor::get_vcf_header() {
if (sample_to_haps.empty()) {
cerr << "Error [vg deconstruct]: No paths found for alt alleles in the graph. Note that "
<< "exhaustive path-free traversal finding is no longer supported, and vg deconstruct "
<< "now only works on embedded paths and GBWT threads" << endl;
<< "now only works on embedded paths and GBWT threads." << endl;
exit(1);
}

Expand All @@ -939,6 +951,13 @@ string Deconstructor::get_vcf_header() {
}
stream << "\">" << endl;
}
if (path_to_sample_phase && cluster_threshold < 1) {
stream << "##FORMAT<ID=TS,Number=1,Type=Float,Descript=\"Similarity between the sample's actual path and its allele\">"
<< endl;
stream << "##FORMAT<ID=TL,Number=1,Type=Integer,Descript=\"Length difference between the sample's actual path and its allele\">"
<< endl;

}
stream << "##INFO=<ID=AC,Number=A,Type=Integer,Description=\"Total number of alternate alleles in called genotypes\">" << endl;
stream << "##INFO=<ID=AF,Number=A,Type=Float,Description=\"Estimated allele frequency in the range (0,1]\">" << endl;
stream << "##INFO=<ID=NS,Number=1,Type=Integer,Description=\"Number of samples with data\">" << endl;
Expand Down Expand Up @@ -1105,6 +1124,7 @@ void Deconstructor::deconstruct(vector<string> ref_paths, const PathPositionHand
bool keep_conflicted,
bool strict_conflicts,
bool long_ref_contig,
double cluster_threshold,
gbwt::GBWT* gbwt) {

this->graph = graph;
Expand All @@ -1117,6 +1137,7 @@ void Deconstructor::deconstruct(vector<string> ref_paths, const PathPositionHand
if (gbwt) {
this->gbwt_reference_samples = gbwtgraph::parse_reference_samples_tag(*gbwt);
}
this->cluster_threshold = cluster_threshold;
this->gbwt = gbwt;

// the need to use nesting is due to a problem with omp tasks and shared state
Expand Down
35 changes: 9 additions & 26 deletions src/deconstructor.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,7 @@ class Deconstructor : public VCFOutputCaller {
bool keep_conflicted,
bool strict_conflicts,
bool long_ref_contig,
double cluster_threshold = 1.0,
gbwt::GBWT* gbwt = nullptr);

private:
Expand Down Expand Up @@ -78,11 +79,12 @@ class Deconstructor : public VCFOutputCaller {
const vector<Traversal>& travs,
const vector<pair<step_handle_t, step_handle_t>>& trav_steps,
int ref_path_idx,
const vector<bool>& use_trav,
const vector<vector<int>>& trav_clusters,
char prev_char, bool use_start) const;

// write traversal path names as genotypes
void get_genotypes(vcflib::Variant& v, const vector<string>& names, const vector<int>& trav_to_allele) const;
void get_genotypes(vcflib::Variant& v, const vector<string>& names, const vector<int>& trav_to_allele,
const vector<pair<double, int64_t>>& trav_to_cluster_info) const;

// given a set of traversals associated with a particular sample, select a set of size <ploidy> for the VCF
// the highest-frequency ALT traversal is chosen
Expand All @@ -96,15 +98,6 @@ class Deconstructor : public VCFOutputCaller {
vector<nid_t> get_context(
step_handle_t start_step,
step_handle_t end_step) const;

// compares node contexts
double context_jaccard(const vector<nid_t>& target,
const vector<nid_t>& query) const;

// specialization for enc_vectors
double context_jaccard(
const dac_vector<>& target,
const vector<nid_t>& query) const;

// the graph
const PathPositionHandleGraph* graph;
Expand Down Expand Up @@ -154,23 +147,13 @@ class Deconstructor : public VCFOutputCaller {

// should we keep conflicted genotypes or not
bool keep_conflicted_genotypes = false;
};

// helpel for measuring set intersectiond and union size
template <typename T>
class count_back_inserter {
size_t &count;
public:
typedef void value_type;
typedef void difference_type;
typedef void pointer;
typedef void reference;
typedef std::output_iterator_tag iterator_category;
count_back_inserter(size_t &count) : count(count) {};
void operator=(const T &){ ++count; }
count_back_inserter &operator *(){ return *this; }
count_back_inserter &operator++(){ return *this; }
// used to merge together similar traversals (to keep allele counts down)
// currently implemented as handle jaccard coefficient. So 1 means only
// merge if identical (which is what deconstruct has always done)
double cluster_threshold = 1.0;
};


}
#endif
32 changes: 24 additions & 8 deletions src/subcommand/deconstruct_main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,7 @@ void help_deconstruct(char** argv){
<< " -K, --keep-conflicted Retain conflicted genotypes in output." << endl
<< " -S, --strict-conflicts Drop genotypes when we have more than one haplotype for any given phase (set by default when using GBWT input)." << endl
<< " -C, --contig-only-ref Only use the CONTIG name (and not SAMPLE#CONTIG#HAPLOTYPE etc) for the reference if possible (ie there is only one reference sample)." << endl
<< " -L, --cluster F Cluster traversals whose (handle) Jaccard coefficient is >= F together (default: 1.0) [experimental]" << endl
<< " -t, --threads N Use N threads" << endl
<< " -v, --verbose Print some status messages" << endl
<< endl;
Expand All @@ -77,6 +78,7 @@ int main_deconstruct(int argc, char** argv){
int context_jaccard_window = 10000;
bool untangle_traversals = false;
bool contig_only_ref = false;
double cluster_threshold = 1.0;

int c;
optind = 2; // force optind past command positional argument
Expand All @@ -97,14 +99,15 @@ int main_deconstruct(int argc, char** argv){
{"all-snarls", no_argument, 0, 'a'},
{"keep-conflicted", no_argument, 0, 'K'},
{"strict-conflicts", no_argument, 0, 'S'},
{"contig-only-ref", no_argument, 0, 'C'},
{"contig-only-ref", no_argument, 0, 'C'},
{"cluster", required_argument, 0, 'L'},
{"threads", required_argument, 0, 't'},
{"verbose", no_argument, 0, 'v'},
{0, 0, 0, 0}
};

int option_index = 0;
c = getopt_long (argc, argv, "hp:P:H:r:g:T:OeKSCd:c:uat:v",
c = getopt_long (argc, argv, "hp:P:H:r:g:T:OeKSCd:c:uaL:t:v",
long_options, &option_index);

// Detect the end of the options.
Expand Down Expand Up @@ -158,7 +161,10 @@ int main_deconstruct(int argc, char** argv){
break;
case 'C':
contig_only_ref = true;
break;
break;
case 'L':
cluster_threshold = max(0.0, min(1.0, parse<double>(optarg)));
break;
case 't':
omp_set_num_threads(parse<int>(optarg));
break;
Expand Down Expand Up @@ -254,13 +260,22 @@ int main_deconstruct(int argc, char** argv){
}

if (refpaths.empty() && refpath_prefixes.empty()) {
bool found_hap;
// No paths specified: use them all
graph->for_each_path_handle([&](path_handle_t path_handle) {
const string& name = graph->get_path_name(path_handle);
if (!Paths::is_alt(name) && PathMetadata::parse_sense(name) != PathSense::HAPLOTYPE) {
refpaths.push_back(name);
}
});
const string& name = graph->get_path_name(path_handle);
if (!Paths::is_alt(name) && PathMetadata::parse_sense(name) != PathSense::HAPLOTYPE) {
refpaths.push_back(name);
} else {
found_hap = true;
}
});

if (!found_hap) {
cerr << "error [vg deconstruct]: All graph paths selected as references (leaving no alts). Please use -P/-p "
<< "to narrow down the reference to a subset of paths, or GBZ/GBWT input that contains haplotype paths" << endl;
return 1;
}
}

// Read the translation
Expand Down Expand Up @@ -344,6 +359,7 @@ int main_deconstruct(int argc, char** argv){
keep_conflicted,
strict_conflicts,
!contig_only_ref,
cluster_threshold,
gbwt_index);
return 0;
}
Expand Down
Loading

1 comment on commit 2fea419

@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.

16 tests passed, 0 tests failed and 0 tests skipped in 17352 seconds

Please sign in to comment.