Attempt to allow offline CRAM

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

common.h

@@ -13,7 +13,7 @@
 #include <vector>
 
 // global definitions/constants
-#define VERSION "0.0.2"
+#define VERSION "0.0.3"
 
 // definitions/constants for argparsing
 #define DEFAULT_NUM_THREADS 1

count.cpp

@@ -58,7 +58,7 @@ void write_ins_counts_json(std::unordered_map<uint32_t, std::unordered_map<std::
     }
 }
 
-counts_t compute_counts(const char* const in_reads_fn, std::string const & ref, uint8_t const min_qual, std::vector<std::pair<uint32_t, uint32_t>> 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<std::pair<uint32_t, uint32_t>> 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()

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<std::pair<uint32_t, uint32_t>> 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<std::pair<uint32_t, uint32_t>> const & min_max_primer_inds, args_t const & user_args);
 
 // compute consensus genome sequence from counts
 std::string compute_consensus(std::vector<std::array<COUNT_T, 5>> const & pos_counts, std::unordered_map<uint32_t, std::unordered_map<std::string, COUNT_T>> & ins_counts, args_t const & user_args);

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<std::pair<uint32_t, uint32_t>> 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);