Skip to content
Snippets Groups Projects
Commit b39ab15e authored by Salson Mikael's avatar Salson Mikael
Browse files

Add an example

parent 39500cb1
No related branches found
No related tags found
No related merge requests found
#include "multicover.hpp"
#include <iostream>
#include <sdsl/rrr_vector.hpp>
#include <sdsl/sd_vector.hpp>
#include <string>
template <class bv>
void print_multicover_info(MultiCover<bv> &cover) {
std::cout << cover.getNbChromosomes() << " chromosomes in this .bam"
<< std::endl << "Number of reads per chromosome" << std::endl;
for (size_t i = 0; i < cover.getNbChromosomes(); i++) {
std::cout << "\t" << cover.getChrName(i)
<< "\t" << cover.getChrCoverage(i).count_reads_between(0, cover.getChrCoverage(i).nb_genome_positions()) << " reads" << std::endl;
}
}
int main(int argc, char **argv) {
if (argc <= 1) {
std::cerr << "Usage: " << argv[0] << " file.bam" << std::endl;
exit(1);
}
MultiCover<sd_vector<> > cover(argv[1]);
print_multicover_info(cover);
// Generate a .bed.gz of the number of reads in every 1,000bp
// non-overlapping range
cover.generateBED(1000, std::string(argv[1])+".bed.gz");
std::cout << "Saving the structure to disk" << std::endl;
cover.save(argv[1]);
std::cout << "Loading the structure that have just been saved" << std::endl;
MultiCover<sd_vector<> > cover2(argv[1], true);
cover2.generateBED(1000, std::string(argv[1])+"-2.bed.gz");
print_multicover_info(cover2);
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment