diff --git a/scripts/scatterplot8.py b/scripts/scatterplot8.py index bee4280..66c9488 100755 --- a/scripts/scatterplot8.py +++ b/scripts/scatterplot8.py @@ -121,16 +121,17 @@ def plot_data(fig, ax, subplot_coords, x, y, c, query_length, ymin, ymax, l_medi all_colors = 'bgrcmyk'; + colors1 = [all_colors[val%len(all_colors)] for val in c[0::2]]; + colors2 = [all_colors[val%len(all_colors)] for val in c[1::2]]; + i = 0; while (i < len(x)): - ax.plot(x[i:(i+2)], y[i:(i+2)], 'k'); + ax.plot(x[i:(i+2)], y[i:(i+2)], color=all_colors[(c[i]) % len(all_colors)]); i += 2; # ax.plot(x, y, 'o'); # ax.scatter(x[0::2], y[0::2], s=10, facecolor='b', lw = 0.2) # ax.scatter(x[1::2], y[1::2], s=10, facecolor='c', lw = 0.2) - colors1 = [all_colors[val%len(all_colors)] for val in c[0::2]]; - colors2 = [all_colors[val%len(all_colors)] for val in c[1::2]]; ax.scatter(x[0::2], y[0::2], s=10, edgecolor=colors1, facecolor=colors1, lw = 0.2) ax.scatter(x[1::2], y[1::2], s=10, edgecolor=colors2, facecolor=colors2, lw = 0.2) diff --git a/src/alignment/alignment.cc b/src/alignment/alignment.cc index 118bfc4..ce2c8af 100644 --- a/src/alignment/alignment.cc +++ b/src/alignment/alignment.cc @@ -9,9 +9,11 @@ #include "libs/opal.h" int AlignRegion(const SingleSequence *read, const Index *index, const ProgramParameters *parameters, const EValueParams *evalue_params, bool extend_to_end, PathGraphEntry *region_results) { - bool align_end_to_end = true; +// bool align_end_to_end = true; // bool spliced_alignment = true; - bool spliced_alignment = false; +// bool spliced_alignment = false; + bool align_end_to_end = parameters->alignment_algorithm != "spliced"; + bool spliced_alignment = parameters->alignment_algorithm == "spliced"; if (parameters->alignment_algorithm == "gotoh") { LogSystem::GetInstance().Log(VERBOSE_LEVEL_ALL_DEBUG, ((int64_t) read->get_sequence_id()) == parameters->debug_read, "Using semiglobal alignment approach.\n", "Alignment"); @@ -25,13 +27,13 @@ int AlignRegion(const SingleSequence *read, const Index *index, const ProgramPar return SemiglobalAlignment(MyersSemiglobalWrapper, read, index, parameters, evalue_params, region_results); - } else if (parameters->alignment_algorithm == "anchor") { + } else if (parameters->alignment_algorithm == "anchor" || parameters->alignment_algorithm == "spliced") { LogSystem::GetInstance().Log(VERBOSE_LEVEL_ALL_DEBUG, ((int64_t) read->get_sequence_id()) == parameters->debug_read, "Using anchored alignment approach.\n", "Alignment"); bool is_linear = region_results->get_region_data().is_split == false || parameters->is_reference_circular == false; LogSystem::GetInstance().Log(VERBOSE_LEVEL_ALL_DEBUG, ((int64_t) read->get_sequence_id()) == parameters->debug_read, "Using Myers' bit-vector algorithm for alignment!\n", "Alignment"); return AnchoredAlignmentNew(MyersNWWrapper, MyersSHWWrapper, read, index, parameters, evalue_params, region_results, align_end_to_end, spliced_alignment); - } else if (parameters->alignment_algorithm == "anchorgotoh") { + } else if (parameters->alignment_algorithm == "anchorgotoh" || parameters->alignment_algorithm == "splicedgotoh") { LogSystem::GetInstance().Log(VERBOSE_LEVEL_ALL_DEBUG, ((int64_t) read->get_sequence_id()) == parameters->debug_read, "Using anchored alignment approach.\n", "Alignment"); bool is_linear = region_results->get_region_data().is_split == false || parameters->is_reference_circular == false; LogSystem::GetInstance().Log(VERBOSE_LEVEL_ALL_DEBUG, ((int64_t) read->get_sequence_id()) == parameters->debug_read, "Using Gotoh's algorithm for alignment!\n", "Alignment"); @@ -39,7 +41,7 @@ int AlignRegion(const SingleSequence *read, const Index *index, const ProgramPar #ifndef RELEASE_VERSION - } else if (parameters->alignment_algorithm == "anchormex") { + } else if (parameters->alignment_algorithm == "anchormex" || parameters->alignment_algorithm == "splicedmex") { LogSystem::GetInstance().Log(VERBOSE_LEVEL_ALL_DEBUG, ((int64_t) read->get_sequence_id()) == parameters->debug_read, "Using anchored alignment approach.\n", "Alignment"); bool is_linear = region_results->get_region_data().is_split == false || parameters->is_reference_circular == false; LogSystem::GetInstance().Log(VERBOSE_LEVEL_ALL_DEBUG, ((int64_t) read->get_sequence_id()) == parameters->debug_read, "Using Match Extend algorithm for alignment!\n", "Alignment"); diff --git a/src/alignment/local_realignment_wrappers.cc b/src/alignment/local_realignment_wrappers.cc index 4651eb7..e00ee06 100644 --- a/src/alignment/local_realignment_wrappers.cc +++ b/src/alignment/local_realignment_wrappers.cc @@ -504,7 +504,7 @@ int MyersSemiglobalWrapper(const int8_t *read_data, int64_t read_length, int myers_return_code = myersCalcEditDistance((const unsigned char *) read_data, read_length, (const unsigned char *) reference_data, reference_length, - alphabet_length, band_width, MYERS_MODE_HW, &score, &positions, &num_positions, + alphabet_length, -1, MYERS_MODE_HW, &score, &positions, &num_positions, true, &alignment, &alignment_length); if (myers_return_code == MYERS_STATUS_ERROR || num_positions == 0 || alignment_length == 0) { @@ -560,7 +560,7 @@ int MyersNWWrapper(const int8_t *read_data, int64_t read_length, int myers_return_code = myersCalcEditDistance((const unsigned char *) read_data, read_length, (const unsigned char *) reference_data, reference_length, - alphabet_length, band_width, MYERS_MODE_NW, &score, &positions, &num_positions, + alphabet_length, -1, MYERS_MODE_NW, &score, &positions, &num_positions, true, &alignment, &alignment_length); if (myers_return_code == MYERS_STATUS_ERROR || num_positions == 0 || alignment_length == 0) { diff --git a/src/graphmap/anchored_mapping.cc b/src/graphmap/anchored_mapping.cc index 1225a4b..cc01244 100644 --- a/src/graphmap/anchored_mapping.cc +++ b/src/graphmap/anchored_mapping.cc @@ -43,8 +43,7 @@ bool CheckDistanceStep(const Vertices& registry_entries, int64_t index_first, in return false; } - -int GenerateClusters(int64_t min_cluster_length, float min_cluster_coverage, std::vector &lcskpp_indices, +int GenerateClusters(int64_t min_num_anchors_in_cluster, int64_t min_cluster_length, float min_cluster_coverage, std::vector &lcskpp_indices, ScoreRegistry* local_score, MappingData* mapping_data, const std::vector indexes, const SingleSequence* read, const ProgramParameters* parameters, std::vector &ret_clusters, std::vector &ret_filtered_lcskpp_indices, std::vector *ret_cluster_ids) { @@ -124,6 +123,14 @@ int GenerateClusters(int64_t min_cluster_length, float min_cluster_coverage, std int64_t covered_bases = ret_clusters[i]->coverage; float cluster_coverage = ((float) covered_bases) / ((float) cluster_length); + // Filter clusters which are too small. +// if ((min_num_anchors_in_cluster > 0 && min_cluster_length > 0 && min_cluster_coverage > 0) && ret_clusters[i]->lcskpp_indices.size() <= min_num_anchors_in_cluster && cluster_length < min_cluster_length) { + if (ret_clusters[i]->lcskpp_indices.size() <= min_num_anchors_in_cluster && + cluster_length < min_cluster_length && + ret_clusters[i]->lcskpp_indices.size() <= min_num_anchors_in_cluster) { + ret_clusters[i]->lcskpp_indices.clear(); + continue; + } ret_filtered_lcskpp_indices.insert(ret_filtered_lcskpp_indices.end(), ret_clusters[i]->lcskpp_indices.begin(), ret_clusters[i]->lcskpp_indices.end()); // printf ("ret_filtered_lcskpp_indices.size() = %ld\n", ret_filtered_lcskpp_indices.size()); @@ -135,7 +142,7 @@ int GenerateClusters(int64_t min_cluster_length, float min_cluster_coverage, std ret_cluster_ids->insert(ret_cluster_ids->end(), cluster_indices.begin(), cluster_indices.end()); } - LogSystem::GetInstance().Log(VERBOSE_LEVEL_ALL_DEBUG, read->get_sequence_id() == parameters->debug_read, FormatString("[Cluster %ld] cluster_length = %ld, covered_bases = %ld\n", i, cluster_length, covered_bases), "L1-PostProcessRegionWithLCS_"); + LogSystem::GetInstance().Log(VERBOSE_LEVEL_ALL_DEBUG, read->get_sequence_id() == parameters->debug_read, FormatString("[Cluster %ld] cluster_length = %ld, covered_bases = %ld\n", i, cluster_length, covered_bases), std::string(__FUNCTION__)); for (int64_t j=0; jlcskpp_indices.size(); j++) { int32_t qpos_start = 0, rpos_start = 0, qpos_end = 0, rpos_end = 0; GetPositionsFromRegistry2(registry_entries, ret_clusters[i]->lcskpp_indices[j], &qpos_start, &rpos_start, &qpos_end, &rpos_end); @@ -143,6 +150,12 @@ int GenerateClusters(int64_t min_cluster_length, float min_cluster_coverage, std } } + for (int64_t i = (ret_clusters.size() - 1); i >= 0; i--) { + if (ret_clusters[i]->lcskpp_indices.size() == 0) { + ret_clusters.erase(ret_clusters.begin() + i); + } + } + return 0; } @@ -462,10 +475,10 @@ int GraphMap::AnchoredPostProcessRegionWithLCS_(ScoreRegistry* local_score, Mapp std::vector first_filtered_clusters; // FilterAnchors(read, local_score, parameters, lcskpp_indices, parameters->error_rate/2.0f + 0.01f, 300, 2, 50, 2, first_filtered_lcskpp_indices, NULL); // FilterAnchors(read, local_score, parameters, lcskpp_indices, parameters->error_rate/2.0f + 0.01f, std::max(read->get_sequence_length() * 0.10f, 200.0f), 2, 50, 2, first_filtered_lcskpp_indices, NULL); - FilterAnchors(read, local_score, parameters, lcskpp_indices, parameters->error_rate/2.0f + 0.01f, 200.0f, 0, 50, 2, first_filtered_lcskpp_indices, NULL); + FilterAnchors(read, local_score, parameters, lcskpp_indices, parameters->error_rate/2 + 0.01f, 200.0f, 0, 50, 2, first_filtered_lcskpp_indices, NULL); // FilterAnchorBreakpoints(read->get_sequence_length() * MIN_CLUSTER_LENGTH_FACTOR, MIN_CLUSTER_COVERAGE_FACTOR, first_filtered_lcskpp_indices, local_score, mapping_data, indexes, read, parameters, clusters, cluster_indices, &cluster_ids); - GenerateClusters(0.0f, 0.0f, first_filtered_lcskpp_indices, local_score, mapping_data, indexes, read, parameters, clusters, cluster_indices, &cluster_ids); + GenerateClusters(2, 50, 20, first_filtered_lcskpp_indices, local_score, mapping_data, indexes, read, parameters, clusters, cluster_indices, &cluster_ids); // Find the L1 parameters (median line and the confidence intervals). float l_diff = read->get_sequence_length() * parameters->error_rate; diff --git a/src/graphmap/filter_anchors.cc b/src/graphmap/filter_anchors.cc index e72d609..8b9b8f7 100644 --- a/src/graphmap/filter_anchors.cc +++ b/src/graphmap/filter_anchors.cc @@ -6,6 +6,7 @@ */ #include "filter_anchors.h" +#include "log_system/log_system.h" //#define DEBUG_TEST1 @@ -27,6 +28,9 @@ int64_t CalcScore(int32_t qpos, int32_t rpos, int32_t next_qpos, int32_t next_rp *score_dist = -mismatch / 2.0; // Divide by two because not necessarily all bases on the diagonal are mismatches, some might be matches. Just a heuristic. int32_t bandwidth = indel_bandwidth_margin * dist; +// int32_t bandwidth = indel_bandwidth_margin * fwd_length; +// printf (" -> fwd_length = %d, dist = %f, mismatch = %f, band = %f\n", fwd_length, dist, mismatch, band); +// fflush(stdout); int32_t x_dist = abs(qpos - next_qpos); @@ -43,6 +47,8 @@ int64_t CalcScore(int32_t qpos, int32_t rpos, int32_t next_qpos, int32_t next_rp return 2; } else if ((fwd_length > 0 && x_dist > fwd_length) && band > bandwidth) { +// printf ("band = %f, bandwidth = %d\n", band, bandwidth); +// fflush(stdout); // All numbers are too high, break the chain. return 3; } @@ -50,40 +56,6 @@ int64_t CalcScore(int32_t qpos, int32_t rpos, int32_t next_qpos, int32_t next_rp return 4; } -int64_t CalcScore1(int32_t qpos, int32_t rpos, int32_t next_qpos, int32_t next_rpos, double indel_bandwidth_margin, int32_t fwd_length) { -// int32_t min_dist = (next_qpos <= qpos) ? (qpos - next_qpos) : -// (next_qpos >= (qpos + seed_length_q)) ? (next_qpos - qpos - seed_length_q) : -// -1; - int32_t min_dist = (next_qpos <= qpos) ? (qpos - next_qpos) : (next_qpos - qpos); - - if (min_dist < 0 || min_dist > fwd_length) { -#ifdef DEBUG_TEST1 - printf (" --[min_dist not satisfied!] min_dist = %ld, fwd_length = %ld\n", min_dist, fwd_length); - fflush(stdout); -#endif - return 1; - } - - int32_t l1 = rpos - qpos, l2 = next_rpos - next_qpos; - double a = abs(l2 - l1); - double band = a * (sqrt(2.0) / 2.0); - double dist = sqrt((double) (rpos - next_rpos) * (rpos - next_rpos) + (qpos - next_qpos) * (qpos - next_qpos)); - int32_t bandwidth = indel_bandwidth_margin * dist; -#ifdef DEBUG_TEST1 - printf (" bandwidth = %ld, band = %f, dist = %f, a = %f\n", bandwidth, band, dist, a); - fflush(stdout); -#endif - if (band > bandwidth) { -#ifdef DEBUG_TEST1 - printf (" --[bandwidth not satisfied!]\n", bandwidth, band, dist, a); - fflush(stdout); -#endif - return 1; - } // Do I need to check the diagonal, or the actual gap? 'a' is the vertical/horizontal distance, and 'band' is the diagonal between two lines. - - return -a; -} - // Used for packed 128-bit hits, mainly for Owler. void GetPositionsFrom128bit(const std::vector &hits, const std::vector &lcskpp_indices, int64_t lcskpp_id, int32_t seed_len, int32_t *qpos_start, int32_t *rpos_start, int32_t *qpos_end, int32_t *rpos_end) { *qpos_start = get128_qpos(hits[lcskpp_indices[lcskpp_id]]); @@ -109,7 +81,7 @@ void GetPositionsFromRegistry2(const Vertices& registry_entries, int64_t vertex_ *rpos_end = registry_entries.reference_ends[vertex_id]; } -int FilterAnchors(const SingleSequence* seq, ScoreRegistry* local_score, const ProgramParameters *parameters, +int FilterAnchors1(const SingleSequence* seq, ScoreRegistry* local_score, const ProgramParameters *parameters, const std::vector &lcskpp_indices, double indel_bandwidth_margin, int32_t max_dist, int32_t lookahead_dist_factor, int64_t min_covered_bases, int32_t cluster_size_cutoff, std::vector &ret_filtered_lcskpp_indices, std::vector *ret_cluster_ids) { @@ -304,3 +276,159 @@ int FilterAnchors(const SingleSequence* seq, ScoreRegistry* local_score, const P return 0; } + + + + + + +int FilterAnchors(const SingleSequence* read, ScoreRegistry* local_score, const ProgramParameters *parameters, + const std::vector &lcskpp_indices, double indel_bandwidth_margin, int32_t max_dist, int32_t lookahead_dist_factor, int64_t min_covered_bases, int32_t cluster_size_cutoff, + std::vector &ret_filtered_lcskpp_indices, std::vector *ret_cluster_ids) { + + int64_t chain_len_lookahead = max_dist; + const Vertices& registry_entries = local_score->get_registry_entries(); + + std::vector > chains; + std::vector chain_cov_bases; + + { + std::vector new_chain; + new_chain.reserve(lcskpp_indices.size()); + chains.push_back(new_chain); + chain_cov_bases.push_back(0); + } + + LOG_DEBUG_SPEC_NO_HEADER("Starting a new chain: %ld\n", chains.size()); +// fflush(stdout); + for (int64_t i=0; i 0) { + LOG_DEBUG_SPEC_NO_HEADER(" - Breaking. CalcScore returned a value > 0 (ret_val = %ld)\n", ret_val); + break; + } + + if (max_score_id == -1 || (score_gap >= max_score_gap && score_mismatch >= max_score_mismatch)) { + LOG_DEBUG_SPEC_NO_HEADER(" - This will be the new max_score_id.", ret_val); + max_score_gap = score_gap; + max_score_mismatch = score_mismatch; + max_score_id = next_id; + } + next_id += 1; + } + + if (max_score_id == -1) { + std::vector new_chain; + chains.push_back(new_chain); + chains.back().reserve(lcskpp_indices.size() - i); + chain_cov_bases.push_back(0); + + int32_t qpos_start = 0, rpos_start = 0, qpos_end = 0, rpos_end = 0; + GetPositionsFromRegistry(registry_entries, lcskpp_indices, i, &qpos_start, &rpos_start, &qpos_end, &rpos_end); + int32_t cov_bases = std::max((qpos_end - qpos_start), (rpos_end - rpos_start)); + chains.back().push_back(lcskpp_indices[i]); + chain_cov_bases.back() += cov_bases; + + LOG_DEBUG_SPEC_NO_HEADER(" Adding anchor %ld: start = [%d, %d], end = [%d, %d]. (Case where max_score_id == -1)\n", i, qpos_start, rpos_start, qpos_end, rpos_end); + LOG_DEBUG_SPEC_NO_HEADER(" End of chain.\n\n"); +// fflush(stdout); + LOG_DEBUG_SPEC_NO_HEADER("Starting a new chain: %ld\n", chains.size()); +// fflush(stdout); + + } else { + int32_t max_qpos_start = 0, max_rpos_start = 0, max_qpos_end = 0, max_rpos_end = 0; + GetPositionsFromRegistry(registry_entries, lcskpp_indices, max_score_id, &max_qpos_start, &max_rpos_start, &max_qpos_end, &max_rpos_end); + int32_t max_cov_bases = std::max((max_qpos_end - max_qpos_start), (max_rpos_end - max_rpos_start)); + + chains.back().push_back(lcskpp_indices[max_score_id]); + chain_cov_bases.back() += max_cov_bases; + i = max_score_id; + + LOG_DEBUG_SPEC_NO_HEADER(" Adding anchor, max_score_id = %ld, start = [%d, %d], end = [%d, %d]\n", max_score_id, max_qpos_start, max_rpos_start, max_qpos_end, max_rpos_end); +// fflush(stdout); + } + } + + ret_filtered_lcskpp_indices.size(); + if (ret_cluster_ids) { ret_cluster_ids->clear(); } + + // Fill the return indices. + for (int64_t i=0; ireserve(num_filtered_indices); } + + // Fill the return indices. + for (int64_t i=0; ipush_back(chains.size() - i); + } + } + } + +// fflush(stdout); +// exit(1); + + return 0; + +} diff --git a/src/graphmap/graphmap.cc b/src/graphmap/graphmap.cc index 8333878..1e8717a 100644 --- a/src/graphmap/graphmap.cc +++ b/src/graphmap/graphmap.cc @@ -80,7 +80,8 @@ void GraphMap::Run(ProgramParameters& parameters) { // double average_num_kmers = ((double) num_kmers_in_genome) / ((double) num_kmers); // parameters.max_num_hits = (int64_t) ceil(average_num_kmers) * 500; int64_t max_seed_count = 0; - ((IndexSpacedHashFast *) this->indexes_[0])->CalcPercentileHits(0.999, ¶meters.max_num_hits, &max_seed_count); +// ((IndexSpacedHashFast *) this->indexes_[0])->CalcPercentileHits(0.9999, ¶meters.max_num_hits, &max_seed_count); + ((IndexSpacedHashFast *) this->indexes_[0])->CalcPercentileHits(0.9999, ¶meters.max_num_hits, &max_seed_count); LOG_ALL("Automatically setting the maximum number of seed hits to: %ld. Maximum seed occurrence in index: %ld.\n", parameters.max_num_hits, max_seed_count); // LogSystem::GetInstance().VerboseLog(VERBOSE_LEVEL_ALL, true, FormatString("Automatically setting the maximum number of kmer hits: %ld\n", parameters.max_num_hits), "Run"); diff --git a/src/graphmap/region_selection.cc b/src/graphmap/region_selection.cc index 1148360..0068d7a 100644 --- a/src/graphmap/region_selection.cc +++ b/src/graphmap/region_selection.cc @@ -755,20 +755,24 @@ int GraphMap::RegionSelectionNoCopyWithDensehash_(int64_t bin_size, MappingData* int64_t *hits = hit_vector[hits_id]; total_num_hits += hit_counts[hits_id]; + int64_t prev_position_bin = -1, prev_reference_index = -1; + for (int64_t j = 0; j < hit_counts[hits_id]; j++) { int64_t position = hits[j]; int64_t local_position = (int64_t) (((uint64_t) position) & MASK_32_BIT); int64_t reference_index = (int64_t) (((uint64_t) position) >> 32); // (raw_position - reference_starting_pos_[(uint64_t) reference_index]); - if (self_overlap == true && - (reference_index % num_fwd_seqs) == read->get_sequence_id()) { + if (self_overlap == true && (reference_index % num_fwd_seqs) == read->get_sequence_id()) { continue; } - if (reference_index < 0) { - LogSystem::GetInstance().Log(VERBOSE_LEVEL_ALL_DEBUG, read->get_sequence_id() == parameters->debug_read, LogSystem::GetInstance().GenerateErrorMessage(ERR_UNEXPECTED_VALUE, "Offending variable: reference_index. reference_index = %ld, y = %ld, j = %ld / (%ld, %ld)\n", reference_index, local_position, j, 0, hit_counts[hits_id]), "SelectRegionsWithHoughAndCircular"); - continue; - } +// if (reference_index < 0 || reference_index >= num_seqs) { +// LogSystem::GetInstance().Log(VERBOSE_LEVEL_ALL_DEBUG, read->get_sequence_id() == parameters->debug_read, LogSystem::GetInstance().GenerateErrorMessage(ERR_UNEXPECTED_VALUE, "Offending variable: reference_index. reference_index = %ld, y = %ld, j = %ld / (%ld, %ld)\n", reference_index, local_position, j, 0, hit_counts[hits_id]), "SelectRegionsWithHoughAndCircular"); +// printf ("Tu sam 123!\n"); +// fflush(stdout); +// exit(1); +// continue; +// } // Convert the absolute coordinates to local coordinates on the hit reference. int64_t x = i; // Coordinate on the read. @@ -785,8 +789,11 @@ int GraphMap::RegionSelectionNoCopyWithDensehash_(int64_t bin_size, MappingData* // Calculate the index of the bin the position belongs to. int64_t position_bin = floor(((float) l_local) * bin_size_inverse); - if (reference_index >= num_seqs || - position_bin >= max_num_bins[reference_index]) { + if (reference_index == prev_reference_index && position_bin == prev_position_bin) { continue; } + prev_reference_index = reference_index; + prev_position_bin = position_bin; + + if (position_bin >= max_num_bins[reference_index]) { continue; } diff --git a/src/index/index_spaced_hash_fast.cc b/src/index/index_spaced_hash_fast.cc index ff3d1fe..dcda85f 100644 --- a/src/index/index_spaced_hash_fast.cc +++ b/src/index/index_spaced_hash_fast.cc @@ -662,15 +662,15 @@ int IndexSpacedHashFast::DeserializeIndex_(FILE* fp_in) { } -#ifndef RELEASE_VERSION - FILE *fp_debug = fopen (FormatString("temp.kmercounts.%s.csv", shape_index_).c_str(), "w"); - std::vector kmer_counts_indices; - OrderedSortArray(kmer_counts_, num_kmers_, kmer_counts_indices); - for (int64_t i1=(kmer_counts_indices.size()-1); i1>=0; i1--) { - fprintf (fp_debug, "%6X\t%ld\n", kmer_counts_indices[i1], kmer_counts_[kmer_counts_indices[i1]]); - } - fclose(fp_debug); -#endif +//#ifndef RELEASE_VERSION +// FILE *fp_debug = fopen (FormatString("temp.kmercounts.%s.csv", shape_index_).c_str(), "w"); +// std::vector kmer_counts_indices; +// OrderedSortArray(kmer_counts_, num_kmers_, kmer_counts_indices); +// for (int64_t i1=(kmer_counts_indices.size()-1); i1>=0; i1--) { +// fprintf (fp_debug, "%6X\t%ld\n", kmer_counts_indices[i1], kmer_counts_[kmer_counts_indices[i1]]); +// } +// fclose(fp_debug); +//#endif return 0; }