diff --git a/Makefile b/Makefile index e051e49..0d2ea95 100644 --- a/Makefile +++ b/Makefile @@ -1,6 +1,6 @@ # use g++ compiler CXX=g++ -CXXFLAGS?=-Wall -pedantic -std=c++11 +CXXFLAGS?=-std=c++11 -fpermissive # flag specifications for release and debug RELEASEFLAGS?=$(CXXFLAGS) -O3 diff --git a/common.h b/common.h index 50081f9..d8beb41 100644 --- a/common.h +++ b/common.h @@ -13,7 +13,7 @@ #include // global definitions/constants -#define VERSION "0.0.2" +#define VERSION "0.0.3" // definitions/constants for argparsing #define DEFAULT_NUM_THREADS 1 diff --git a/count.cpp b/count.cpp index 8329628..11b5808 100644 --- a/count.cpp +++ b/count.cpp @@ -58,7 +58,7 @@ void write_ins_counts_json(std::unordered_map> const & min_max_primer_inds) { +counts_t compute_counts(const char* const in_reads_fn, std::string const & ref, uint8_t const min_qual, std::vector> const & min_max_primer_inds, args_t const & user_args) { // open reference FASTA file and CRAM/BAM/SAM file htsFile* reads = hts_open(in_reads_fn, "r"); if(!reads) { @@ -75,6 +75,11 @@ counts_t compute_counts(const char* const in_reads_fn, std::string const & ref, std::cerr << "CRAM/BAM/SAM has " << header->n_targets << " references, but it should have exactly 1: " << in_reads_fn << std::endl; exit(1); } + // if CRAM file, load reference in htslib + if(reads->format.format == cram) { + cram_load_reference((reads->fp).cram, user_args.in_ref_fn); + } + // prepare helper variables for computing counts counts_t counts; // counts_t object to store the counts themselves bam1_t* src = bam_init1(); // holds the current alignment record, which is read by sam_read1() diff --git a/count.h b/count.h index 6d7041f..8952d0b 100644 --- a/count.h +++ b/count.h @@ -1,5 +1,6 @@ #ifndef COUNT_H #define COUNT_H +#include "htslib/cram/cram.h" #include "common.h" #include "argparse.h" @@ -10,7 +11,7 @@ struct counts_t { }; // compute position and insertion counts -counts_t compute_counts(const char* const in_reads_fn, std::string const & ref, uint8_t const min_qual, std::vector> const & min_max_primer_inds); +counts_t compute_counts(const char* const in_reads_fn, std::string const & ref, uint8_t const min_qual, std::vector> const & min_max_primer_inds, args_t const & user_args); // compute consensus genome sequence from counts std::string compute_consensus(std::vector> const & pos_counts, std::unordered_map> & ins_counts, args_t const & user_args); diff --git a/main.cpp b/main.cpp index ecf8c65..a6d391c 100644 --- a/main.cpp +++ b/main.cpp @@ -8,7 +8,7 @@ int main(int argc, char** argv) { args_t user_args = parse_args(argc, argv); std::string ref = read_fasta(user_args.in_ref_fn); std::vector> min_max_primer_inds = find_overlapping_primers(ref.length(), user_args.primer_bed_fn, user_args.primer_offset); - counts_t counts = compute_counts(user_args.in_reads_fn, ref, user_args.min_qual, min_max_primer_inds); + counts_t counts = compute_counts(user_args.in_reads_fn, ref, user_args.min_qual, min_max_primer_inds, user_args); std::string consensus = compute_consensus(counts.pos_counts, counts.ins_counts, user_args); write_consensus_fasta(consensus, user_args); write_ins_counts_json(counts.ins_counts, user_args.out_ins_counts_fn);