Introduce FeatureBarcodeMatrix. (b98c2a2d) · Commits · github_fork / Chromap

Makefile

+1 −1

Original line number	Diff line number	Diff line
		@@ -2,7 +2,7 @@ CXX=g++
		CXXFLAGS=-std=c++11 -Wall -O3 -fopenmp -msse4.1
		LDFLAGS=-lm -lz

		cpp_source=sequence_batch.cc index.cc candidate_processor.cc ksw.cc output_tools.cc chromap.cc
		cpp_source=sequence_batch.cc index.cc candidate_processor.cc feature_barcode_matrix.cc ksw.cc output_tools.cc chromap.cc
		src_dir=src
		objs_dir=objs
		objs+=$(patsubst %.cc,$(objs_dir)/%.o,$(cpp_source))

src/chromap.cc

+20 −321

Original line number	Diff line number	Diff line
		@@ -11,6 +11,7 @@
		#include <limits>
		#include <random>
		#include <sstream>
		#include <type_traits>

		#include "candidate_processor.h"
		#include "cxxopts.hpp"
		@@ -367,321 +368,6 @@ bool Chromap<MappingRecord>::CorrectBarcodeAt(
		}
		}

		template <typename MappingRecord>
		uint32_t Chromap<MappingRecord>::CallPeaks(uint16_t coverage_threshold,
		uint32_t num_reference_sequences,
		const SequenceBatch &reference) {
		return 0;
		}

		template <>
		uint32_t Chromap<PairedEndMappingWithBarcode>::CallPeaks(
		uint16_t coverage_threshold, uint32_t num_reference_sequences,
		const SequenceBatch &reference) {
		double real_start_time = GetRealTime();
		std::vector<std::vector<PairedEndMappingWithBarcode>> &mappings =
		allocate_multi_mappings_
		? allocated_mappings_on_diff_ref_seqs_
		: (remove_pcr_duplicates_ ? deduped_mappings_on_diff_ref_seqs_
		: mappings_on_diff_ref_seqs_);
		// Build pileup
		for (uint32_t ri = 0; ri < num_reference_sequences; ++ri) {
		pileup_on_diff_ref_seqs_.emplace_back(std::vector<uint16_t>());
		pileup_on_diff_ref_seqs_[ri].assign(reference.GetSequenceLengthAt(ri), 0);
		for (size_t mi = 0; mi < mappings[ri].size(); ++mi) {
		for (uint16_t pi = 0; pi < mappings[ri][mi].fragment_length_; ++pi) {
		++pileup_on_diff_ref_seqs_[ri]
		[mappings[ri][mi].GetStartPosition() + pi];
		}
		}
		}
		std::cerr << "Built pileup in " << GetRealTime() - real_start_time << "s.\n";
		real_start_time = GetRealTime();
		// Call and save peaks
		tree_extras_on_diff_ref_seqs_.clear();
		tree_info_on_diff_ref_seqs_.clear();
		tree_extras_on_diff_ref_seqs_.reserve(num_reference_sequences);
		tree_info_on_diff_ref_seqs_.reserve(num_reference_sequences);
		uint32_t peak_count = 0;
		for (uint32_t ri = 0; ri < num_reference_sequences; ++ri) {
		tree_extras_on_diff_ref_seqs_.emplace_back(std::vector<uint32_t>());
		tree_extras_on_diff_ref_seqs_[ri].reserve(
		reference.GetSequenceLengthAt(ri) / 100);
		peaks_on_diff_ref_seqs_.emplace_back(std::vector<Peak>());
		uint32_t peak_start_position = 0;
		uint16_t peak_length = 0;
		for (size_t pi = 0; pi < reference.GetSequenceLengthAt(ri); ++pi) {
		if (pileup_on_diff_ref_seqs_[ri][pi] >= coverage_threshold) {
		if (peak_length == 0) { // start a new peak
		peak_start_position = pi;
		}
		++peak_length; // extend the peak
		} else if (peak_length > 0) { // save the previous peak
		// TODO(Haowen): improve peak calling
		peaks_on_diff_ref_seqs_[ri].emplace_back(
		Peak{peak_start_position, peak_length, peak_count});
		tree_extras_on_diff_ref_seqs_[ri].emplace_back(0);
		feature_barcode_matrix_writer_.OutputPeaks(peak_start_position,
		peak_length, ri, reference);
		++peak_count;
		peak_length = 0;
		}
		}
		BuildAugmentedTreeForPeaks(ri);
		}
		std::cerr << "Call peaks and built peak augmented tree in "
		<< GetRealTime() - real_start_time << "s.\n";
		// Output feature matrix
		return peak_count;
		}

		template <typename MappingRecord>
		void Chromap<MappingRecord>::OutputFeatureMatrix(
		uint32_t num_sequences, const SequenceBatch &reference) {}

		template <>
		void Chromap<PairedEndMappingWithBarcode>::OutputFeatureMatrix(
		uint32_t num_sequences, const SequenceBatch &reference) {
		uint32_t num_peaks = 0;
		if (cell_by_bin_) {
		feature_barcode_matrix_writer_.OutputPeaks(bin_size_, num_sequences,
		reference);
		for (uint32_t i = 0; i < num_sequences; ++i) {
		uint32_t ref_seq_length = reference.GetSequenceLengthAt(i);
		num_peaks += ref_seq_length / bin_size_;
		if (ref_seq_length % bin_size_ != 0) {
		++num_peaks;
		}
		}
		} else {
		num_peaks = CallPeaks(depth_cutoff_to_call_peak_, num_sequences, reference);
		}
		std::vector<std::vector<PairedEndMappingWithBarcode>> &mappings =
		allocate_multi_mappings_
		? allocated_mappings_on_diff_ref_seqs_
		: (remove_pcr_duplicates_ ? deduped_mappings_on_diff_ref_seqs_
		: mappings_on_diff_ref_seqs_);
		double real_start_time = GetRealTime();
		// First pass to index barcodes
		uint32_t barcode_index = 0;
		for (uint32_t rid = 0; rid < num_sequences; ++rid) {
		for (uint32_t mi = 0; mi < mappings[rid].size(); ++mi) {
		uint64_t barcode_key = mappings[rid][mi].cell_barcode_;
		khiter_t barcode_index_table_iterator =
		kh_get(k64_seq, barcode_index_table_, barcode_key);
		if (barcode_index_table_iterator == kh_end(barcode_index_table_)) {
		int khash_return_code;
		barcode_index_table_iterator = kh_put(k64_seq, barcode_index_table_,
		barcode_key, &khash_return_code);
		assert(khash_return_code != -1 && khash_return_code != 0);
		kh_value(barcode_index_table_, barcode_index_table_iterator) =
		barcode_index;
		++barcode_index;
		feature_barcode_matrix_writer_.AppendBarcodeOutput(barcode_key);
		}
		}
		}
		std::cerr << "Index and output barcodes in "
		<< GetRealTime() - real_start_time << "s.\n";
		real_start_time = GetRealTime();
		// Second pass to generate matrix
		khash_t(kmatrix) *matrix = kh_init(kmatrix);
		std::vector<uint32_t> overlapped_peak_indices;
		for (uint32_t rid = 0; rid < num_sequences; ++rid) {
		for (uint32_t mi = 0; mi < mappings[rid].size(); ++mi) {
		uint64_t barcode_key = mappings[rid][mi].cell_barcode_;
		khiter_t barcode_index_table_iterator =
		kh_get(k64_seq, barcode_index_table_, barcode_key);
		uint64_t barcode_index =
		kh_value(barcode_index_table_, barcode_index_table_iterator);
		overlapped_peak_indices.clear();
		if (cell_by_bin_) {
		GetNumOverlappedBins(rid, mappings[rid][mi].GetStartPosition(),
		mappings[rid][mi].GetEndPosition() -
		mappings[rid][mi].GetStartPosition(),
		num_sequences, reference, overlapped_peak_indices);
		} else {
		GetNumOverlappedPeaks(rid, mappings[rid][mi], overlapped_peak_indices);
		}
		size_t num_overlapped_peaks = overlapped_peak_indices.size();
		for (size_t pi = 0; pi < num_overlapped_peaks; ++pi) {
		uint32_t peak_index = overlapped_peak_indices[pi];
		uint64_t entry_index = (barcode_index << 32) \| peak_index;
		khiter_t matrix_iterator = kh_get(kmatrix, matrix, entry_index);
		if (matrix_iterator == kh_end(matrix)) {
		int khash_return_code;
		matrix_iterator =
		kh_put(kmatrix, matrix, entry_index, &khash_return_code);
		assert(khash_return_code != -1 && khash_return_code != 0);
		kh_value(matrix, matrix_iterator) = 1;
		} else {
		kh_value(matrix, matrix_iterator) += 1;
		}
		}
		}
		}
		std::cerr << "Generate feature matrix in " << GetRealTime() - real_start_time
		<< "s.\n";
		// Output matrix
		real_start_time = GetRealTime();
		feature_barcode_matrix_writer_.WriteMatrixOutputHead(
		num_peaks, kh_size(barcode_index_table_), kh_size(matrix));
		uint64_t key;
		uint32_t value;
		std::vector<std::pair<uint64_t, uint32_t>> feature_matrix;
		feature_matrix.reserve(kh_size(matrix));
		// kh_foreach(matrix, key, value,
		// output_tools_->AppendMatrixOutput((uint32_t)key, (uint32_t)(key >> 32),
		// value));
		kh_foreach(matrix, key, value, feature_matrix.emplace_back(key, value));
		kh_destroy(kmatrix, matrix);
		std::sort(feature_matrix.begin(), feature_matrix.end());
		for (size_t i = 0; i < feature_matrix.size(); ++i) {
		feature_barcode_matrix_writer_.AppendMatrixOutput(
		(uint32_t)feature_matrix[i].first,
		(uint32_t)(feature_matrix[i].first >> 32), feature_matrix[i].second);
		}
		std::cerr << "Output feature matrix in " << GetRealTime() - real_start_time
		<< "s.\n";
		}

		template <typename MappingRecord>
		void Chromap<MappingRecord>::GetNumOverlappedBins(
		uint32_t rid, uint32_t start_position, uint16_t mapping_length,
		uint32_t num_sequences, const SequenceBatch &reference,
		std::vector<uint32_t> &overlapped_peak_indices) {
		uint32_t bin_index = 0;
		for (uint32_t i = 0; i < rid; ++i) {
		uint32_t ref_seq_length = reference.GetSequenceLengthAt(i);
		bin_index += ref_seq_length / bin_size_;
		if (ref_seq_length % bin_size_ != 0) {
		++bin_index;
		}
		}
		bin_index += start_position / bin_size_;
		overlapped_peak_indices.emplace_back(bin_index);
		uint32_t max_num_overlapped_bins = mapping_length / bin_size_ + 2;
		for (uint32_t i = 0; i < max_num_overlapped_bins; ++i) {
		if (start_position + mapping_length - 1 >=
		(bin_index + 1 + i) * bin_size_) {
		overlapped_peak_indices.emplace_back(bin_index + 1 + i);
		}
		}
		}

		template <typename MappingRecord>
		void Chromap<MappingRecord>::BuildAugmentedTreeForPeaks(uint32_t ref_id) {
		// std::sort(mappings.begin(), mappings.end(), IntervalLess());
		int max_level = 0;
		size_t i, last_i = 0; // last_i points to the rightmost node in the tree
		uint32_t last = 0; // last is the max value at node last_i
		int k;
		std::vector<Peak> &peaks = peaks_on_diff_ref_seqs_[ref_id];
		std::vector<uint32_t> &extras = tree_extras_on_diff_ref_seqs_[ref_id];
		if (peaks.size() == 0) {
		max_level = -1;
		}
		for (i = 0; i < peaks.size(); i += 2) {
		last_i = i;
		// last = mappings[i].max = mappings[i].en; // leaves (i.e. at level 0)
		last = extras[i] =
		peaks[i].start_position + peaks[i].length; // leaves (i.e. at level 0)
		}
		for (k = 1; 1LL << k <= (int64_t)peaks.size();
		++k) { // process internal nodes in the bottom-up order
		size_t x = 1LL << (k - 1);
		size_t i0 = (x << 1) - 1;
		size_t step = x << 2; // i0 is the first node
		for (i = i0; i < peaks.size();
		i += step) { // traverse all nodes at level k
		uint32_t el = extras[i - x]; // max value of the left child
		uint32_t er =
		i + x < peaks.size() ? extras[i + x] : last; // of the right child
		uint32_t e = peaks[i].start_position + peaks[i].length;
		e = e > el ? e : el;
		e = e > er ? e : er;
		extras[i] = e; // set the max value for node i
		}
		last_i =
		last_i >> k & 1
		? last_i - x
		: last_i +
		x; // last_i now points to the parent of the original last_i
		if (last_i < peaks.size() &&
		extras[last_i] > last) // update last accordingly
		last = extras[last_i];
		}
		max_level = k - 1;
		tree_info_on_diff_ref_seqs_.emplace_back(max_level, peaks.size());
		}

		template <typename MappingRecord>
		uint32_t Chromap<MappingRecord>::GetNumOverlappedPeaks(
		uint32_t ref_id, const MappingRecord &mapping,
		std::vector<uint32_t> &overlapped_peak_indices) {
		int t = 0;
		StackCell stack[64];
		// out.clear();
		overlapped_peak_indices.clear();
		int num_overlapped_peaks = 0;
		int max_level = tree_info_on_diff_ref_seqs_[ref_id].first;
		uint32_t num_tree_nodes = tree_info_on_diff_ref_seqs_[ref_id].second;
		std::vector<Peak> &peaks = peaks_on_diff_ref_seqs_[ref_id];
		std::vector<uint32_t> &extras = tree_extras_on_diff_ref_seqs_[ref_id];
		// uint32_t interval_start = mapping.fragment_start_position;
		uint32_t interval_start =
		mapping.GetStartPosition() > (uint32_t)multi_mapping_allocation_distance_
		? mapping.GetStartPosition() - multi_mapping_allocation_distance_
		: 0;
		uint32_t interval_end =
		mapping.GetEndPosition() + (uint32_t)multi_mapping_allocation_distance_;
		stack[t++] = StackCell(max_level, (1LL << max_level) - 1,
		0); // push the root; this is a top down traversal
		while (
		t) { // the following guarantees that numbers in out[] are always sorted
		StackCell z = stack[--t];
		if (z.k <=
		3) { // we are in a small subtree; traverse every node in this subtree
		size_t i, i0 = z.x >> z.k << z.k, i1 = i0 + (1LL << (z.k + 1)) - 1;
		if (i1 >= num_tree_nodes) {
		i1 = num_tree_nodes;
		}
		for (i = i0; i < i1 && peaks[i].start_position < interval_end; ++i) {
		if (interval_start <
		peaks[i].start_position +
		peaks[i].length) { // if overlap, append to out[]
		// out.push_back(i);
		overlapped_peak_indices.emplace_back(peaks[i].index);
		++num_overlapped_peaks;
		}
		}
		} else if (z.w == 0) { // if left child not processed
		size_t y =
		z.x - (1LL << (z.k - 1)); // the left child of z.x; NB: y may be out
		// of range (i.e. y>=a.size())
		stack[t++] = StackCell(
		z.k, z.x,
		1); // re-add node z.x, but mark the left child having been processed
		if (y >= num_tree_nodes \|\|
		extras[y] > interval_start) // push the left child if y is out of
		// range or may overlap with the query
		stack[t++] = StackCell(z.k - 1, y, 0);
		} else if (z.x < num_tree_nodes &&
		peaks[z.x].start_position <
		interval_end) { // need to push the right child
		if (interval_start < peaks[z.x].start_position + peaks[z.x].length) {
		// out.push_back(z.x); // test if z.x overlaps the query; if yes, append
		// to out[]
		overlapped_peak_indices.emplace_back(peaks[z.x].index);
		++num_overlapped_peaks;
		}
		stack[t++] = StackCell(z.k - 1, z.x + (1LL << (z.k - 1)),
		0); // push the right child
		}
		}
		return num_overlapped_peaks;
		}

		template <typename MappingRecord>
		uint32_t Chromap<MappingRecord>::LoadPairedEndReadsWithBarcodes(
		SequenceBatch &read_batch1, SequenceBatch &read_batch2,
		@@ -1653,12 +1339,25 @@ void Chromap<MappingRecord>::MapPairedEndReads() {
		OutputMappings(num_reference_sequences, reference, mappings);
		}

		if (!is_bulk_data_ && !matrix_output_prefix_.empty()) {
		feature_barcode_matrix_writer_.InitializeMatrixOutput(
		matrix_output_prefix_);
		OutputFeatureMatrix(num_reference_sequences, reference);
		feature_barcode_matrix_writer_.FinalizeMatrixOutput();
		}
		// Temporarily disable feature matrix output. Do not delete the following
		// commented code.
		//if (!is_bulk_data_ && !matrix_output_prefix_.empty()) {
		// if constexpr (std::is_same<MappingRecord,
		// PairedEndMappingWithBarcode>::value) {
		// FeatureBarcodeMatrix feature_barcode_matrix(
		// cell_by_bin_, bin_size_, multi_mapping_allocation_distance_,
		// depth_cutoff_to_call_peak_);
		// std::vector<std::vector<PairedEndMappingWithBarcode>> &mappings =
		// allocate_multi_mappings_
		// ? allocated_mappings_on_diff_ref_seqs_
		// : (remove_pcr_duplicates_ ? deduped_mappings_on_diff_ref_seqs_
		// : mappings_on_diff_ref_seqs_);

		// feature_barcode_matrix.OutputFeatureMatrix(num_reference_sequences,
		// reference, mappings,
		// matrix_output_prefix_);
		// }
		//}
		}

		output_tools_.FinalizeMappingOutput();

src/chromap.h

+4 −68

Original line number	Diff line number	Diff line
		#ifndef CHROMAP_H_
		#define CHROMAP_H_

		#include <sys/resource.h>
		#include <sys/time.h>

		#include <memory>
		#include <random>
		#include <string>
		#include <tuple>
		#include <vector>

		#include "feature_barcode_matrix_writer.h"
		#include "feature_barcode_matrix.h"
		#include "index.h"
		#include "index_parameters.h"
		#include "khash.h"
		@@ -25,17 +22,6 @@
		#define CHROMAP_VERSION "0.1.4-r285"

		namespace chromap {
		struct uint128_t {
		uint64_t first;
		uint64_t second;
		};

		struct StackCell {
		size_t x; // node
		int k, w; // k: level; w: 0 if left child hasn't been processed
		StackCell(){};
		StackCell(int k_, size_t x_, int w_) : x(x_), k(k_), w(w_){};
		};

		struct _mm_history {
		unsigned int timestamp;
		@@ -45,36 +31,6 @@ struct _mm_history {
		uint32_t repetitive_seed_length;
		};

		struct Peak {
		uint32_t start_position;
		uint16_t length;
		uint32_t index;
		};

		KHASH_MAP_INIT_INT64(k128, uint128_t);
		KHASH_MAP_INIT_INT64(k64_seq, uint64_t);
		KHASH_SET_INIT_INT(k32_set);
		KHASH_MAP_INIT_INT64(kmatrix, uint32_t);

		struct BarcodeWithQual {
		uint32_t corrected_base_index1;
		char correct_base1;
		uint32_t corrected_base_index2;
		char correct_base2;
		double score;
		bool operator>(const BarcodeWithQual &b) const {
		return std::tie(score, corrected_base_index1, correct_base1,
		corrected_base_index2, correct_base2) >
		std::tie(b.score, b.corrected_base_index1, b.correct_base1,
		b.corrected_base_index2, b.correct_base2);
		}
		};

		#define SortMappingWithoutBarcode(m) \
		(((((m).fragment_start_position_ << 16) \| (m).fragment_length_) << 8) \| \
		(m).mapq_)
		//#define SortMappingWithoutBarcode(m) (m)

		class ChromapDriver {
		public:
		ChromapDriver() {}
		@@ -95,7 +51,6 @@ class Chromap {
		barcode_lookup_table_ = NULL;
		barcode_whitelist_lookup_table_ = NULL;
		barcode_histogram_ = NULL;
		barcode_index_table_ = NULL;
		}

		// For mapping
		@@ -158,7 +113,6 @@ class Chromap {
		barcode_lookup_table_ = kh_init(k64_seq);
		barcode_whitelist_lookup_table_ = kh_init(k64_seq);
		barcode_histogram_ = kh_init(k64_seq);
		barcode_index_table_ = kh_init(k64_seq);

		NUM_VPU_LANES_ = 0;
		if (error_threshold_ < 8) {
		@@ -181,9 +135,7 @@ class Chromap {
		if (barcode_histogram_ != NULL) {
		kh_destroy(k64_seq, barcode_histogram_);
		}
		if (barcode_index_table_ != NULL) {
		kh_destroy(k64_seq, barcode_index_table_);
		}

		if (barcode_lookup_table_ != NULL) {
		kh_destroy(k64_seq, barcode_lookup_table_);
		}
		@@ -430,18 +382,7 @@ class Chromap {
		bool CorrectBarcodeAt(uint32_t barcode_index, SequenceBatch &barcode_batch,
		uint64_t &num_barcode_in_whitelist,
		uint64_t &num_corrected_barcode);
		uint32_t CallPeaks(uint16_t coverage_threshold,
		uint32_t num_reference_sequences,
		const SequenceBatch &reference);
		void OutputFeatureMatrix(uint32_t num_sequences,
		const SequenceBatch &reference);
		void GetNumOverlappedBins(uint32_t rid, uint32_t start_position,
		uint16_t mapping_length, uint32_t num_sequences,
		const SequenceBatch &reference,
		std::vector<uint32_t> &overlapped_peak_indices);
		uint32_t GetNumOverlappedPeaks(
		uint32_t ref_id, const MappingRecord &mapping,
		std::vector<uint32_t> &overlapped_peak_indices);

		void BuildAugmentedTreeForPeaks(uint32_t ref_id);
		void OutputMappingsInVector(
		uint8_t mapq_threshold, uint32_t num_reference_sequences,
		@@ -546,7 +487,6 @@ class Chromap {
		// (max_level, # nodes)
		std::vector<std::pair<int, uint32_t> > tree_info_on_diff_ref_seqs_;
		OutputTools<MappingRecord> output_tools_;
		FeatureBarcodeMatrixWriter feature_barcode_matrix_writer_;
		// For mapping stats.
		uint64_t num_candidates_ = 0;
		uint64_t num_mappings_ = 0;
		@@ -561,13 +501,9 @@ class Chromap {
		uint64_t num_corrected_barcode_ = 0;
		uint32_t barcode_length_ = 0;
		khash_t(k64_seq) * barcode_histogram_;
		khash_t(k64_seq) * barcode_index_table_;
		bool skip_barcode_check_ = false;

		// For peak calling
		std::vector<std::vector<uint16_t> > pileup_on_diff_ref_seqs_;
		std::vector<std::vector<Peak> > peaks_on_diff_ref_seqs_;
		};

		} // namespace chromap

		#endif // CHROMAP_H_

src/feature_barcode_matrix.cc

0 → 100644

+325 −0

File added.

Preview size limit exceeded, changes collapsed.

src/feature_barcode_matrix.h

0 → 100644

+87 −0

Original line number	Diff line number	Diff line
		#ifndef FEATURE_BARCODE_MATRIX_H_
		#define FEATURE_BARCODE_MATRIX_H_

		#include <assert.h>

		#include <algorithm>
		#include <cinttypes>
		#include <cstring>
		#include <functional>
		#include <iostream>
		#include <string>
		#include <vector>

		#include "bed_mapping.h"
		#include "feature_barcode_matrix_writer.h"
		#include "khash.h"
		#include "utils.h"

		namespace chromap {

		struct Peak {
		uint32_t start_position;
		uint16_t length;
		uint32_t index;
		};

		class FeatureBarcodeMatrix {
		public:
		FeatureBarcodeMatrix(bool cell_by_bin, int bin_size, int overlap_distance,
		uint16_t depth_cutoff_to_call_peak)
		: cell_by_bin_(cell_by_bin),
		bin_size_(bin_size),
		overlap_distance_(overlap_distance),
		depth_cutoff_to_call_peak_(depth_cutoff_to_call_peak) {
		barcode_index_table_ = kh_init(k64_seq);
		}

		~FeatureBarcodeMatrix() {
		if (barcode_index_table_ != NULL) {
		kh_destroy(k64_seq, barcode_index_table_);
		}
		}

		void OutputFeatureMatrix(
		uint32_t num_sequences, const SequenceBatch &reference,
		const std::vector<std::vector<PairedEndMappingWithBarcode>> &mappings,
		const std::string &matrix_output_prefix);

		private:
		void BuildAugmentedTreeForPeaks(uint32_t ref_id);

		uint32_t GetNumOverlappedPeaks(
		uint32_t ref_id, const PairedEndMappingWithBarcode &mapping,
		std::vector<uint32_t> &overlapped_peak_indices);

		void GetNumOverlappedBins(uint32_t rid, uint32_t start_position,
		uint16_t mapping_length, uint32_t num_sequences,
		const SequenceBatch &reference,
		std::vector<uint32_t> &overlapped_peak_indices);

		uint32_t CallPeaks(
		uint16_t coverage_threshold, uint32_t num_reference_sequences,
		const SequenceBatch &reference,
		const std::vector<std::vector<PairedEndMappingWithBarcode>> &mappings);

		const bool cell_by_bin_;
		const int bin_size_;
		const int overlap_distance_;
		const uint16_t depth_cutoff_to_call_peak_;

		khash_t(k64_seq) * barcode_index_table_;
		// (max_level, # nodes)
		std::vector<std::pair<int, uint32_t>> tree_info_on_diff_ref_seqs_;

		// max
		std::vector<std::vector<uint32_t>> tree_extras_on_diff_ref_seqs_;

		// For peak calling.
		std::vector<std::vector<uint16_t>> pileup_on_diff_ref_seqs_;
		std::vector<std::vector<Peak>> peaks_on_diff_ref_seqs_;

		FeatureBarcodeMatrixWriter feature_barcode_matrix_writer_;
		};

		} // namespace chromap

		#endif // FEATURE_BARCODE_MATRIX_H_

Admin message