C++ KmerCounter::clear_kmer方法代码示例

Path_n_count_pair IRKE::inchworm(KmerCounter &kcounter,
                                 char direction,
                                 kmer_int_type_t kmer,
                                 Kmer_visitor &visitor,
                                 float min_connectivity)
{

    // cout << "inchworm" << endl;

    Path_n_count_pair entire_path;
    entire_path.second = 0; // init cumulative path coverage

    unsigned int inchworm_round = 0;

    unsigned long num_total_kmers = kcounter.size();

    Kmer_visitor eliminator(kcounter.get_kmer_length(), DOUBLE_STRANDED_MODE);

    while (true) {


        if (IRKE_COMMON::__DEVEL_rand_fracture) {

            // terminate extension with probability of __DEVEL_rand_fracture_prob

            float prob_to_fracture = rand() / (float) RAND_MAX;
            //cerr << "prob: " << prob_to_fracture << endl;

            if (prob_to_fracture <= IRKE_COMMON::__DEVEL_rand_fracture_prob) {

                // cerr << "Fracturing at iworm round: " << inchworm_round << " given P: " << prob_to_fracture << endl;

                return (entire_path);
            }
        }

        inchworm_round++;
        eliminator.clear();

        if (inchworm_round > num_total_kmers) {
            throw (string("Error, inchworm rounds have exceeded the number of possible seed kmers"));
        }

        if (IRKE_COMMON::MONITOR >= 3) {
            cerr << endl << "Inchworm round(" << string(1, direction) << "): " << inchworm_round << " searching kmer: "
                << kmer << endl;
            string kmer_str = kcounter.get_kmer_string(kmer);
            cerr << kcounter.describe_kmer(kmer_str) << endl;
        }

        visitor.erase(kmer); // seed kmer must be not visited already.

        Kmer_Occurence_Pair kmer_pair(kmer, kcounter.get_kmer_count(kmer));
        Path_n_count_pair best_path = inchworm_step(kcounter,
                                                    direction,
                                                    kmer_pair,
                                                    visitor,
                                                    eliminator,
                                                    inchworm_round,
                                                    0,
                                                    min_connectivity,
                                                    MAX_RECURSION);

        vector<kmer_int_type_t> &kmer_list = best_path.first;
        unsigned int num_kmers = kmer_list.size();

        if ((IRKE_COMMON::__DEVEL_zero_kmer_on_use && num_kmers >= 1) || best_path.second > 0) {
            // append info to entire path in reverse order, so starts just after seed kmer

            int first_index = num_kmers - 1;
            int last_index = 0;
            if (CRAWL) {
                last_index = first_index - CRAWL_LENGTH + 1;
                if (last_index < 0) {
                    last_index = 0;
                }
            }

            for (int i = first_index; i >= last_index; i--) {
                kmer_int_type_t kmer_extend = kmer_list[i];
                entire_path.first.push_back(kmer_extend);
                visitor.add(kmer_extend);
                //entire_path.second += kcounter.get_kmer_count(kmer_extend);

                // selected here, zero out:


                if (IRKE_COMMON::__DEVEL_zero_kmer_on_use) {
                    kcounter.clear_kmer(kmer_extend);
                }

            }

            kmer = entire_path.first[entire_path.first.size() - 1];

            entire_path.second += best_path.second;

        }
        else {
            // no extension possible
//.........这里部分代码省略.........

//.........这里部分代码省略.........
            for (unsigned int i = 0; i < reverse_path.size(); i++) {
                cerr  << "\tReverse path kmer: " << kcounter.get_kmer_string(reverse_path[i]) << endl; 
            }
        }
        
		
		unsigned int total_counts = selected_path_n_pair_forward.second + selected_path_n_pair_reverse.second + kcounter.get_kmer_count(kmer); 
		
		vector<kmer_int_type_t>& reverse_path = selected_path_n_pair_reverse.first;
		
		vector<kmer_int_type_t> joined_path = _join_forward_n_reverse_paths(reverse_path, kmer, forward_path);
		
		// report sequence reconstructed from path.
		
		vector<unsigned int> assembly_base_coverage;
		string sequence = reconstruct_path_sequence(kcounter, joined_path, assembly_base_coverage);
		
		unsigned int avg_cov =  static_cast<unsigned int> ( (float)total_counts/(sequence.length()-kcounter.get_kmer_length() +1) + 0.5);
		
		/*
		  cout << "Inchworm-reconstructed sequence, length: " << sequence.length() 
		  << ", avgCov: " << avg_cov
		  << " " << sequence << endl;
		*/
		
		
		
		if (sequence.length() >= MIN_ASSEMBLY_LENGTH && avg_cov >= MIN_ASSEMBLY_COVERAGE) {
			
			INCHWORM_ASSEMBLY_COUNTER++;
			
			stringstream headerstream;
			
			
			headerstream << ">a" << INCHWORM_ASSEMBLY_COUNTER << ";" << avg_cov 
						 << " K: " << kmer_length
						 << " length: " << sequence.length();
			
			string header = headerstream.str();
			
            sequence = add_fasta_seq_line_breaks(sequence, 60);
            
			cout << header << endl << sequence << endl;
			
			if (WRITE_COVERAGE) {
				
				coverage_writer << header << endl;
				
				for (unsigned int i = 0; i < assembly_base_coverage.size(); i++) {
					coverage_writer << assembly_base_coverage[i];
					if ( (i+1) % 30 == 0) {
						coverage_writer << endl;
					}
					else {
						coverage_writer << " ";
					}
				}
				coverage_writer << endl;
			}
			
		}
		
		// remove path
		for (unsigned int i = 0; i < joined_path.size(); i++) {
			
			kmer_int_type_t kmer = joined_path[i];
			
			/*
			  if (DEBUG) {
			  cout << "\tpruning kmer: " << kmer << endl;
			  }
			*/
			
				kcounter.clear_kmer(kmer);
		}
		
		/*
		  if (DEBUG) {
		  cout << "done pruning kmers." << endl;
		  }
		*/
		
	}
	
	if (IRKE_COMMON::MONITOR) {
		cerr << endl;
	}
	
	if (WRITE_COVERAGE) {
		coverage_writer.close();
	}

    
    // drop sorted kmer list as part of cleanup
    clear_sorted_kmers_list();
    
	
	return; // end of runIRKE
	
}

//.........这里部分代码省略.........
        // report sequence reconstructed from path.

        vector<unsigned int> assembly_base_coverage;

        string sequence = reconstruct_path_sequence(kcounter, joined_path, assembly_base_coverage);

        unsigned int avg_cov =
            static_cast<unsigned int> ((float) total_counts / (sequence.length() - kcounter.get_kmer_length() + 1)
                + 0.5);

        /*
          cout << "Inchworm-reconstructed sequence, length: " << sequence.length()
          << ", avgCov: " << avg_cov
          << " " << sequence << endl;
        */

        size_t contig_length = sequence.length();

        if (contig_length >= MIN_ASSEMBLY_LENGTH && avg_cov >= MIN_ASSEMBLY_COVERAGE) {


            *(tmpfiles[myTid].fh) << total_counts << endl
                << avg_cov << endl
                << kmer_count << endl
                << sequence << endl;

        }

        // remove path

        if (IRKE_COMMON::__DEVEL_zero_kmer_on_use) {

            // dont forget the seed. The forward/reverse path kmers already cleared.
            kcounter.clear_kmer(kmer);

        } else {

            for (unsigned int i = 0; i < joined_path.size(); i++) {

                kmer_int_type_t kmer = joined_path[i];
                kcounter.clear_kmer(kmer);
            }

        }

    }

    if (IRKE_COMMON::MONITOR) {
        cerr << endl;
    }

    if (WRITE_COVERAGE) {
        coverage_writer.close();
    }

    // drop sorted kmer list as part of cleanup
    clear_sorted_kmers_list();


    //------------------------------------------------------------------------------      
    // examine the contigs generated by the individual threads, remove redundancies.
    //------------------------------------------------------------------------------

    map<unsigned long long, bool> seen_contig_already;

    for (unsigned int i = 0; i < tmpfiles.size(); i++) {