Skip to content

Add competitive --bam and --cram CLI options, CRAM output with reference, fallback to BAM if reference missing, update help, enable BAM/CRAM by default #499

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
BenjaminDEMAILLE wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Add competitive --bam and --cram CLI options, CRAM output with reference, fallback to BAM if reference missing, update help, enable BAM/CRAM by default #499

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down
40 changes: 34 additions & 6 deletions src/ProcessReads.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -649,24 +649,52 @@ void MasterProcessor::update(const std::vector<uint32_t>& c, const std::vector<R
void MasterProcessor::processAln(const EMAlgorithm& em, bool useEM = true) {
// open bamfile and fetch header
std::string bamfn = opt.output + "/pseudoalignments.bam";
bool is_cram = false;
// Mutually exclusive BAM/CRAM flags
if (opt.force_bam && opt.force_cram) {
std::cerr << "Error: --bam and --cram cannot be used together." << std::endl;
exit(1);
}
if (opt.force_cram) {
if (opt.reference_fasta.empty()) {
std::cerr << "[kallisto] Warning: --cram specified but --cram-reference missing. Falling back to BAM output." << std::endl;
bamfn = opt.output + "/pseudoalignments.bam";
is_cram = false;
} else {
bamfn = opt.output + "/pseudoalignments.cram";
is_cram = true;
}
} else if (opt.force_bam) {
bamfn = opt.output + "/pseudoalignments.bam";
is_cram = false;
} else if (opt.reference_fasta.size() > 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);
}
}
Expand Down
4 changes: 4 additions & 0 deletions src/common.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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),
Expand Down
12 changes: 9 additions & 3 deletions src/main.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -427,7 +428,7 @@ void ParseOptionsTCCQuant(int argc, char **argv, ProgramOptions& opt) {
break;
}

switch (c) {
switch (c) {
case 0:
break;
case 't': {
Expand Down Expand Up @@ -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;
}
}
Expand Down Expand Up @@ -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() {
Expand Down