{

_max_distance = max_distance;

_max_read_frac_neighbor = max_read_frac_neighbor;

{

_graph.reserve(H.numReads());

// Assign node properties

for (auto i = 0; i < H.numReads(); ++i)

{

auto c = H.getRead(i).count();

_graph.emplace_back(i, c); // NEW: add read count

}

// Compute distances and add a link only if it is

// bounded by _max_distance

for (auto iter = H._HASHcounter.begin(); iter != H._HASHcounter.end(); ++iter)

{

auto start = *iter; // first pos. in _HASHvalues with the current fingerprint

auto jter = iter; ++jter;

auto end = (jter != H._HASHcounter.end()) ? *jter : H._HASHvalues.size(); // past-to-last pos. in _HASHvalues with the current fingerprint

for (auto i = start; i < end; ++i)

{

if (i >= H._HASHvalues.size()) std::cerr << "i >= H._HASHvalues.size() [" << i << "," << H._HASHvalues.size() << "]" << std::endl;

uint32_t r1 = H._HASHvalues.at(i).read_id;

uint32_t f1 = H._HASHvalues.at(i).block_id;

for (auto j = i+1; j < end; ++j)

{

if (j >= H._HASHvalues.size()) std::cerr << "j >= H._HASHvalues.size() [" << j << "," << H._HASHvalues.size() << "]" << std::endl;

uint32_t f2 = H._HASHvalues.at(j).block_id;

if (f1 != f2) continue; // only consider the sequences witn the same fingerprint on the same block

uint32_t r2 = H._HASHvalues.at(j).read_id;

uint32_t dist = H.ComputeHammingDistance (r1, r2);

if (dist > _max_distance) continue;

// if j no not exist among the neighbors of i yet, then add it and add the distance as well

// if j exists already, then check the distance

// if the distance is larger, then update it, otherwise do nothing

if (r1 >= _graph.size()) std::cerr << "r1 >= _graph.size() [" << r1 << "," << _graph.size() << "]" << std::endl;

auto n1 = _graph.at(r1).neighbors.begin();

auto d1 = _graph.at(r1).distances.begin();

while (n1 != _graph.at(r1).neighbors.end() and *n1 != r2)

{

++n1;

++d1;

}

if (r2 >= _graph.size()) std::cerr << "r2 >= _graph.size() [" << r2 << "," << _graph.size() << "]" << std::endl;

if (n1 == _graph.at(r1).neighbors.end())

{

_graph.at(r1).neighbors.push_back(r2);

_graph.at(r1).distances.push_back(dist);

_graph.at(r2).neighbors.push_back(r1);

_graph.at(r2).distances.push_back(dist);

}

else

{

if (*d1 <= dist) continue;

*d1 = dist;

auto n2 = _graph.at(r1).neighbors.begin();

auto d2 = _graph.at(r1).distances.begin();

while (n2 != _graph.at(r2).neighbors.end() and *n2 != r1)

{

++n2;

++d2;

}

if (n2 == _graph.at(r2).neighbors.end())

{

std::cerr << "the graph is not undirected!" << std::endl;

return;

}

*d2 = dist;

}

}

}

}

{

for (std::vector <Node>::iterator N = _graph.begin(); N != _graph.end(); ++N)

N->color = WHITE;

{

for (std::vector <Node>::iterator N = _graph.begin(); N != _graph.end(); ++N)

N->color = BLACK;

{

bool all_black = true;

for (std::vector <Node>::iterator N = _graph.begin(); N != _graph.end(); ++N)

all_black = (all_black and (N->color == BLACK));

return all_black;

{

std::vector <uint32_t> v;

v.reserve(_graph.size());

for (auto iter = _graph.begin(); iter != _graph.end(); ++iter)

v.push_back(iter->neighbors.size());

return v;

{

std::vector <uint32_t> v = compute_degree();

std::vector < std::vector <uint32_t> > w;

auto max_iter = std::max_element(v.begin(),v.end());

if (max_iter == v.end()) return w;

size_t max = *max_iter;

w.resize(max+1);

auto i = 0;

for (auto iter = v.begin(); iter != v.end(); ++iter)

{

if (w.size() == w.capacity())

w.reserve(w.size() + MAX_NEIGHBORS);

w.at(*iter).push_back(i);

++i;

}

return w;

{

std::vector <uint32_t> v = compute_degree();

std::vector < std::map <uint32_t, uint32_t> > w;

auto max_iter = std::max_element(v.begin(),v.end());

if (max_iter == v.end()) return w;

size_t max = *max_iter;

w.resize(max+1);

auto i = 0;

for (auto iter = v.begin(); iter != v.end(); ++iter)

{

w.at(*iter).insert(std::make_pair(i,i));

++i;

}

return w;

{

// node attributes

std::string outname(out_dir); outname.append("/node_attr.tsv");

std::ofstream outfile; outfile.open(outname.c_str());

if (not outfile.is_open()) { std::cerr << "Error while opening file " << outname << " for writing" << std::endl; return; }

outfile << "id num_neighbors counts" << std::endl;

for (auto ii = _graph.begin(); ii != _graph.end(); ++ii)

{

outfile << H.getRead(ii->node_idx).id() << " "

<< ii->neighbors.size() << " "

<< H.getRead(ii->node_idx).count() << std::endl;

}

outfile.close();

// edge attributes

outname = out_dir; outname.append("/edge_attr.tsv");

outfile.open(outname.c_str());

if (not outfile.is_open()) { std::cerr << "Error while opening file " << outname << " for writing" << std::endl; return; }

outfile << "id1 id2 count1 count2" << std::endl;

for (auto ii = _graph.begin(); ii != _graph.end(); ++ii)

{

for (auto jj = ii->neighbors.begin(); jj != ii->neighbors.end(); ++jj)

{

if (ii->node_idx < _graph.at(*jj).node_idx)

{

auto idx1 = (H.getRead(ii->node_idx).count() >= H.getRead(_graph.at(*jj).node_idx).count()) ? ii->node_idx : _graph.at(*jj).node_idx;

auto idx2 = (H.getRead(ii->node_idx).count() >= H.getRead(_graph.at(*jj).node_idx).count()) ? _graph.at(*jj).node_idx : ii->node_idx;

outfile << H.getRead(idx1).id() << " "

<< H.getRead(idx2).id() << " "

<< H.getRead(idx1).count() << " "

<< H.getRead(idx2).count() << std::endl;

}

}

}

outfile.close();

// neighbor attributes

outname = out_dir; outname.append("/neighbor_attr.tsv");

outfile.open(outname.c_str());

if (not outfile.is_open()) { std::cerr << "Error while opening file " << outname << " for writing" << std::endl; return; }

outfile << "id top_neigh_count top2_neigh_count" << std::endl;

for (auto ii = _graph.begin(); ii != _graph.end(); ++ii)

{

if (ii->neighbors.size() > 1)

{

// store the read count for all neighbors of ii

std::vector<uint32_t> neigh_counts;

neigh_counts.reserve(ii->neighbors.size());

for (auto jj = ii->neighbors.begin(); jj != ii->neighbors.end(); ++jj)

neigh_counts.push_back(H.getRead(*jj).count());

// find the top and the 2nd top read count neighbors

auto it_max_count = std::max_element(neigh_counts.begin(), neigh_counts.end());

uint32_t max_count = *it_max_count;

neigh_counts.erase(it_max_count);

auto it_sec_max_count = std::max_element(neigh_counts.begin(), neigh_counts.end());

uint32_t sec_max_count = *it_sec_max_count;

// check that the two neighbors of ii are not mutual neighbors

auto max_pos = it_max_count - neigh_counts.begin();

auto sec_max_pos = it_sec_max_count - neigh_counts.begin();

if (max_pos <= sec_max_pos) ++sec_max_pos;

auto jj = ii->neighbors.at(max_pos);

auto kk = ii->neighbors.at(sec_max_pos);

// print the non-mutual top count neighbors of ii

std::vector<uint32_t> nn = _graph.at(jj).neighbors;

if (std::find(nn.begin(), nn.end(), kk) == nn.end())

outfile << H.getRead(ii->node_idx).id() << " "

<< max_count << " "

<< sec_max_count << std::endl;

}

}

outfile.close();

// degree frequency

outname = out_dir; outname.append("/degree_freq.tsv");

outfile.open(outname.c_str());

if (not outfile.is_open()) { std::cerr << "Error while opening file " << outname << " for writing" << std::endl; return; }

outfile << "degree num_nodes" << std::endl;

std::vector < std::vector<uint32_t> > w = degree_sets();

auto i = 0;

for (auto iter = w.begin(); iter != w.end(); ++iter)

outfile << i++ << " " << iter->size() << std::endl;

outfile.close();

{

// set all nodes to WHITE = not visited

reset_color();

auto deg_nodes = compute_degree(); // idx = node idx, value = degree

auto degrees = degree_maps(); // idx = degree, value = map of nodes with the same degree (specified by idx)

// scan nodes by non-increasing degree order

uint32_t stars = 0;

uint32_t out_of_stars = 0; // initialize the count of the edges linking two distict stars

for (auto nodes = degrees.rbegin(); nodes != degrees.rend(); ++nodes)

{

while (not nodes->empty())

{

auto ii = nodes->begin()->first; // first node in the map

if (ii >= _graph.size())

{

std::cerr << "Node " << ii << " is outside of node idx span" << std::endl;

continue;

}

if (_graph.at(ii).color == BLACK) continue;

// visit node ii

_graph.at(ii).color = GRAY;

// visit ii's neighbors jj

std::vector<uint32_t> ne;

ne.reserve(_graph.at(ii).neighbors.size());

for (auto jj = _graph.at(ii).neighbors.begin(); jj != _graph.at(ii).neighbors.end(); ++jj)

{

if (_graph.at(*jj).color == BLACK) continue; // the edge (ii,jj) has already been counted in out_of_stars

_graph.at(*jj).color = GRAY;

_graph.at(*jj).parent = ii; // THIS CONTAINS THE STAR/GROUP INFORMATION!!!

ne.push_back(*jj);

}

// scan the neighbors of jj

std::vector <uint32_t> nne;

nne.reserve(MAX_NEIGHBORS);

for (auto jj = ne.begin(); jj != ne.end(); ++jj)

{

auto v = _graph.at(*jj).neighbors;

for (auto kk = v.begin(); kk != v.end(); ++kk)

{

// new neighbor

if (_graph.at(*kk).color == WHITE and std::find(nne.begin(), nne.end(), *kk) == nne.end())

nne.push_back(*kk);

}

}

out_of_stars += nne.size();

// update the degree of the neighbors of the neighbors

for (auto kk = nne.begin(); kk != nne.end(); ++kk)

{

if (*kk >= deg_nodes.size())

{

std::cerr << "Node " << *kk << " is outside of deg nodes span" << std::endl;

continue;

}

auto current_degree = deg_nodes.at(*kk);

uint32_t edges_to_remove = 0;

auto v = _graph.at(*kk).neighbors;

for (auto hh = v.begin(); hh != v.end(); ++hh)

edges_to_remove += (_graph.at(*hh).color == GRAY);

if (edges_to_remove > current_degree)

{

std::cerr << "The number of edges to remove from node " << *kk << " is greater than its degree ("

<< edges_to_remove << " > " << current_degree << ")" << std::endl;

continue;

}

auto new_degree = current_degree - edges_to_remove;

// find the current slot and remove kk

auto nodes_n = degrees.begin() + current_degree;

auto iter = nodes_n->find(*kk);

if (iter == nodes_n->end())

{

std::cerr << "Node " << *kk << " was not found on degree slot " << current_degree << std::endl;

continue;

}

nodes_n->erase(iter);

// find the new slot and add kk to it

nodes_n = degrees.begin() + new_degree;

nodes_n->insert(std::make_pair(*kk,*kk));

deg_nodes.at(*kk) = new_degree; // IMPORTANT!! otherwise the node won't be found in the correct degrees slot

}

// remove ii from the degree slot

nodes->erase(nodes->begin());

// remove ii's neighbors from their degree slot

for (auto jj = ne.begin(); jj != ne.end(); ++jj)

{

if (*jj >= deg_nodes.size())

{

std::cerr << "Node " << *jj << " is outside of deg nodes span" << std::endl;

continue;

}

if (deg_nodes.at(*jj) >= degrees.size())

{

std::cerr << "Node " << *jj << " degree is greater than maximum degree ("

<< deg_nodes.at(*jj) << " >= " << degrees.size() << ")" << std::endl;

}

auto nodes_n = degrees.begin() + deg_nodes.at(*jj); // iterator to the map containing the nodes with degree = degree of j

auto iter = nodes_n->find(*jj);

if (iter == nodes_n->end())

{

std::cerr << "Node " << *jj << " was not found on degree slot " << deg_nodes.at(*jj) << std::endl;

continue;

}

nodes_n->erase(iter);

}

// end the visit of ii and jj

// N.B.: kk are NOT marked as visited yet

_graph.at(ii).color = BLACK;

for (auto jj = ne.begin(); jj != ne.end(); ++jj)

_graph.at(*jj).color = BLACK;

++stars;

}

}

std::cout << "number of groups: " << stars << std::endl;

std::cout << "number of edges connecting two groups: " << out_of_stars << std::endl;

{

std::map < uint32_t, std::vector <uint32_t> > m;

for (auto ii = _graph.begin(); ii != _graph.end(); ++ii)

{

auto cc = H.getRead(ii->node_idx).count();

if (m.find(cc) == m.end())

{

std::vector <uint32_t> v;

m.insert(std::make_pair(cc, v));

}

m.find(cc)->second.push_back(ii->node_idx);

}

return m;

{

// set all nodes to WHITE = not visited

reset_color();

auto counts = count_map(H); // first = count, second = vector of node indexes

// scan nodes by non-increasing degree order

uint32_t stars = 0;

uint32_t out_of_stars = 0; // initialize the count of the edges linking two distict stars

for (auto cc = counts.rbegin(); cc != counts.rend(); ++cc)

{

// scan node indexes

for (auto ii = cc->second.begin(); ii != cc->second.end(); ++ii)

{

if (*ii >= _graph.size())

{

std::cerr << "Node " << *ii << " is outside of node idx span" << std::endl;

continue;

}

if (_graph.at(*ii).color == BLACK) continue;

// visit node ii

_graph.at(*ii).color = GRAY;

// visit ii's neighbors jj

std::vector<uint32_t> ne;

ne.reserve(_graph.at(*ii).neighbors.size());

for (auto jj = _graph.at(*ii).neighbors.begin(); jj != _graph.at(*ii).neighbors.end(); ++jj)

{

if (_graph.at(*jj).color == BLACK) continue; // the edge (ii,jj) has already been counted in out_of_stars

// NEW: check that the count of the neighbor is below _max_read_frac_neighbor fraction

if (_max_read_frac_neighbor * _graph.at(*ii).count >= _graph.at(*jj).count)

{

_graph.at(*jj).color = GRAY;

_graph.at(*jj).parent = *ii; // THIS CONTAINS THE STAR/GROUP INFORMATION!!!

ne.push_back(*jj);

} // else {

}

// scan the neighbors of jj

std::vector <uint32_t> nne;

nne.reserve(MAX_NEIGHBORS);

for (auto jj = ne.begin(); jj != ne.end(); ++jj)

{

auto v = _graph.at(*jj).neighbors;

for (auto kk = v.begin(); kk != v.end(); ++kk)

{

// new neighbor

if (_graph.at(*kk).color == WHITE and std::find(nne.begin(), nne.end(), *kk) == nne.end())

nne.push_back(*kk);

}

}

out_of_stars += nne.size();

// end the visit of ii and jj

// N.B.: kk are NOT marked as visited yet

_graph.at(*ii).color = BLACK;

for (auto jj = ne.begin(); jj != ne.end(); ++jj)

_graph.at(*jj).color = BLACK;

++stars;

}

}

std::cout << "number of groups: " << stars << std::endl;

std::cout << "number of edges connecting two groups: " << out_of_stars << std::endl;

{

// collect stars

std::map < uint32_t, std::map <uint32_t, uint32_t> > groups; // key = parent idx, value = (idx,idx)

for (auto ii = _graph.begin(); ii != _graph.end(); ++ii)

{

if (groups.find(ii->parent) == groups.end())

{

std::map <uint32_t, uint32_t> m;

groups.insert(std::make_pair(ii->parent, m));

}

auto map_iter = groups.find(ii->parent);

map_iter->second.insert(std::make_pair(ii->node_idx,ii->node_idx));

}

// sort stars

std::map < uint32_t, std::vector <uint32_t> > groups_sorted; // key = group count, value = list of parents identifying the groups

for (auto ii = groups.begin(); ii != groups.end(); ++ii)

{

auto c = 0;

for (auto jj = ii->second.begin(); jj != ii->second.end(); ++jj)

c += H.getRead(jj->first).count();

if (groups_sorted.find(c) == groups_sorted.end())

{

std::vector <uint32_t> v;

groups_sorted.insert(std::make_pair(c,v));

}

groups_sorted.find(c)->second.push_back(ii->first);

}

// select min group count

uint32_t c;

bool found = false;

for (auto ii = groups_sorted.rbegin(); ii != groups_sorted.rend(); ++ii)

{

for (auto jj = ii->second.begin(); jj != ii->second.end(); ++jj)

{

double tot_degree = _graph.at(*jj).neighbors.size();

double in_degree = groups.find(*jj)->second.size() - 1;

found = (in_degree < MIN_FRAC_IN_DEGREE * tot_degree);

if (found) break;

}

if (found) break;

c = ii->first;

}

std::cout << "Min group count: " << c << std::endl;

// compute group size frequency

uint32_t max_size_all = 0;

for (auto ii = groups.begin(); ii != groups.end(); ++ii)

{

if (max_size_all < ii->second.size())

max_size_all = ii->second.size();

}

std::vector <uint32_t> group_sizes_all; group_sizes_all.resize(max_size_all, 0);

for (auto ii = groups.begin(); ii != groups.end(); ++ii)

{

if (ii->second.size()-1 >= group_sizes_all.size())

{

std::cerr << "ii->second.size()-1 >= group_sizes_all.size() (" << (ii->second.size()-1) << " > " << group_sizes_all.size() << ")" << std::endl;

return;

}

group_sizes_all.at(ii->second.size()-1)++;

}

// compute selected group size frequency

uint32_t max_size = 0;

for (auto ii = groups_sorted.rbegin(); ii != groups_sorted.rend(); ++ii)

{

if (ii->first < c) break;

for (auto jj = ii->second.begin(); jj != ii->second.end(); ++jj)

{

auto kk = groups.find(*jj);

if (max_size < kk->second.size())

max_size = kk->second.size();

}

}

std::vector <uint32_t> group_sizes; group_sizes.resize(max_size, 0);

for (auto ii = groups_sorted.rbegin(); ii != groups_sorted.rend(); ++ii)

{

if (ii->first < c) break;

for (auto jj = ii->second.begin(); jj != ii->second.end(); ++jj)

{

auto kk = groups.find(*jj);

if (kk != groups.end())

{

if (kk->second.size()-1 >= group_sizes.size())

{

std::cerr << "kk->second.size()-1 >= group_sizes.size() (" << kk->second.size()-1 << " > " << group_sizes.size() << ")" << std::endl;

return;

}

group_sizes.at(kk->second.size()-1)++;

}

}

}

// print groups: all and selected

std::string outname(out_dir); outname.append("/groups_all.tsv");

std::ofstream outfile_all; outfile_all.open(outname.c_str());

if (not outfile_all.is_open()) { std::cerr << "Error while opening file " << outname << " for writing" << std::endl; return; }

outname = out_dir; outname.append("/groups.tsv");

std::ofstream outfile; outfile.open(outname.c_str());

if (not outfile.is_open()) { std::cerr << "Error while opening file " << outname << " for writing" << std::endl; return; }

outfile_all << "group size hub hub_count group_count out_hub_degree" << std::endl;

outfile << "group size hub hub_count group_count out_hub_degree" << std::endl;

size_t group_id = 1;

for (auto ii = groups_sorted.rbegin(); ii != groups_sorted.rend(); ++ii)

{

for (auto jj = ii->second.begin(); jj != ii->second.end(); ++jj)

{

std::ostringstream ss;

ss << group_id++ << " "

<< groups.find(*jj)->second.size() << " "

<< H.getRead(*jj).id() << " "

<< H.getRead(*jj).count() << " "

<< ii->first << " "

<< (_graph.at(*jj).neighbors.size() - (groups.find(*jj)->second.size() - 1));

outfile_all << ss.str() << std::endl;

if (ii->first >= c) outfile << ss.str() << std::endl;

}

}

outfile.close();

// print node and parent information (parent = hub of the group, node = element of the group)

outname = out_dir; outname.append("/groups_info.tsv");

outfile.open(outname.c_str());

if (not outfile.is_open()) { std::cerr << "Error while opening file " << outname << " for writing" << std::endl; return; }

outfile << "node hub dist count" << std::endl;

for (auto ii = _graph.begin(); ii != _graph.end(); ++ii)

{

outfile << H.getRead(ii->node_idx).id() << " "

<< H.getRead(ii->parent).id() << " "

<< H.ComputeHammingDistance(ii->node_idx, ii->parent) << " "

<< H.getRead(ii->node_idx).count() << std::endl;

}

outfile.close();

// collect the edges connecting nodes belonging to different groups

outname = out_dir; outname.append("/out_of_groups.tsv");

outfile.open(outname.c_str());

if (not outfile.is_open()) { std::cerr << "Error while opening file " << outname << " for writing" << std::endl; return; }

outfile << "node1 node2 hub1 hub2" << std::endl;

for (auto ii = _graph.begin(); ii != _graph.end(); ++ii)

{

for (auto jj = ii->neighbors.begin(); jj != ii->neighbors.end(); ++jj)

{

if (ii->node_idx < *jj) continue;

if (ii->parent != _graph.at(*jj).parent)

outfile << H.getRead(ii->node_idx).id() << " "

<< H.getRead(*jj).id() << " "

<< H.getRead(ii->parent).id() << " "

<< H.getRead(_graph.at(*jj).parent).id() << std::endl;

}

}

outfile.close();

// all groups size frequency

outname = out_dir; outname.append("/group_all_size_freq.tsv");

outfile.open(outname.c_str());

if (not outfile.is_open()) { std::cerr << "Error while opening file " << outname << " for writing" << std::endl; return; }

outfile << "size num_groups" << std::endl;

size_t i = 1;

for (auto ii = group_sizes_all.begin(); ii != group_sizes_all.end(); ++ii)

outfile << i++ << " " << *ii << std::endl;

outfile.close();

// selected groups size frequency

outname = out_dir; outname.append("/group_size_freq.tsv");

outfile.open(outname.c_str());

if (not outfile.is_open()) { std::cerr << "Error while opening file " << outname << " for writing" << std::endl; return; }

outfile << "size num_groups" << std::endl;

i = 1;

for (auto ii = group_sizes.begin(); ii != group_sizes.end(); ++ii)

outfile << i++ << " " << *ii << std::endl;

outfile.close();