xsample03.cpp

Seach for human mutation (SNP) in within chr1.

Seach for human mutation (SNP) in within chr1. Benchmark comaprison of different methods

See also
  bm::sparse_vector
  bm::rsc_sparse_vector
  bm::sparse_vector_scanner

/*
Copyright(c) 2018 Anatoliy Kuznetsov(anatoliy_kuznetsov at yahoo.com)
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
For more information please visit: http://bitmagic.io
*/
/** \example xsample03.cpp
Seach for human mutation (SNP) in within chr1.
Benchmark comaprison of different methods
\sa bm::sparse_vector
\sa bm::rsc_sparse_vector
\sa bm::sparse_vector_scanner
*/
/*! \file xsample03.cpp
\brief Example: SNP search in human genome
Brief description of used method:
1. Parse SNP chromosome report and extract information about SNP number and
location in the chromosome
2. Use this information to build bit-transposed sparse_vector<>
where vector position matches chromosome position and SNP ids (aka rsid)
is kept as a bit-transposed matrix
3. Build rank-select compressed sparse vector, dropping all NULL columns
(this data format is pretty sparse, since number of SNPs is significantly
less than number of chromosome bases (1:5 or less)
Use memory report to understand memory footprint for each form of storage
4. Run benchmarks searching for 500 randomly selected SNPs using
- bm::sparse_vector<>
- bm::rsc_sparse_vector<>
- std::vector<pair<unsigned, unsigned> >
This example should be useful for construction of compressed columnar
tables with parallel search capabilities.
*/
#include <iostream>
#include <sstream>
#include <chrono>
#include <regex>
#include <time.h>
#include <stdio.h>
#include <vector>
#include <chrono>
#include <map>
#include <utility>
using namespace std;
#include "bm.h"
#include "bmalgo.h"
#include "bmserial.h"
#include "bmrandom.h"
#include "bmsparsevec.h"
#include "bmsparsevec_compr.h"
#include "bmsparsevec_algo.h"
#include "bmsparsevec_serial.h"
#include "bmalgo_similarity.h"
#include "bmsparsevec_util.h"
#include "bmdbg.h"
#include "bmtimer.h"
#include "bmundef.h" /* clear the pre-proc defines from BM */
static
void show_help()
{
std::cerr
<< "BitMagic SNP Search Sample Utility (c) 2018" << std::endl
<< "-isnp file-name -- input set file (SNP FASTA) to parse" << std::endl
<< "-svout spase vector output -- sparse vector name to save" << std::endl
<< "-rscout rs-compressed spase vector output -- name to save" << std::endl
<< "-svin sparse vector input -- sparse vector file name to load " << std::endl
<< "-rscin rs-compressed sparse vector input -- file name to load " << std::endl
<< "-diag -- run diagnostics" << std::endl
<< "-timing -- collect timings" << std::endl
;
}
// Arguments
//
std::string sv_out_name;
std::string rsc_out_name;
std::string sv_in_name;
std::string rsc_in_name;
std::string isnp_name;
bool is_diag = false;
bool is_timing = false;
bool is_bench = false;
static
int parse_args(int argc, char *argv[])
{
for (int i = 1; i < argc; ++i)
{
std::string arg = argv[i];
if ((arg == "-h") || (arg == "--help"))
{
show_help();
return 0;
}
if (arg == "-svout" || arg == "--svout")
{
if (i + 1 < argc)
{
sv_out_name = argv[++i];
}
else
{
std::cerr << "Error: -svout requires file name" << std::endl;
return 1;
}
continue;
}
if (arg == "-rscout" || arg == "--rscout")
{
if (i + 1 < argc)
{
rsc_out_name = argv[++i];
}
else
{
std::cerr << "Error: -rscout requires file name" << std::endl;
return 1;
}
continue;
}
if (arg == "-svin" || arg == "--svin")
{
if (i + 1 < argc)
{
sv_in_name = argv[++i];
}
else
{
std::cerr << "Error: -svin requires file name" << std::endl;
return 1;
}
continue;
}
if (arg == "-rscin" || arg == "--rscin")
{
if (i + 1 < argc)
{
rsc_in_name = argv[++i];
}
else
{
std::cerr << "Error: -rscin requires file name" << std::endl;
return 1;
}
continue;
}
if (arg == "-isnp" || arg == "--isnp" || arg == "-snp" || arg == "--snp")
{
if (i + 1 < argc)
{
isnp_name = argv[++i];
}
else
{
std::cerr << "Error: -isnp requires file name" << std::endl;
return 1;
}
continue;
}
if (arg == "-diag" || arg == "--diag" || arg == "-d" || arg == "--d")
is_diag = true;
if (arg == "-timing" || arg == "--timing" || arg == "-t" || arg == "--t")
is_timing = true;
if (arg == "-bench" || arg == "--bench" || arg == "-b" || arg == "--b")
is_bench = true;
} // for i
return 0;
}
// Global types
//
typedef bm::sparse_vector<unsigned, bm::bvector<> > sparse_vector_u32;
typedef bm::rsc_sparse_vector<unsigned, sparse_vector_u32 > rsc_sparse_vector_u32;
typedef std::vector<std::pair<unsigned, unsigned> > vector_pairs;
// ----------------------------------------------------------------------------
bm::chrono_taker<>::duration_map_type timing_map;
// SNP report format parser (extraction and transformation)
// Parser extracts SNP id (rsid) and positio on genome to build
// sparse vector where index (position in vector) metches position on the
// chromosome 1.
//
static
int load_snp_report(const std::string& fname, sparse_vector_u32& sv)
{
bm::chrono_taker tt1(cout, "1. parse input SNP chr report", 1, &timing_map);
std::ifstream fin(fname.c_str(), std::ios::in);
if (!fin.good())
return -1;
unsigned rs_id, rs_pos;
size_t idx;
std::string line;
std::string delim = " \t";
std::regex reg("\\s+");
std::sregex_token_iterator it_end;
bm::bvector<> bv_rs;
bv_rs.init();
unsigned rs_cnt = 0;
for (unsigned i = 0; std::getline(fin, line); ++i)
{
if (line.empty() ||
!isdigit(line.front())
)
continue;
// regex based tokenizer
std::sregex_token_iterator it(line.begin(), line.end(), reg, -1);
std::vector<std::string> line_vec(it, it_end);
if (line_vec.empty())
continue;
// parse columns of interest
try
{
rs_id = unsigned(std::stoul(line_vec.at(0), &idx));
if (bv_rs.test(rs_id))
{
continue;
}
rs_pos = unsigned(std::stoul(line_vec.at(11), &idx));
bv_rs.set_bit_no_check(rs_id);
sv.set(rs_pos, rs_id);
++rs_cnt;
}
catch (std::exception& /*ex*/)
{
continue; // detailed disgnostics commented out
// error detected, because some columns are sometimes missing
// just ignore it
//
/*
std::cerr << ex.what() << "; ";
std::cerr << "rs=" << line_vec.at(0) << " pos=" << line_vec.at(11) << std::endl;
continue;
*/
}
if (rs_cnt % (4 * 1024) == 0)
std::cout << "\r" << rs_cnt << " / " << i; // PROGRESS report
} // for i
std::cout << std::endl;
std::cout << "SNP count=" << rs_cnt << std::endl;
return 0;
}
// Generate random subset of random values from a sparse vector
//
static
void generate_random_subset(const sparse_vector_u32& sv, std::vector<unsigned>& vect, unsigned count)
{
const sparse_vector_u32::bvector_type* bv_null = sv.get_null_bvector();
bm::random_subset<bm::bvector<> > rand_sampler;
bm::bvector<> bv_sample;
rand_sampler.sample(bv_sample, *bv_null, count);
bm::bvector<>::enumerator en = bv_sample.first();
for (; en.valid(); ++en)
{
unsigned idx = *en;
unsigned v = sv[idx];
vect.push_back(v);
}
}
// build std::vector of pairs (position to rs)
//
static
void build_vector_pairs(const sparse_vector_u32& sv, vector_pairs& vp)
{
sparse_vector_u32::const_iterator it = sv.begin();
sparse_vector_u32::const_iterator it_end = sv.end();
for (; it != it_end; ++it)
{
if (!it.is_null())
{
std::pair<unsigned, unsigned> pos2rs = std::make_pair(it.pos(), it.value());
vp.push_back(pos2rs);
}
}
}
// O(N) -- O(N/2) linear search in vector of pairs (position - rsid)
//
static
bool search_vector_pairs(const vector_pairs& vp, unsigned rs_id, unsigned& pos)
{
for (unsigned i = 0; i < vp.size(); ++i)
{
if (vp[i].second == rs_id)
{
pos = vp[i].first;
return true;
}
}
return false;
}
// SNP search benchmark
// Search for SNPs using different data structures (Bitmagic and STL)
//
// An extra step is verification, to make sure all methods are consistent
//
static
void run_benchmark(const sparse_vector_u32& sv, const rsc_sparse_vector_u32& csv)
{
const unsigned rs_sample_count = 2000;
std::vector<unsigned> rs_vect;
generate_random_subset(sv, rs_vect, rs_sample_count);
if (rs_vect.empty())
{
std::cerr << "Benchmark subset empty!" << std::endl;
return;
}
// build traditional sparse vector
vector_pairs vp;
build_vector_pairs(sv, vp);
// search result bit-vectors
//
bm::bvector<> bv_found1;
bm::bvector<> bv_found2;
bm::bvector<> bv_found3;
bv_found1.init(); bv_found2.init(); bv_found3.init();// pre-initialize vectors
if (!sv.empty())
{
bm::chrono_taker tt1(cout, "09. rs search (sv)", unsigned(rs_vect.size()), &timing_map);
// scanner class
// it's important to keep scanner class outside the loop to avoid
// unnecessary re-allocs and construction costs
//
bm::sparse_vector_scanner<sparse_vector_u32> scanner;
for (unsigned i = 0; i < rs_vect.size(); ++i)
{
unsigned rs_id = rs_vect[i];
unsigned rs_pos;
bool found = scanner.find_eq(sv, rs_id, rs_pos);
if (found)
{
bv_found1.set_bit_no_check(rs_pos);
}
else
{
std::cout << "Error: rs_id = " << rs_id << " not found!" << std::endl;
}
} // for
}
if (!csv.empty())
{
bm::chrono_taker tt1(cout, "10. rs search (rsc_sv)", unsigned(rs_vect.size()), &timing_map);
bm::sparse_vector_scanner<rsc_sparse_vector_u32> scanner; // scanner class
for (unsigned i = 0; i < rs_vect.size(); ++i)
{
unsigned rs_id = rs_vect[i];
unsigned rs_pos;
bool found = scanner.find_eq(csv, rs_id, rs_pos);
if (found)
{
bv_found2.set_bit_no_check(rs_pos);
}
else
{
std::cout << "rs_id = " << rs_id << " not found!" << std::endl;
}
} // for
}
if (vp.size())
{
bm::chrono_taker tt1(cout, "11. rs search (std::vector<>)", unsigned(rs_vect.size()), &timing_map);
for (unsigned i = 0; i < rs_vect.size(); ++i)
{
unsigned rs_id = rs_vect[i];
unsigned rs_pos;
bool found = search_vector_pairs(vp, rs_id, rs_pos);
if (found)
{
bv_found3.set_bit_no_check(rs_pos);
}
else
{
std::cout << "rs_id = " << rs_id << " not found!" << std::endl;
}
} // for
}
// compare results from various methods (check integrity)
int res = bv_found1.compare(bv_found2);
if (res != 0)
{
std::cerr << "Error: search discrepancy (sparse search) detected!" << std::endl;
}
res = bv_found1.compare(bv_found3);
if (res != 0)
{
std::cerr << "Error: search discrepancy (std::vector<>) detected!" << std::endl;
}
}
int main(int argc, char *argv[])
{
if (argc < 3)
{
show_help();
return 1;
}
sparse_vector_u32 sv(bm::use_null);
rsc_sparse_vector_u32 csv;
try
{
auto ret = parse_args(argc, argv);
if (ret != 0)
return ret;
if (!isnp_name.empty())
{
auto res = load_snp_report(isnp_name, sv);
if (res != 0)
{
return res;
}
}
if (!sv_in_name.empty())
{
bm::chrono_taker tt1(cout, "02. Load sparse vector", 1, &timing_map);
file_load_svector(sv, sv_in_name);
}
// load rs-compressed sparse vector
if (!rsc_in_name.empty())
{
bm::chrono_taker tt1(cout, "03. Load rsc sparse vector", 1, &timing_map);
file_load_svector(csv, rsc_in_name);
}
// if rs-compressed vector is not available - build it on the fly
if (csv.empty() && !sv.empty())
{
sparse_vector_u32 sv2(bm::use_null);
{
bm::chrono_taker tt1(cout, "04. rs compress sparse vector", 1, &timing_map);
csv.load_from(sv);
}
{
bm::chrono_taker tt1(cout, "05. rs de-compress sparse vector", 1, &timing_map);
csv.load_to(sv2);
}
if (!sv.equal(sv2)) // diagnostics check (just in case)
{
std::cerr << "Error: rs-compressed vector check failed!" << std::endl;
return 1;
}
}
// save SV vector
if (!sv_out_name.empty() && !sv.empty())
{
bm::chrono_taker tt1(cout, "07. Save sparse vector", 1, &timing_map);
sv.optimize();
file_save_svector(sv, sv_out_name);
}
// save RS sparse vector
if (!rsc_out_name.empty() && !csv.empty())
{
bm::chrono_taker tt1(cout, "08. Save RS sparse vector", 1, &timing_map);
csv.optimize();
file_save_svector(csv, rsc_out_name);
}
if (is_diag) // print memory diagnostics
{
if (!sv.empty())
{
std::cout << std::endl
<< "sparse vector statistics:"
<< std::endl;
bm::print_svector_stat(std::cout, sv, true);
}
if (!csv.empty())
{
std::cout << std::endl
<< "RS compressed sparse vector statistics:"
<< std::endl;
bm::print_svector_stat(std::cout, csv, true);
}
}
if (is_bench) // run set of benchmarks
{
run_benchmark(sv, csv);
}
if (is_timing) // print all collected timings
{
std::cout << std::endl << "Performance:" << std::endl;
bm::chrono_taker<>::print_duration_map(cout, timing_map, bm::chrono_taker<>::ct_ops_per_sec);
}
}
catch (std::exception& ex)
{
std::cerr << "Error:" << ex.what() << std::endl;
return 1;
}
return 0;
}
bm.h
Compressed bit-vector bvector<> container, set algebraic methods, traversal iterators.
bmalgo.h
Algorithms for bvector<> (main include).
bmalgo_similarity.h
bmrandom.h
Generation of random subset.
bmserial.h
Serialization / compression of bvector<>. Set theoretical operations on compressed BLOBs.
bmsparsevec.h
Sparse constainer sparse_vector<> for integer types using bit-transposition transform.
bmsparsevec_algo.h
Algorithms for bm::sparse_vector.
bmsparsevec_compr.h
Compressed sparse container rsc_sparse_vector<> for integer types.
bmsparsevec_serial.h
Serialization for sparse_vector<>.
bmsparsevec_util.h
bmtimer.h
Timing utilities for benchmarking (internal).
bmundef.h
pre-processor un-defines to avoid global space pollution (internal)
bm::base_sparse_vector::get_null_bvector
const bvector_type * get_null_bvector() const BMNOEXCEPT
Get bit-vector of assigned values or NULL (if not constructed that way).
Definition bmbmatrix.h:451
bm::bvector::enumerator
Constant iterator designed to enumerate "ON" bits.
Definition bm.h:603
bm::bvector::iterator_base::valid
bool valid() const BMNOEXCEPT
Checks if iterator is still valid.
Definition bm.h:283
bm::bvector
Bitvector Bit-vector container with runtime compression of bits.
Definition bm.h:115
bm::bvector::test
bool test(size_type n) const BMNOEXCEPT
returns true if bit n is set and false is bit n is 0.
Definition bm.h:1502
bm::bvector::first
enumerator first() const
Returns enumerator pointing on the first non-zero bit.
Definition bm.h:1871
bm::bvector::init
void init()
Explicit post-construction initialization. Must be caled to make sure safe use of *_no_check() method...
Definition bm.h:2280
bm::bvector::compare
int compare(const bvector< Alloc > &bvect) const BMNOEXCEPT
Lexicographical comparison with a bitvector.
Definition bm.h:3744
bm::bvector::set_bit_no_check
void set_bit_no_check(size_type n)
Set bit without checking preconditions (size, etc).
Definition bm.h:4646
bm::chrono_taker
Utility class to collect performance measurements and statistics.
Definition bmtimer.h:41
bm::chrono_taker::duration_map_type
std::map< std::string, statistics > duration_map_type
test name to duration map
Definition bmtimer.h:66
bm::chrono_taker::ct_ops_per_sec
@ ct_ops_per_sec
Definition bmtimer.h:60
bm::chrono_taker::print_duration_map
static void print_duration_map(TOut &tout, const duration_map_type &dmap, format fmt=ct_time)
Definition bmtimer.h:150
bm::random_subset
Definition bmrandom.h:57
bm::random_subset::sample
void sample(BV &bv_out, const BV &bv_in, size_type sample_count)
Get random subset of input vector.
Definition bmrandom.h:140
bm::rsc_sparse_vector
Rank-Select compressed sparse vector.
Definition bmsparsevec_compr.h:59
bm::rsc_sparse_vector::load_from
void load_from(const sparse_vector_type &sv_src)
Load compressed vector from a sparse vector (with NULLs).
Definition bmsparsevec_compr.h:1369
bm::rsc_sparse_vector::load_to
void load_to(sparse_vector_type &sv) const
Exort compressed vector to a sparse vector (with NULLs).
Definition bmsparsevec_compr.h:1410
bm::rsc_sparse_vector::optimize
void optimize(bm::word_t *temp_block=0, typename bvector_type::optmode opt_mode=bvector_type::opt_compress, statistics *stat=0)
run memory optimization for all vector slices
Definition bmsparsevec_compr.h:1644
bm::rsc_sparse_vector::empty
bool empty() const BMNOEXCEPT
return true if vector is empty
Definition bmsparsevec_compr.h:379
bm::sparse_vector_scanner
algorithms for sparse_vector scan/search
Definition bmsparsevec_algo.h:582
bm::sparse_vector_scanner::find_eq
void find_eq(const SV &sv, value_type value, bvector_type &bv_out)
find all sparse vector elements EQ to search value
Definition bmsparsevec_algo.h:3071
bm::sparse_vector
succinct sparse vector with runtime compression using bit-slicing / transposition method
Definition bmsparsevec.h:87
bm::sparse_vector::equal
bool equal(const sparse_vector< Val, BV > &sv, bm::null_support null_able=bm::use_null) const BMNOEXCEPT
check if another sparse vector has the same content and size
Definition bmsparsevec.h:2300
bm::sparse_vector::push_back
void push_back(value_type v)
push value back into vector
Definition bmsparsevec.h:1883
bm::sparse_vector::set
void set(size_type idx, value_type v)
set specified element with bounds checking and automatic resize
Definition bmsparsevec.h:1849
bm::sparse_vector::end
const_iterator end() const BMNOEXCEPT
Provide const iterator access to the end.
Definition bmsparsevec.h:575
bm::sparse_vector::empty
bool empty() const BMNOEXCEPT
return true if vector is empty
Definition bmsparsevec.h:734
bm::sparse_vector< unsigned, bm::bvector<> >::bvector_type
bm::bvector<> bvector_type
Definition bmsparsevec.h:90
bm::sparse_vector::begin
const_iterator begin() const BMNOEXCEPT
Provide const iterator access to container content.
Definition bmsparsevec.h:2310
bm::sparse_vector::optimize
void optimize(bm::word_t *temp_block=0, typename bvector_type::optmode opt_mode=bvector_type::opt_compress, typename sparse_vector< Val, BV >::statistics *stat=0)
run memory optimization for all vector planes
Definition bmsparsevec.h:2148
bm::use_null
@ use_null
support "non-assigned" or "NULL" logic
Definition bmconst.h:230
timing_map
bm::chrono_taker ::duration_map_type timing_map
Definition sample11.cpp:46
main
int main(void)
Definition sample1.cpp:43
sparse_vector_u32
bm::sparse_vector< unsigned, bm::bvector<> > sparse_vector_u32
Definition xsample02.cpp:68
parse_args
static int parse_args(int argc, char *argv[])
Definition xsample03.cpp:110
vector_pairs
std::vector< std::pair< unsigned, unsigned > > vector_pairs
Definition xsample03.cpp:206
is_bench
bool is_bench
Definition xsample03.cpp:107
sv_in_name
std::string sv_in_name
Definition xsample03.cpp:102
rsc_in_name
std::string rsc_in_name
Definition xsample03.cpp:103
build_vector_pairs
static void build_vector_pairs(const sparse_vector_u32 &sv, vector_pairs &vp)
Definition xsample03.cpp:312
search_vector_pairs
static bool search_vector_pairs(const vector_pairs &vp, unsigned rs_id, unsigned &pos)
Definition xsample03.cpp:330
show_help
static void show_help()
Definition xsample03.cpp:81
isnp_name
std::string isnp_name
Definition xsample03.cpp:104
is_diag
bool is_diag
Definition xsample03.cpp:105
rsc_sparse_vector_u32
bm::rsc_sparse_vector< unsigned, sparse_vector_u32 > rsc_sparse_vector_u32
Definition xsample03.cpp:205
generate_random_subset
static void generate_random_subset(const sparse_vector_u32 &sv, std::vector< unsigned > &vect, unsigned count)
Definition xsample03.cpp:292
sv_out_name
std::string sv_out_name
Definition xsample03.cpp:100
rsc_out_name
std::string rsc_out_name
Definition xsample03.cpp:101
load_snp_report
static int load_snp_report(const std::string &fname, sparse_vector_u32 &sv)
Definition xsample03.cpp:218
run_benchmark
static void run_benchmark(const sparse_vector_u32 &sv, const rsc_sparse_vector_u32 &csv)
Definition xsample03.cpp:349
is_timing
bool is_timing
Definition xsample03.cpp:106