Commit e9a08938 authored by Li's avatar Li
Browse files

Preliminary implementaion of minimizer cache

parent ecf36897

src/chromap.cc

+133 −41
Original line number Diff line number Diff line
@@ -12,8 +12,16 @@ 


#include "cxxopts.hpp"

#include "ksw.h"

#include "mmcache.hpp"



namespace chromap {

struct _mm_history

{

	std::vector<std::pair<uint64_t, uint64_t> > minimizers ;

	std::vector<struct _candidate> positive_candidates ;

	std::vector<struct _candidate> negative_candidates ;

} ;



template <typename MappingRecord>

void Chromap<MappingRecord>::TrimAdapterForPairedEndRead(uint32_t pair_index, SequenceBatch *read_batch1, SequenceBatch *read_batch2) {

  const char *read1 = read_batch1->GetSequenceAt(pair_index);
@@ -370,10 +378,10 @@ void Chromap<MappingRecord>::MapPairedEndReads() { 
  positive_hits2.reserve(max_seed_frequencies_[0]);

  negative_hits1.reserve(max_seed_frequencies_[0]);

  negative_hits2.reserve(max_seed_frequencies_[0]);

  std::vector<uint64_t> positive_candidates1;

  std::vector<uint64_t> positive_candidates2;

  std::vector<uint64_t> negative_candidates1;

  std::vector<uint64_t> negative_candidates2;

  std::vector<struct _candidate> positive_candidates1;

  std::vector<struct _candidate> positive_candidates2;

  std::vector<struct _candidate> negative_candidates1;

  std::vector<struct _candidate> negative_candidates2;

  positive_candidates1.reserve(max_seed_frequencies_[0]);

  positive_candidates2.reserve(max_seed_frequencies_[0]);

  negative_candidates1.reserve(max_seed_frequencies_[0]);
@@ -434,15 +442,20 @@ void Chromap<MappingRecord>::MapPairedEndReads() { 
        index.GenerateCandidates(error_threshold_, minimizers2, &positive_hits2, &negative_hits2, &positive_candidates2, &negative_candidates2);

        uint32_t current_num_candidates2 = positive_candidates2.size() + negative_candidates2.size();

        if (current_num_candidates1 > 0 && current_num_candidates2 > 0) {

          positive_candidates1.swap(positive_hits1);

          /*positive_candidates1.swap(positive_hits1);

          negative_candidates1.swap(negative_hits1);

          positive_candidates2.swap(positive_hits2);

          negative_candidates2.swap(negative_hits2);

          negative_candidates2.swap(negative_hits2);*/

	  std::vector<struct _candidate> raw_pos1(positive_candidates1) ;

	  std::vector<struct _candidate> raw_neg1(negative_candidates1) ;

	  std::vector<struct _candidate> raw_pos2(positive_candidates2) ;

	  std::vector<struct _candidate> raw_neg2(negative_candidates2) ;

	  

          positive_candidates1.clear();

          positive_candidates2.clear();

          negative_candidates1.clear();

          negative_candidates2.clear();

          ReduceCandidatesForPairedEndRead(positive_hits1, negative_hits1, positive_hits2, negative_hits2, &positive_candidates1, &negative_candidates1, &positive_candidates2, &negative_candidates2);

          ReduceCandidatesForPairedEndRead(raw_pos1, raw_neg1, raw_pos2, raw_neg2, &positive_candidates1, &negative_candidates1, &positive_candidates2, &negative_candidates2);

          thread_num_candidates += positive_candidates1.size() + positive_candidates2.size() + negative_candidates1.size() + negative_candidates2.size();

          positive_mappings1.clear();

          positive_mappings2.clear();
@@ -452,9 +465,9 @@ void Chromap<MappingRecord>::MapPairedEndReads() { 
          int num_best_mappings1, num_second_best_mappings1;

          int min_num_errors2, second_min_num_errors2;

          int num_best_mappings2, num_second_best_mappings2;

          VerifyCandidates(read_batch1, pair_index, reference, positive_candidates1, negative_candidates1, &positive_mappings1, &negative_mappings1, &min_num_errors1, &num_best_mappings1, &second_min_num_errors1, &num_second_best_mappings1);

          VerifyCandidates(read_batch1, pair_index, reference, minimizers1, positive_candidates1, negative_candidates1, &positive_mappings1, &negative_mappings1, &min_num_errors1, &num_best_mappings1, &second_min_num_errors1, &num_second_best_mappings1);

          uint32_t current_num_mappings1 = positive_mappings1.size() + negative_mappings1.size();

          VerifyCandidates(read_batch2, pair_index, reference, positive_candidates2, negative_candidates2, &positive_mappings2, &negative_mappings2, &min_num_errors2, &num_best_mappings2, &second_min_num_errors2, &num_second_best_mappings2);

          VerifyCandidates(read_batch2, pair_index, reference, minimizers2, positive_candidates2, negative_candidates2, &positive_mappings2, &negative_mappings2, &min_num_errors2, &num_best_mappings2, &second_min_num_errors2, &num_second_best_mappings2);

          uint32_t current_num_mappings2 = positive_mappings2.size() + negative_mappings2.size();

          if (current_num_mappings1 > 0 && current_num_mappings2 > 0) {

            int min_sum_errors, second_min_sum_errors;
@@ -559,22 +572,22 @@ void Chromap<MappingRecord>::OutputMappings(uint32_t num_reference_sequences, co 
}



template <typename MappingRecord>

void Chromap<MappingRecord>::ReduceCandidatesForPairedEndReadOnOneDirection(const std::vector<uint64_t> &candidates1, const std::vector<uint64_t> &candidates2, std::vector<uint64_t> *filtered_candidates1, std::vector<uint64_t> *filtered_candidates2) {

void Chromap<MappingRecord>::ReduceCandidatesForPairedEndReadOnOneDirection(const std::vector<struct _candidate> &candidates1, const std::vector<struct _candidate> &candidates2, std::vector<struct _candidate> *filtered_candidates1, std::vector<struct _candidate> *filtered_candidates2) {

  uint32_t i1 = 0;

  uint32_t i2 = 0;

  uint32_t mapping_positions_distance = max_insert_size_;

  uint32_t previous_end_i2 = i2;

  while (i1 < candidates1.size() && i2 < candidates2.size()) {

    if (candidates1[i1] > candidates2[i2] + mapping_positions_distance) {

    if (candidates1[i1].refPos > candidates2[i2].refPos + mapping_positions_distance) {

      ++i2;

    } else if (candidates2[i2] > candidates1[i1] + mapping_positions_distance) {

    } else if (candidates2[i2].refPos > candidates1[i1].refPos + mapping_positions_distance) {

      ++i1;

    } else {

      // ok, find a pair, we store current ni2 somewhere and keep looking until we go out of the range, 

      // then we go back and then move to next pi1 and keep doing the similar thing. 

      filtered_candidates1->emplace_back(candidates1[i1]);

      uint32_t current_i2 = i2;

      while (current_i2 < candidates2.size() && candidates2[current_i2] <= candidates1[i1] + mapping_positions_distance) {

      while (current_i2 < candidates2.size() && candidates2[current_i2].refPos <= candidates1[i1].refPos + mapping_positions_distance) {

        if (current_i2 >= previous_end_i2) {

          filtered_candidates2->emplace_back(candidates2[current_i2]);

        }
@@ -587,7 +600,7 @@ void Chromap<MappingRecord>::ReduceCandidatesForPairedEndReadOnOneDirection(cons 
}



template <typename MappingRecord>

void Chromap<MappingRecord>::ReduceCandidatesForPairedEndRead(const std::vector<uint64_t> &positive_candidates1, const std::vector<uint64_t> &negative_candidates1, const std::vector<uint64_t> &positive_candidates2, const std::vector<uint64_t> &negative_candidates2, std::vector<uint64_t> *filtered_positive_candidates1, std::vector<uint64_t> *filtered_negative_candidates1, std::vector<uint64_t> *filtered_positive_candidates2, std::vector<uint64_t> *filtered_negative_candidates2) {

void Chromap<MappingRecord>::ReduceCandidatesForPairedEndRead(const std::vector<struct _candidate> &positive_candidates1, const std::vector<struct _candidate> &negative_candidates1, const std::vector<struct _candidate> &positive_candidates2, const std::vector<struct _candidate> &negative_candidates2, std::vector<struct _candidate> *filtered_positive_candidates1, std::vector<struct _candidate> *filtered_negative_candidates1, std::vector<struct _candidate> *filtered_positive_candidates2, std::vector<struct _candidate> *filtered_negative_candidates2) {

  ReduceCandidatesForPairedEndReadOnOneDirection(positive_candidates1, negative_candidates2, filtered_positive_candidates1, filtered_negative_candidates2);

  ReduceCandidatesForPairedEndReadOnOneDirection(negative_candidates1, positive_candidates2, filtered_negative_candidates1, filtered_positive_candidates2);

}
@@ -871,6 +884,8 @@ void Chromap<MappingRecord>::MapSingleEndReads() { 
  SequenceBatch read_batch_for_loading(read_batch_size_);

  SequenceBatch barcode_batch(read_batch_size_);

  SequenceBatch barcode_batch_for_loading(read_batch_size_);

  mm_cache mm_to_candidates_cache(1000000) ;

  struct _mm_history *mm_history = new struct _mm_history[read_batch_size_];

  read_batch_for_loading.InitializeLoading(read_file1_path_);

  if (!is_bulk_data_) {

    barcode_batch_for_loading.InitializeLoading(barcode_file_path_);
@@ -911,15 +926,15 @@ void Chromap<MappingRecord>::MapSingleEndReads() { 
  static uint64_t thread_num_mapped_reads = 0; 

  static uint64_t thread_num_uniquely_mapped_reads = 0; 

#pragma omp threadprivate(thread_num_candidates, thread_num_mappings, thread_num_mapped_reads, thread_num_uniquely_mapped_reads)

#pragma omp parallel default(none) shared(reference, index, read_batch, barcode_batch, read_batch_for_loading, barcode_batch_for_loading, std::cerr, num_loaded_reads_for_loading, num_loaded_reads, num_reference_sequences, mappings_on_diff_ref_seqs_for_diff_threads, mappings_on_diff_ref_seqs_for_diff_threads_for_saving) num_threads(num_threads_) reduction(+:num_candidates_, num_mappings_, num_mapped_reads_, num_uniquely_mapped_reads_)

#pragma omp parallel default(none) shared(reference, index, read_batch, barcode_batch, read_batch_for_loading, barcode_batch_for_loading, std::cerr, num_loaded_reads_for_loading, num_loaded_reads, num_reference_sequences, mappings_on_diff_ref_seqs_for_diff_threads, mappings_on_diff_ref_seqs_for_diff_threads_for_saving, mm_to_candidates_cache, mm_history) num_threads(num_threads_) reduction(+:num_candidates_, num_mappings_, num_mapped_reads_, num_uniquely_mapped_reads_)

  {

  std::vector<std::pair<uint64_t, uint64_t> > minimizers;

  std::vector<uint64_t> positive_hits;

  std::vector<uint64_t> negative_hits;

  positive_hits.reserve(max_seed_frequencies_[0]);

  negative_hits.reserve(max_seed_frequencies_[0]);

  std::vector<uint64_t> positive_candidates;

  std::vector<uint64_t> negative_candidates;

  std::vector<struct _candidate> positive_candidates;

  std::vector<struct _candidate> negative_candidates;

  positive_candidates.reserve(max_seed_frequencies_[0]);

  negative_candidates.reserve(max_seed_frequencies_[0]);

  std::vector<std::pair<int, uint64_t> > positive_mappings;
@@ -949,7 +964,26 @@ void Chromap<MappingRecord>::MapSingleEndReads() { 
        negative_hits.clear();

        positive_candidates.clear();

        negative_candidates.clear();

        index.GenerateCandidates(error_threshold_, minimizers, &positive_hits, &negative_hits, &positive_candidates, &negative_candidates);

	if (mm_to_candidates_cache.Query(minimizers, positive_candidates, negative_candidates, 

		read_batch.GetSequenceLengthAt(read_index) ) == -1)

		index.GenerateCandidates(error_threshold_, minimizers, &positive_hits, &negative_hits, 

					&positive_candidates, &negative_candidates);

	/*else

	{

		printf("successful cache load.%s\n", read_batch.GetSequenceNameAt(read_index)) ;

	}*/

	/*for (uint32_t i = 0 ; i < positive_candidates.size() ; ++i)

	{

		printf("LI_DEBUG +: %d\n", int(positive_candidates[i].refPos)) ;

	}

	for (uint32_t i = 0 ; i < negative_candidates.size() ; ++i)

	{

		printf("LI_DEBUG -: %d\n", int(negative_candidates[i].refPos)) ;

	}*/

	mm_history[read_index].minimizers = minimizers ;

	mm_history[read_index].positive_candidates = positive_candidates ;

	mm_history[read_index].negative_candidates = negative_candidates ;



        uint32_t current_num_candidates = positive_candidates.size() + negative_candidates.size(); 

        //std::cerr << "Generated candidates!\n";

        if (current_num_candidates > 0) {
@@ -958,7 +992,7 @@ void Chromap<MappingRecord>::MapSingleEndReads() { 
          negative_mappings.clear();

          int min_num_errors, second_min_num_errors;

          int num_best_mappings, num_second_best_mappings;

          VerifyCandidates(read_batch, read_index, reference, positive_candidates, negative_candidates, &positive_mappings, &negative_mappings, &min_num_errors, &num_best_mappings, &second_min_num_errors, &num_second_best_mappings);

          VerifyCandidates(read_batch, read_index, reference, minimizers, positive_candidates, negative_candidates, &positive_mappings, &negative_mappings, &min_num_errors, &num_best_mappings, &second_min_num_errors, &num_second_best_mappings);

          uint32_t current_num_mappings = positive_mappings.size() + negative_mappings.size();

          //std::cerr << "Verified candidates!\n";

          if (current_num_mappings > 0) {
@@ -974,6 +1008,10 @@ void Chromap<MappingRecord>::MapSingleEndReads() { 
        }

      }

    }

    for (uint32_t read_index = 0; read_index < num_loaded_reads; ++read_index) {

    	mm_to_candidates_cache.Update(mm_history[read_index].minimizers, mm_history[read_index].positive_candidates,

				mm_history[read_index].negative_candidates) ;

    }

#pragma omp taskwait

    num_loaded_reads = num_loaded_reads_for_loading;

    read_batch_for_loading.SwapSequenceBatch(read_batch);
@@ -993,6 +1031,7 @@ void Chromap<MappingRecord>::MapSingleEndReads() { 
    num_uniquely_mapped_reads_ += thread_num_uniquely_mapped_reads;

  } // end of updating shared mapping stats

  } // end of openmp parallel region

  delete[] mm_history ;

  read_batch_for_loading.FinalizeLoading();

  if (!is_bulk_data_) {

    barcode_batch_for_loading.FinalizeLoading();
@@ -1417,23 +1456,26 @@ void Chromap<MappingRecord>::AllocateMultiMappings(uint32_t num_reference_sequen 
}



template <typename MappingRecord>

void Chromap<MappingRecord>::VerifyCandidatesOnOneDirectionUsingSIMD(Direction candidate_direction, const SequenceBatch &read_batch, uint32_t read_index, const SequenceBatch &reference, const std::vector<uint64_t> &candidates, std::vector<std::pair<int, uint64_t> > *mappings, int *min_num_errors, int *num_best_mappings, int *second_min_num_errors, int *num_second_best_mappings) {

void Chromap<MappingRecord>::VerifyCandidatesOnOneDirectionUsingSIMD(Direction candidate_direction, const SequenceBatch &read_batch, uint32_t read_index, const SequenceBatch &reference, const std::vector<struct _candidate> &candidates, std::vector<std::pair<int, uint64_t> > *mappings, int *min_num_errors, int *num_best_mappings, int *second_min_num_errors, int *num_second_best_mappings) {

  const char *read = read_batch.GetSequenceAt(read_index);

  uint32_t read_length = read_batch.GetSequenceLengthAt(read_index);

  const std::string &negative_read = read_batch.GetNegativeSequenceAt(read_index); 



  size_t num_candidates = candidates.size();

  uint64_t valid_candidates[NUM_VPU_LANES_];

  struct _candidate valid_candidates[NUM_VPU_LANES_];

  const char *valid_candidate_starts[NUM_VPU_LANES_];

  uint32_t valid_candidate_index = 0;

  size_t candidate_index = 0;

  uint32_t mmCntThreshold = 0 ;

  //int16_t mapping_edit_distances_16[NUM_VPU_LANES_];

  //int16_t mapping_end_positions_16[NUM_VPU_LANES_]; 

  //int32_t mapping_edit_distances_32[NUM_VPU_LANES_];

  //int32_t mapping_end_positions_32[NUM_VPU_LANES_]; 

  while (candidate_index < num_candidates) {

    uint32_t rid = candidates[candidate_index] >> 32;

    uint32_t position = candidates[candidate_index];

    if (candidates[candidate_index].mmCnt < mmCntThreshold)

    	break ;

    uint32_t rid = candidates[candidate_index].refPos >> 32;

    uint32_t position = candidates[candidate_index].refPos ;

    if (candidate_direction == kNegative) {

      position = position - read_length + 1;

    }
@@ -1471,9 +1513,9 @@ void Chromap<MappingRecord>::VerifyCandidatesOnOneDirectionUsingSIMD(Direction c 
              (*num_second_best_mappings)++;

            }

            if (candidate_direction == kPositive) {

              mappings->emplace_back((uint8_t)mapping_edit_distances[mi], valid_candidates[mi] - error_threshold_ + mapping_end_positions[mi]);

              mappings->emplace_back((uint8_t)mapping_edit_distances[mi], valid_candidates[mi].refPos - error_threshold_ + mapping_end_positions[mi]);

            } else {

              mappings->emplace_back((uint8_t)mapping_edit_distances[mi], valid_candidates[mi] - read_length + 1 - error_threshold_ + mapping_end_positions[mi]); 

              mappings->emplace_back((uint8_t)mapping_edit_distances[mi], valid_candidates[mi].refPos - read_length + 1 - error_threshold_ + mapping_end_positions[mi]); 

            }

          }

        }
@@ -1501,21 +1543,28 @@ void Chromap<MappingRecord>::VerifyCandidatesOnOneDirectionUsingSIMD(Direction c 
              (*num_second_best_mappings)++;

            }

            if (candidate_direction == kPositive) {

              mappings->emplace_back((uint8_t)mapping_edit_distances[mi], valid_candidates[mi] - error_threshold_ + mapping_end_positions[mi]);

              mappings->emplace_back((uint8_t)mapping_edit_distances[mi], valid_candidates[mi].refPos - error_threshold_ + mapping_end_positions[mi]);

            } else {

              mappings->emplace_back((uint8_t)mapping_edit_distances[mi], valid_candidates[mi] - read_length + 1 - error_threshold_ + mapping_end_positions[mi]); 

              mappings->emplace_back((uint8_t)mapping_edit_distances[mi], valid_candidates[mi].refPos - read_length + 1 - error_threshold_ + mapping_end_positions[mi]); 

            }

          } else {

	  	mmCntThreshold = valid_candidates[mi].mmCnt ;

	  }

        }

      }

      valid_candidate_index = 0;



      // Check whether we should stop early. Assuming the candidates are sorted 

      if ( GetMAPQ( 0, 0, read_length + error_threshold_, *min_num_errors, *num_best_mappings, 

      	*second_min_num_errors, *second_min_num_errors ) < 30 )

      	break ;

    }

    ++candidate_index;

  }



  for (uint32_t ci = 0; ci < valid_candidate_index; ++ci) {

    uint32_t rid = valid_candidates[ci] >> 32;

    uint32_t position = valid_candidates[ci];

    uint32_t rid = valid_candidates[ci].refPos >> 32;

    uint32_t position = valid_candidates[ci].refPos ;

    if (candidate_direction == kNegative) {

      position = position - read_length + 1;

    }
@@ -1541,23 +1590,23 @@ void Chromap<MappingRecord>::VerifyCandidatesOnOneDirectionUsingSIMD(Direction c 
        (*num_second_best_mappings)++;

      }

      if (candidate_direction == kPositive) {

        mappings->emplace_back(num_errors, valid_candidates[ci] - error_threshold_ + mapping_end_position);

        mappings->emplace_back(num_errors, valid_candidates[ci].refPos - error_threshold_ + mapping_end_position);

      } else {

        mappings->emplace_back(num_errors, valid_candidates[ci] - read_length + 1 - error_threshold_ + mapping_end_position); 

        mappings->emplace_back(num_errors, valid_candidates[ci].refPos - read_length + 1 - error_threshold_ + mapping_end_position); 

      }

    }

  }

}



template <typename MappingRecord>

void Chromap<MappingRecord>::VerifyCandidatesOnOneDirection(Direction candidate_direction, const SequenceBatch &read_batch, uint32_t read_index, const SequenceBatch &reference, const std::vector<uint64_t> &candidates, std::vector<std::pair<int, uint64_t> > *mappings, int *min_num_errors, int *num_best_mappings, int *second_min_num_errors, int *num_second_best_mappings) {

void Chromap<MappingRecord>::VerifyCandidatesOnOneDirection(Direction candidate_direction, const SequenceBatch &read_batch, uint32_t read_index, const SequenceBatch &reference, const std::vector<struct _candidate> &candidates, std::vector<std::pair<int, uint64_t> > *mappings, int *min_num_errors, int *num_best_mappings, int *second_min_num_errors, int *num_second_best_mappings) {

  const char *read = read_batch.GetSequenceAt(read_index);

  uint32_t read_length = read_batch.GetSequenceLengthAt(read_index);

  const std::string &negative_read = read_batch.GetNegativeSequenceAt(read_index); 



  for (uint32_t ci = 0; ci < candidates.size(); ++ci) {

    uint32_t rid = candidates[ci] >> 32;

    uint32_t position = candidates[ci];

    uint32_t rid = candidates[ci].refPos >> 32;

    uint32_t position = candidates[ci].refPos ;

    if (candidate_direction == kNegative) {

      position = position - read_length + 1;

    }
@@ -1583,29 +1632,72 @@ void Chromap<MappingRecord>::VerifyCandidatesOnOneDirection(Direction candidate_ 
        (*num_second_best_mappings)++;

      }

      if (candidate_direction == kPositive) {

        mappings->emplace_back(num_errors, candidates[ci] - error_threshold_ + mapping_end_position);

        mappings->emplace_back(num_errors, candidates[ci].refPos - error_threshold_ + mapping_end_position);

      } else {

        mappings->emplace_back(num_errors, candidates[ci] - read_length + 1 - error_threshold_ + mapping_end_position); 

        mappings->emplace_back(num_errors, candidates[ci].refPos - read_length + 1 - error_threshold_ + mapping_end_position); 

      }

    }

  }

}



template <typename MappingRecord>

void Chromap<MappingRecord>::VerifyCandidates(const SequenceBatch &read_batch, uint32_t read_index, const SequenceBatch &reference, const std::vector<uint64_t> &positive_candidates, const std::vector<uint64_t> &negative_candidates, std::vector<std::pair<int, uint64_t> > *positive_mappings, std::vector<std::pair<int, uint64_t> > *negative_mappings, int *min_num_errors, int *num_best_mappings, int *second_min_num_errors, int *num_second_best_mappings) {

void Chromap<MappingRecord>::VerifyCandidates(const SequenceBatch &read_batch, uint32_t read_index, const SequenceBatch &reference, const std::vector<std::pair<uint64_t, uint64_t> > &minimizers, const std::vector<struct _candidate> &positive_candidates, const std::vector<struct _candidate> &negative_candidates, std::vector<std::pair<int, uint64_t> > *positive_mappings, std::vector<std::pair<int, uint64_t> > *negative_mappings, int *min_num_errors, int *num_best_mappings, int *second_min_num_errors, int *num_second_best_mappings) {

  *min_num_errors = error_threshold_ + 1;

  *num_best_mappings = 0;

  *second_min_num_errors = error_threshold_ + 1;

  *num_second_best_mappings = 0;

  

  // Directly obtain the mapping in ideal case.

  uint32_t i ;

  int maxCnt = 0 ;

  int maxTag = 0 ;

  //printf("LI_DEBUG: %u %u\n", positive_candidates.size() + negative_candidates.size(), minimizers.size()) ;

  for ( i = 0 ; i < positive_candidates.size() ; ++i ) {

  	//printf("+ %u %u\n", i, positive_candidates[i].mmCnt) ;

  	if (positive_candidates[i].mmCnt == minimizers.size()) {

		maxTag = i << 1 ;

		++maxCnt ;

	}

  }

  for ( i = 0 ; i < negative_candidates.size() ; ++i ) {

  	//printf("- %u %u\n", i, negative_candidates[i].mmCnt) ;

  	if (negative_candidates[i].mmCnt == minimizers.size()) {

		maxTag = (i << 1)|1 ;

		++maxCnt ;

	}

  }

  

  if (maxCnt == 1) {

      Direction candidate_direction = (maxTag & 1) ? kNegative : kPositive ;

      uint32_t ci = maxTag >> 1 ;

      *num_best_mappings = 1 ;

      *num_second_best_mappings = 0 ;

      *min_num_errors = 0 ;

      if (candidate_direction == kPositive) {

        uint32_t read_length = read_batch.GetSequenceLengthAt(read_index);

        positive_mappings->emplace_back(0, positive_candidates[ci].refPos + read_length - 1);

      } else {

        negative_mappings->emplace_back(0, negative_candidates[ci].refPos); 

      }

      //printf("Saved %d\n", positive_candidates.size() + negative_candidates.size() ) ;

      return ;

  }

  //printf("Notsaved %d\n", positive_candidates.size() + negative_candidates.size()) ;



  // Use more sophicated approach to obtain the mapping

  if (positive_candidates.size() < (size_t)NUM_VPU_LANES_) {

    VerifyCandidatesOnOneDirection(kPositive, read_batch, read_index, reference, positive_candidates, positive_mappings, min_num_errors, num_best_mappings, second_min_num_errors, num_second_best_mappings);

  } else {

    VerifyCandidatesOnOneDirectionUsingSIMD(kPositive, read_batch, read_index, reference, positive_candidates, positive_mappings, min_num_errors, num_best_mappings, second_min_num_errors, num_second_best_mappings);

    std::vector<struct _candidate> sorted_candidates(positive_candidates) ; 

    std::sort(sorted_candidates.begin(), sorted_candidates.end()) ;

    VerifyCandidatesOnOneDirectionUsingSIMD(kPositive, read_batch, read_index, reference, sorted_candidates, positive_mappings, min_num_errors, num_best_mappings, second_min_num_errors, num_second_best_mappings);

  }

  if (negative_candidates.size() < (size_t)NUM_VPU_LANES_) {

    VerifyCandidatesOnOneDirection(kNegative, read_batch, read_index, reference, negative_candidates, negative_mappings, min_num_errors, num_best_mappings, second_min_num_errors, num_second_best_mappings);

  } else {

    VerifyCandidatesOnOneDirectionUsingSIMD(kNegative, read_batch, read_index, reference, negative_candidates, negative_mappings, min_num_errors, num_best_mappings, second_min_num_errors, num_second_best_mappings);

    std::vector<struct _candidate> sorted_candidates(negative_candidates) ; 

    std::sort(sorted_candidates.begin(), sorted_candidates.end()) ;

    VerifyCandidatesOnOneDirectionUsingSIMD(kNegative, read_batch, read_index, reference, sorted_candidates, negative_mappings, min_num_errors, num_best_mappings, second_min_num_errors, num_second_best_mappings);

  }

}

src/chromap.h

+6 −10

File changed.

Preview size limit exceeded, changes collapsed.

src/index.cc

+19 −4

File changed.

Preview size limit exceeded, changes collapsed.

src/index.h

+22 −2

File changed.

Preview size limit exceeded, changes collapsed.

src/mmcache.hpp

0 → 100644
+161 −0

File added.

Preview size limit exceeded, changes collapsed.

CRA Git | Maintained and supported by SUSTech CRA and CCSE