diff --git a/CMakeLists.txt b/CMakeLists.txt index e97263a6..f8983067 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -23,7 +23,7 @@ add_compile_definitions("MAX_KMER_SIZE=${MAX_KMER_SIZE}") option(USE_HDF5 "Compile with HDF5 support" OFF) #OFF by default -option(USE_BAM "Compile with HTSLIB support" OFF) # OFF by default +option(USE_BAM "Compile with HTSLIB support" ON) # ON by défaut if(USE_HDF5) add_compile_definitions("USE_HDF5=ON") diff --git a/src/ProcessReads.cpp b/src/ProcessReads.cpp index 1d1e4fff..d7ad7f01 100755 --- a/src/ProcessReads.cpp +++ b/src/ProcessReads.cpp @@ -649,24 +649,52 @@ void MasterProcessor::update(const std::vector& c, const std::vector 0) { + bamfn = opt.output + "/pseudoalignments.cram"; + is_cram = true; + } else if (bamfn.size() >= 5 && bamfn.substr(bamfn.size()-5) == ".cram") { + is_cram = true; + } + if (opt.pseudobam) { if (opt.genomebam) { bamh = createPseudoBamHeaderGenome(model); - bamfps = new htsFile*[numSortFiles]; for (int i = 0; i < numSortFiles; i++) { - bamfps[i] = sam_open((opt.output + "/tmp." + std::to_string(i) + ".bam").c_str(), "wb1"); + std::string tmpFile = opt.output + "/tmp." + std::to_string(i) + (is_cram ? ".cram" : ".bam"); + bamfps[i] = sam_open(tmpFile.c_str(), is_cram ? "wc" : "wb1"); int r = sam_hdr_write(bamfps[i], bamh); + if (is_cram && opt.reference_fasta.size() > 0) { + hts_set_fai_filename(bamfps[i], opt.reference_fasta.c_str()); + } } - } else { bamh = createPseudoBamHeaderTrans(index); - bamfp = sam_open(bamfn.c_str(), "wb"); + bamfp = sam_open(bamfn.c_str(), is_cram ? "wc" : "wb"); int r = sam_hdr_write(bamfp, bamh); + if (is_cram && opt.reference_fasta.size() > 0) { + hts_set_fai_filename(bamfp, opt.reference_fasta.c_str()); + } } - if (opt.threads > 1 && !opt.genomebam) { - // makes no sens to use threads on unsorted bams hts_set_threads(bamfp, opt.threads); } } diff --git a/src/common.h b/src/common.h index 3709021f..37e24d7a 100755 --- a/src/common.h +++ b/src/common.h @@ -91,6 +91,9 @@ struct BUSOptions { }; struct ProgramOptions { + bool force_bam = false; // --bam + bool force_cram = false; // --cram + std::string bam_or_cram; // "bam" or "cram" to force output format bool verbose; bool aa; bool distinguish; @@ -159,6 +162,7 @@ struct ProgramOptions { std::string genemap; std::string priors; std::string tmp_dir; + std::string reference_fasta; // chemin vers la référence pour CRAM ProgramOptions() : verbose(false), diff --git a/src/main.cpp b/src/main.cpp index 3d1097fb..8b2c0fc1 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -415,7 +415,8 @@ void ParseOptionsTCCQuant(int argc, char **argv, ProgramOptions& opt) { {"gtf", required_argument, 0, 'G'}, {"bootstrap-samples", required_argument, 0, 'b'}, {"seed", required_argument, 0, 'd'}, - {"priors", required_argument, 0, 'p'}, + {"priors", required_argument, 0, 'p'}, + {"cram-reference", required_argument, 0, 1001}, // nouvelle option pour CRAM {0,0,0,0} }; int c; @@ -427,7 +428,7 @@ void ParseOptionsTCCQuant(int argc, char **argv, ProgramOptions& opt) { break; } - switch (c) { + switch (c) { case 0: break; case 't': { @@ -487,6 +488,10 @@ void ParseOptionsTCCQuant(int argc, char **argv, ProgramOptions& opt) { opt.priors = optarg; break; } + case 1001: { + opt.reference_fasta = optarg; + break; + } default: break; } } @@ -2130,7 +2135,8 @@ void usageBus() { << " --aa Align to index generated from a FASTA-file containing amino acid sequences" << endl << " --inleaved Specifies that input is an interleaved FASTQ file" << endl << " --batch-barcodes Records both batch and extracted barcode in BUS file" << endl - << " --verbose Print out progress information every 1M proccessed reads" << endl; + << " --verbose Print out progress information every 1M proccessed reads" << endl + << " --cram-reference=FASTA Output CRAM using the provided FASTA reference (must be the same as for the index, enables .cram output)" << endl; } void usageIndex() {