Skip to content
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

Deconstruct option to cluster similar alleles together #4301

Merged
merged 4 commits into from
May 29, 2024
Merged

Deconstruct option to cluster similar alleles together #4301

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
3e766a7
start trav clustering; move jaccard code alongside
glennhickey May 22, 2024
cb1c4d1
start traversal-clustering logic
glennhickey May 23, 2024
9d01127
add some tags for alt cluster info
glennhickey May 24, 2024
f491fdb
move to weighted jac; add small test
glennhickey May 27, 2024
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
117 changes: 69 additions & 48 deletions src/deconstructor.cpp
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading