Add native CRAM output support (unsorted and coordinate-sorted) via HTSlib

source/BAMoutput.cpp

@@ -1,3 +1,4 @@
+
 #include "BAMoutput.h"
 #include <sys/stat.h>
 #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;
@@ -49,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) {
 

source/BAMoutput.h

@@ -5,15 +5,21 @@
 #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);
     void coordOneAlign (char *bamIn, uint bamSize, uint iRead);
     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;
 };

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);

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 <string> bamBinNamesV;
-            for (uint32 ibin=0; ibin<nBins; ibin++) {
-
-                bamBinNamesV.push_back(P.outBAMsortTmpDir+"/b"+std::to_string((uint) ibin));
-                struct stat buffer;
-                if (stat (bamBinNamesV.back().c_str(), &buffer) != 0) {//check if file exists
-                    bamBinNamesV.pop_back();
-                };
-            };
-            for (uint32 ibin=0; ibin<bamBinNamesV.size(); ibin++) {
+            // Output format detection
+            std::string outFormat = P.outSAMtype.empty() ? "BAM" : P.outSAMtype[0];
+            if (outFormat == "CRAM" || outFormat == "SAM" || outFormat == "BAM") {
+                // Open output with hts_open
+                const char* mode = (outFormat == "CRAM") ? "wc" : (outFormat == "SAM" ? "w" : "wb");
+                htsFile* htsOut = hts_open(P.outBAMfileCoordName.c_str(), mode);
+                if (!htsOut) {
+                    ostringstream errOut;
+                    errOut << "EXITING because of fatal ERROR: could not open output file: " << P.outBAMfileCoordName << "\n";
+                    exitWithError(errOut.str(), std::cerr, P.inOut->logMain, 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<nBins; ibin++) {
+                    std::string binFile = P.outBAMsortTmpDir+"/b"+std::to_string((uint) ibin);
+                    struct stat buffer;
+                    if (stat(binFile.c_str(), &buffer) != 0) continue;
+                    BGZF* in = bgzf_open(binFile.c_str(), "r");
+                    if (!in) continue;
+                    bam1_t* b = bam_init1();
+                    while (bam_read1(in, b) >= 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 <string> bamBinNamesV;
+                for (uint32 ibin=0; ibin<nBins; ibin++) {
+                    bamBinNamesV.push_back(P.outBAMsortTmpDir+"/b"+std::to_string((uint) ibin));
+                    struct stat buffer;
+                    if (stat (bamBinNamesV.back().c_str(), &buffer) != 0) {
+                        bamBinNamesV.pop_back();
+                    }
+                }
+                for (uint32 ibin=0; ibin<bamBinNamesV.size(); ibin++) {
                     bamBinNames[ibin] = (char*) bamBinNamesV.at(ibin).c_str();
-            };
-            bam_cat(bamBinNamesV.size(), bamBinNames, 0, P.outBAMfileCoordName.c_str());
+                }
+                bam_cat(bamBinNamesV.size(), bamBinNames, 0, P.outBAMfileCoordName.c_str());
+            }
         };
     };    
 };