Convert the submodule star-sys/STAR to a git subtree

star-sys/STAR/source/BAMbinSortByCoordinate.cpp

@@ -0,0 +1,74 @@
+#include "BAMbinSortByCoordinate.h"
+#include "ErrorWarning.h"
+#include "serviceFuns.cpp"
+#include "BAMfunctions.h"
+
+void BAMbinSortByCoordinate(uint32 iBin, uint binN, uint binS, uint nThreads, string dirBAMsort, Parameters &P, Genome &mapGen) {
+
+    if (binS==0) return; //nothing to do for empty bins
+    //allocate arrays
+    char *bamIn=new char[binS+1];
+    uint *startPos=new uint[binN*3];
+
+    uint bamInBytes=0;
+    //load all aligns
+    for (uint it=0; it<nThreads; it++) {
+        string bamInFile=dirBAMsort+to_string(it)+"/"+to_string((uint) iBin);
+        ifstream bamInStream;
+        bamInStream.open(bamInFile.c_str(),std::ios::binary | std::ios::ate);//open at the end to get file size
+        int64 s1=bamInStream.tellg();
+        if (s1>0)         {
+            bamInStream.seekg(std::ios::beg);
+            bamInStream.read(bamIn+bamInBytes,s1);//read the whole file
+        } else if (s1<0) {
+            ostringstream errOut;
+            errOut << "EXITING because of FATAL ERROR: failed reading from temporary file: " << dirBAMsort+to_string(it)+"/"+to_string((uint) iBin);
+            exitWithError(errOut.str(),std::cerr, P.inOut->logMain, 1, P);
+        };
+        bamInBytes += bamInStream.gcount();
+        bamInStream.close();
+        remove(bamInFile.c_str());
+    };
+    if (bamInBytes!=binS) {
+        ostringstream errOut;
+        errOut << "EXITING because of FATAL ERROR: number of bytes expected from the BAM bin does not agree with the actual size on disk: ";
+        errOut << "Expected bin size=" <<binS <<" ; size on disk="<< bamInBytes <<" ; bin number="<< iBin <<"\n";
+        exitWithError(errOut.str(),std::cerr, P.inOut->logMain, 1, P);
+    };
+
+    //extract coordinates
+
+    for (uint ib=0,ia=0;ia<binN;ia++) {
+        uint32 *bamIn32=(uint32*) (bamIn+ib);
+        startPos[ia*3]  =( ((uint) bamIn32[1]) << 32) | ( (uint)bamIn32[2] );
+        startPos[ia*3+2]=ib;
+        ib+=bamIn32[0]+sizeof(uint32);//note that size of the BAM record does not include the size record itself
+        startPos[ia*3+1]=*( (uint*) (bamIn+ib) ); //read order
+        ib+=sizeof(uint);
+    };
+
+    //sort
+    qsort((void*) startPos, binN, sizeof(uint)*3, funCompareArrays<uint,3>);
+
+    BGZF *bgzfBin;
+    bgzfBin=bgzf_open((dirBAMsort+"/b"+to_string((uint) iBin)).c_str(),("w"+to_string((long long) P.outBAMcompression)).c_str());
+    if (bgzfBin==NULL) {
+        ostringstream errOut;
+        errOut <<"EXITING because of fatal ERROR: could not open temporary bam file: " << dirBAMsort+"/b"+to_string((uint) iBin) << "\n";
+        errOut <<"SOLUTION: check that the disk is not full, increase the max number of open files with Linux command ulimit -n before running STAR";
+        exitWithError(errOut.str(), std::cerr, P.inOut->logMain, EXIT_CODE_PARAMETER, P);
+    };
+
+    outBAMwriteHeader(bgzfBin,P.samHeaderSortedCoord,mapGen.chrNameAll,mapGen.chrLengthAll);
+    //send ordered aligns to bgzf one-by-one
+    for (uint ia=0;ia<binN;ia++) {
+        char* ib=bamIn+startPos[ia*3+2];
+        bgzf_write(bgzfBin,ib, *((uint32*) ib)+sizeof(uint32) );
+    };
+
+    bgzf_flush(bgzfBin);
+    bgzf_close(bgzfBin);
+    //release memory
+    delete [] bamIn;
+    delete [] startPos;
+};

star-sys/STAR/source/BAMbinSortByCoordinate.h

@@ -0,0 +1,11 @@
+#ifndef CODE_BAMbinSortByCoordinate
+#define CODE_BAMbinSortByCoordinate
+#include "IncludeDefine.h"
+#include "Parameters.h"
+#include "Genome.h"
+
+#include SAMTOOLS_BGZF_H
+
+void BAMbinSortByCoordinate(uint32 iBin, uint binN, uint binS, uint nThreads, string dirBAMsort, Parameters &P, Genome &mapGen);
+
+#endif

star-sys/STAR/source/BAMbinSortUnmapped.cpp

@@ -0,0 +1,80 @@
+#include "BAMbinSortUnmapped.h"
+#include "ErrorWarning.h"
+#include "BAMfunctions.h"
+
+void BAMbinSortUnmapped(uint32 iBin, uint nThreads, string dirBAMsort, Parameters &P, Genome &mapGen) {
+
+    BGZF *bgzfBin;
+    bgzfBin=bgzf_open((dirBAMsort+"/b"+to_string((uint) iBin)).c_str(),("w"+to_string((long long) P.outBAMcompression)).c_str());
+    if (bgzfBin==NULL) {
+        ostringstream errOut;
+        errOut <<"EXITING because of fatal ERROR: could not open temporary bam file: " << dirBAMsort+"/b"+to_string((uint) iBin) << "\n";
+        errOut <<"SOLUTION: check that the disk is not full, increase the max number of open files with Linux command ulimit -n before running STAR";
+        exitWithError(errOut.str(), std::cerr, P.inOut->logMain, EXIT_CODE_PARAMETER, P);
+    };
+
+    outBAMwriteHeader(bgzfBin,P.samHeaderSortedCoord,mapGen.chrNameAll,mapGen.chrLengthAll);
+
+
+    vector<string> bamInFile;
+    std::map <uint,uint> startPos;
+
+    for (uint it=0; it<nThreads; it++) {//files from all threads, and BySJout
+        bamInFile.push_back(dirBAMsort+to_string(it)+"/"+to_string((uint) iBin));
+        bamInFile.push_back(dirBAMsort+to_string(it)+"/"+to_string((uint) iBin)+".BySJout");
+    };
+    vector<uint32> bamSize(bamInFile.size(),0);//record sizes
+
+    //allocate arrays
+    char **bamIn=new char* [bamInFile.size()];
+    ifstream *bamInStream = new ifstream [bamInFile.size()];
+
+    for (uint it=0; it<bamInFile.size(); it++) {//initialize
+        bamIn[it] = new char [BAMoutput_oneAlignMaxBytes];
+
+        bamInStream[it].open(bamInFile.at(it).c_str());//opean all files
+
+        bamInStream[it].read(bamIn[it],sizeof(int32));//read BAM record size
+        if (bamInStream[it].good()) {
+            bamSize[it]=((*(uint32*)bamIn[it])+sizeof(int32));//true record size +=4 (4 bytes for uint-iRead)
+            bamInStream[it].read(bamIn[it]+sizeof(int32),bamSize.at(it)-sizeof(int32)+sizeof(uint));//read the rest of the record, including last uint = iRead
+            startPos[*(uint*)(bamIn[it]+bamSize.at(it))]=it;//startPos[iRead]=it : record the order of the files to output
+        } else {//nothing to do here, file is empty, do not record it
+        };
+    };
+
+    //send ordered aligns to bgzf one-by-one
+    while (startPos.size()>0) {
+        uint it=startPos.begin()->second;
+        uint startNext=startPos.size()>1 ? (++startPos.begin())->first : (uint) -1;
+
+        while (true) {
+            bgzf_write(bgzfBin, bamIn[it], bamSize.at(it));
+            bamInStream[it].read(bamIn[it],sizeof(int32));//read record size
+            if (bamInStream[it].good()) {
+                 bamSize[it]=((*(uint32*)bamIn[it])+sizeof(int32));
+                 bamInStream[it].read(bamIn[it]+sizeof(int32),bamSize.at(it)-sizeof(int32)+sizeof(uint));//read the rest of the record, including la$
+                 uint iRead=*(uint*)(bamIn[it]+bamSize.at(it));
+                 if (iRead>startNext) {//this read from this chunk is > than a read from another chunk
+                     startPos[iRead]=it;
+                     break;
+                 };
+            } else {//nothing to do here, reached the end of the file
+                 break;
+            };
+        };
+        startPos.erase(startPos.begin());
+    };
+
+    bgzf_flush(bgzfBin);
+    bgzf_close(bgzfBin);
+
+
+    for (uint it=0; it<bamInFile.size(); it++) {//destroy at the end
+        bamInStream[it].close();
+        remove(bamInFile.at(it).c_str());
+        delete [] bamIn[it];
+    };
+    delete [] bamIn;
+    delete [] bamInStream;
+};

star-sys/STAR/source/BAMbinSortUnmapped.h

@@ -0,0 +1,11 @@
+#ifndef CODE_BAMbinSortUnmapped
+#define CODE_BAMbinSortUnmapped
+#include "IncludeDefine.h"
+#include "Parameters.h"
+#include "Genome.h"
+
+#include SAMTOOLS_BGZF_H
+
+void BAMbinSortUnmapped(uint32 iBin, uint nThreads, string dirBAMsort, Parameters &P, Genome &mapGen);
+
+#endif

star-sys/STAR/source/BAMfunctions.cpp

@@ -0,0 +1,112 @@
+#include <stdexcept>
+#include "BAMfunctions.h"
+#include "htslib/htslib/kstring.h"
+
+string bam_cigarString (bam1_t *b) {//output CIGAR string
+//    kstring_t strK;
+//    kstring_t *str=&strK;
+    const bam1_core_t *c = &b->core;
+
+    string cigarString("");
+    if ( c->n_cigar > 0 ) {
+      uint32_t *cigar = bam_get_cigar(b);
+      for (int i = 0; i < c->n_cigar; ++i) {
+        cigarString+=to_string((uint)bam_cigar_oplen(cigar[i]))+bam_cigar_opchr(cigar[i]);
+      };
+    };
+
+
+//	if (c->n_cigar) { // cigar
+//		for (int i = 0; i < c->n_cigar; ++i) {
+//			kputw(bam_cigar_oplen(cigar[i]), str);
+//			kputc(bam_cigar_opchr(cigar[i]), str);
+//		}
+//	} else kputc('*', str);
+//
+//    string cigarString (str->s,str->l);
+    return cigarString;
+};
+
+int bam_read1_fromArray(char *bamChar, bam1_t *b) //modified from samtools bam_read1 to assign BAM record in mmemry to bam structure
+{
+	bam1_core_t *c = &b->core;
+	int32_t block_len; //, ret, i;
+// // 	uint32_t x[8];
+// // 	if ((ret = bgzf_read(fp, &block_len, 4)) != 4) {
+// // 		if (ret == 0) return -1; // normal end-of-file
+// // 		else return -2; // truncated
+// // 	}
+ 	uint32_t *x;
+
+    uint32_t *bamU32=(uint32_t*) bamChar;
+    block_len=bamU32[0];
+
+// // 	if (bgzf_read(fp, x, 32) != 32) return -3;
+// // 	if (fp->is_be) {
+// // 		ed_swap_4p(&block_len);
+// // 		for (i = 0; i < 8; ++i) ed_swap_4p(x + i);
+// // 	}
+    x=bamU32+1;
+
+	c->tid = x[0]; c->pos = x[1];
+	c->bin = x[2]>>16; c->qual = x[2]>>8&0xff; c->l_qname = x[2]&0xff;
+	c->flag = x[3]>>16; c->n_cigar = x[3]&0xffff;
+	c->l_qseq = x[4];
+	c->mtid = x[5]; c->mpos = x[6]; c->isize = x[7];
+	b->l_data = block_len - 32;
+	if (b->l_data < 0 || c->l_qseq < 0) return -4;
+	if ((char *)bam_get_aux(b) - (char *)b->data > b->l_data)
+		return -4;
+	if (b->m_data < b->l_data) {
+		b->m_data = b->l_data;
+		kroundup32(b->m_data);
+		b->data = (uint8_t*)realloc(b->data, b->m_data);
+		if (!b->data)
+			return -4;
+	}
+// // 	if (bgzf_read(fp, b->data, b->l_data) != b->l_data) return -4;
+// // 	//b->l_aux = b->l_data - c->n_cigar * 4 - c->l_qname - c->l_qseq - (c->l_qseq+1)/2;
+// // 	if (fp->is_be) swap_data(c, b->l_data, b->data, 0);
+    b->data=(uint8_t*) bamChar+4*9;
+
+	return 4 + block_len;
+}
+
+
+void outBAMwriteHeader (BGZF* fp, const string &samh, const vector <string> &chrn, const vector <uint> &chrl) {
+    throw std::runtime_error("Unimplemented!");
+    //bgzf_write(fp,"BAM\001",4);
+    int32 hlen=samh.size();
+    //bgzf_write(fp,(char*) &hlen,sizeof(hlen));
+    //bgzf_write(fp,samh.c_str(),hlen);
+    int32 nchr=(int32) chrn.size();
+    //bgzf_write(fp,(char*) &nchr,sizeof(nchr));
+    for (int32 ii=0;ii<nchr;ii++) {
+        int32 rlen = (int32) (chrn.at(ii).size()+1);
+        int32 slen = (int32) chrl[ii];
+        //bgzf_write(fp,(char*) &rlen,sizeof(rlen));
+        //bgzf_write(fp,chrn.at(ii).data(),rlen); //this includes \0 at the end of the string
+        //bgzf_write(fp,(char*) &slen,sizeof(slen));
+    };
+    //bgzf_flush(fp);
+};
+
+template <class TintType>
+TintType bamAttributeInt(const char *bamAux, const char *attrName) {//not tested!!!
+    const char *attrStart=strstr(bamAux,attrName);
+    if (attrStart==NULL) return (TintType) -1;
+    switch (attrStart[2]) {
+        case ('c'):
+            return (TintType) *(int8_t*)(attrStart+3);
+        case ('s'):
+            return (TintType) *(int16_t*)(attrStart+3);
+        case ('i'):
+            return (TintType) *(int32_t*)(attrStart+3);
+        case ('C'):
+            return (TintType) *(uint8_t*)(attrStart+3);
+        case ('S'):
+            return (TintType) *(uint16_t*)(attrStart+3);
+        case ('I'):
+            return (TintType) *(uint32_t*)(attrStart+3);
+    };
+};

star-sys/STAR/source/BAMfunctions.h

@@ -0,0 +1,10 @@
+#ifndef DEF_BAMfunctions
+#define DEF_BAMfunctions
+
+#include "IncludeDefine.h"
+#include SAMTOOLS_BGZF_H
+#include SAMTOOLS_SAM_H
+        void outBAMwriteHeader (BGZF* fp, const string &samh, const vector <string> &chrn, const vector <uint> &chrl);
+        int bam_read1_fromArray(char *bamChar, bam1_t *b);
+        string bam_cigarString (bam1_t *b);
+#endif