From 3e90996150e4dcb27846cf6b45946616096b8697 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Mon, 27 Oct 2025 21:42:54 +0100 Subject: [PATCH 1/4] feat(cram): initial CRAM output support in BAMoutput (htsFile constructor and header) --- source/BAMoutput.cpp | 17 +++++++++++++++++ source/BAMoutput.h | 13 ++++++++++--- 2 files changed, 27 insertions(+), 3 deletions(-) diff --git a/source/BAMoutput.cpp b/source/BAMoutput.cpp index e61faa64..074d11cc 100644 --- a/source/BAMoutput.cpp +++ b/source/BAMoutput.cpp @@ -1,3 +1,4 @@ + #include "BAMoutput.h" #include #include "GlobalVariables.h" @@ -6,6 +7,22 @@ #include "ThreadControl.h" #include "streamFuns.h" +BAMoutput::BAMoutput(htsFile *htsOutIn, Parameters &Pin) : P(Pin) { + // allocate BAM array with one bin, streamed directly into htsFile (SAM/BAM/CRAM) + bamArraySize = P.chunkOutBAMsizeBytes; + bamArray = new char[bamArraySize]; + binBytes1 = 0; + htsOut = htsOutIn; + bgzfBAM = nullptr; + binSize = 0; + binStream = nullptr; + binStart = nullptr; + binBytes = nullptr; + binTotalBytes = nullptr; + binTotalN = nullptr; + nBins = 0; +} + BAMoutput::BAMoutput (int iChunk, string tmpDir, Parameters &Pin) : P(Pin){//allocate bam array nBins=P.outBAMcoordNbins; diff --git a/source/BAMoutput.h b/source/BAMoutput.h index bc939eac..57e72abb 100644 --- a/source/BAMoutput.h +++ b/source/BAMoutput.h @@ -5,7 +5,12 @@ #include SAMTOOLS_BGZF_H #include "Parameters.h" -class BAMoutput {// +extern "C" { +#include "htslib/htslib/hts.h" +#include "htslib/htslib/sam.h" +} + +class BAMoutput { public: //sorted output BAMoutput (int iChunk, string tmpDir, Parameters &Pin); @@ -13,7 +18,8 @@ class BAMoutput {// void coordBins (); void coordFlush (); //unsorted output - BAMoutput (BGZF *bgzfBAMin, Parameters &Pin); + BAMoutput (BGZF *bgzfBAMin, Parameters &Pin); // legacy BAM + BAMoutput (htsFile *htsOut, Parameters &Pin); // new: generic (SAM/BAM/CRAM) void unsortedOneAlign (char *bamIn, uint bamSize, uint bamSize2); void unsortedFlush (); void coordUnmappedPrepareBySJout(); @@ -29,7 +35,8 @@ class BAMoutput {// char **binStart; //pointers to starts of the bins uint64 *binBytes, binBytes1;//number of bytes currently written to each bin ofstream **binStream;//output streams for each bin - BGZF *bgzfBAM; + BGZF *bgzfBAM; // legacy BAM + htsFile *htsOut; // new: generic (SAM/BAM/CRAM) Parameters &P; string bamDir; }; From 284a8414ee76387ad8f333b97ba0c7ebc2ef9042 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Mon, 27 Oct 2025 21:44:25 +0100 Subject: [PATCH 2/4] feat(cram): detect output format and open with hts_open for CRAM/SAM/BAM output --- source/ReadAlignChunk.cpp | 49 ++++++++++++++++++++++++++++----------- 1 file changed, 35 insertions(+), 14 deletions(-) diff --git a/source/ReadAlignChunk.cpp b/source/ReadAlignChunk.cpp index 8f7f180b..c5b319da 100644 --- a/source/ReadAlignChunk.cpp +++ b/source/ReadAlignChunk.cpp @@ -39,29 +39,50 @@ ReadAlignChunk::ReadAlignChunk(Parameters& Pin, Genome &genomeIn, Transcriptome chunkOutBAMtotal=0; }; + // Detect output format and open with hts_open if needed + auto get_hts_mode = [](const std::string& format) -> const char* { + if (format == "CRAM") return "wc"; + if (format == "BAM") return "wb"; + if (format == "SAM") return "w"; + return nullptr; + }; + + std::string outFormat = P.outSAMtype.empty() ? "BAM" : P.outSAMtype[0]; + if (P.outBAMunsorted) { - chunkOutBAMunsorted = new BAMoutput (P.inOut->outBAMfileUnsorted, P); + if (outFormat == "CRAM" || outFormat == "SAM" || outFormat == "BAM") { + htsFile* htsOut = hts_open(P.outBAMfileUnsortedName.c_str(), get_hts_mode(outFormat)); + chunkOutBAMunsorted = new BAMoutput(htsOut, P); + } else { + chunkOutBAMunsorted = new BAMoutput(P.inOut->outBAMfileUnsorted, P); + } RA->outBAMunsorted = chunkOutBAMunsorted; } else { - chunkOutBAMunsorted=NULL; - RA->outBAMunsorted=NULL; - }; + chunkOutBAMunsorted = NULL; + RA->outBAMunsorted = NULL; + } if (P.outBAMcoord) { - chunkOutBAMcoord = new BAMoutput (iChunk, P.outBAMsortTmpDir, P); + // Coordinate sorted output is handled via temp bins, keep as is for now + chunkOutBAMcoord = new BAMoutput(iChunk, P.outBAMsortTmpDir, P); RA->outBAMcoord = chunkOutBAMcoord; } else { - chunkOutBAMcoord=NULL; - RA->outBAMcoord=NULL; - }; - - if ( P.quant.trSAM.bamYes ) { - chunkOutBAMquant = new BAMoutput (P.inOut->outQuantBAMfile,P); + chunkOutBAMcoord = NULL; + RA->outBAMcoord = NULL; + } + + if (P.quant.trSAM.bamYes) { + if (outFormat == "CRAM" || outFormat == "SAM" || outFormat == "BAM") { + htsFile* htsOut = hts_open(P.outQuantBAMfileName.c_str(), get_hts_mode(outFormat)); + chunkOutBAMquant = new BAMoutput(htsOut, P); + } else { + chunkOutBAMquant = new BAMoutput(P.inOut->outQuantBAMfile, P); + } RA->outBAMquant = chunkOutBAMquant; } else { - chunkOutBAMquant=NULL; - RA->outBAMquant=NULL; - }; + chunkOutBAMquant = NULL; + RA->outBAMquant = NULL; + } if (P.outSJ.yes) { chunkOutSJ = new OutSJ (P.limitOutSJcollapsed, P, mapGen); From 7e6641c7270054abc3715225f4d439c581db43e2 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Mon, 27 Oct 2025 21:44:47 +0100 Subject: [PATCH 3/4] feat(cram): enable sam_write1 for htsFile (CRAM/SAM/BAM) output in BAMoutput --- source/BAMoutput.cpp | 37 ++++++++++++++++++++++++------------- 1 file changed, 24 insertions(+), 13 deletions(-) diff --git a/source/BAMoutput.cpp b/source/BAMoutput.cpp index 074d11cc..79eea1be 100644 --- a/source/BAMoutput.cpp +++ b/source/BAMoutput.cpp @@ -66,30 +66,41 @@ BAMoutput::BAMoutput (BGZF *bgzfBAMin, Parameters &Pin) : P(Pin){//allocate BAM nBins=0; }; -void BAMoutput::unsortedOneAlign (char *bamIn, uint bamSize, uint bamSize2) {//record one alignment to the buffer, write buffer if needed - - if (bamSize==0) return; //no output, could happen if one of the mates is not mapped - - if (binBytes1+bamSize2 > bamArraySize) {//write out this buffer +void BAMoutput::unsortedOneAlign (char *bamIn, uint bamSize, uint bamSize2) { + if (bamSize==0) return; + if (htsOut) { + // Write directly using sam_write1 for each record + bam1_t b; + memset(&b, 0, sizeof(bam1_t)); + bam_read1_fromArray(bamIn, &b); + sam_write1(htsOut, NULL, &b); // TODO: pass correct header if available + if (b.data) b.data = NULL; // avoid double free + return; + } + // Legacy BGZF mode + if (binBytes1+bamSize2 > bamArraySize) { if (g_threadChunks.threadBool) pthread_mutex_lock(&g_threadChunks.mutexOutSAM); bgzf_write(bgzfBAM,bamArray,binBytes1); if (g_threadChunks.threadBool) pthread_mutex_unlock(&g_threadChunks.mutexOutSAM); - - binBytes1=0;//rewind the buffer - }; - + binBytes1=0; + } memcpy(bamArray+binBytes1, bamIn, bamSize); binBytes1 += bamSize; +} -}; -void BAMoutput::unsortedFlush () {//flush all alignments +void BAMoutput::unsortedFlush () { + if (htsOut) { + // nothing to flush for htsFile, records are written immediately + binBytes1=0; + return; + } if (g_threadChunks.threadBool) pthread_mutex_lock(&g_threadChunks.mutexOutSAM); bgzf_write(bgzfBAM,bamArray,binBytes1); if (g_threadChunks.threadBool) pthread_mutex_unlock(&g_threadChunks.mutexOutSAM); - binBytes1=0;//rewind the buffer -}; + binBytes1=0; +} void BAMoutput::coordOneAlign (char *bamIn, uint bamSize, uint iRead) { From 979423910a03015b845221e25e89669d448edb2b Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Mon, 27 Oct 2025 21:45:52 +0100 Subject: [PATCH 4/4] feat(cram): enable coordinate-sorted CRAM/SAM/BAM output using sam_write1 after sorting --- source/bamSortByCoordinate.cpp | 59 ++++++++++++++++++++++++++-------- 1 file changed, 45 insertions(+), 14 deletions(-) diff --git a/source/bamSortByCoordinate.cpp b/source/bamSortByCoordinate.cpp index 4d7ef11e..b1b38ab5 100644 --- a/source/bamSortByCoordinate.cpp +++ b/source/bamSortByCoordinate.cpp @@ -77,21 +77,52 @@ void bamSortByCoordinate (Parameters &P, ReadAlignChunk **RAchunk, Genome &genom }; }; - //concatenate all BAM files, using bam_cat - char **bamBinNames = new char* [nBins]; - vector bamBinNamesV; - for (uint32 ibin=0; ibinlogMain, EXIT_CODE_PARAMETER, P); + } + // Write header + bam_hdr_t* header = sam_hdr_parse(P.samHeaderSortedCoord.size(), P.samHeaderSortedCoord.c_str()); + sam_hdr_write(htsOut, header); + // For each bin, read and write records + for (uint32 ibin=0; ibin= 0) { + sam_write1(htsOut, header, b); + } + bam_destroy1(b); + bgzf_close(in); + } + sam_hdr_destroy(header); + hts_close(htsOut); + } else { + // Default: BAM output using bam_cat + char **bamBinNames = new char* [nBins]; + vector bamBinNamesV; + for (uint32 ibin=0; ibin