Fixed a bug causing rare seg-faults with for --peOverlap* options and chimeric detection.

Author: alexdobin

CHANGES.md

@@ -1,3 +1,6 @@
+* Fixed a bug causing rare seg-faults with for --peOverlap* options and chimeric detection.
+* Implemented SAM attributes sS and sQ to output sequence and quality of the entire barcode read, and sM to output the barcode match status.
+* Implemented complex barcodes in STARsolo with --soloType CB_UMI_Complex, --soloCBmatchWLtype --soloAdapterSequence, --soloAdapterMismatchesNmax, --soloCBposition,--soloUMIposition
 * Fixed a problem in STARsolo with unmapped reads counts in Solo.out/*.stats files.
 * Fixed a bug in STARsolo with counting reads for splice junctions. Solo.out/matrixSJ.mtx output is slighlty changed.
 * Implemented CB/UB/GX/GN BAM tags for STARsolo.

bin/Linux_x86_64/STAR

@@ -4,7 +4,7 @@
 #include "BAMfunctions.h"
 #include "SequenceFuns.h"
 
-void BAMbinSortByCoordinate(uint32 iBin, uint binN, uint binS, uint nThreads, string dirBAMsort, Parameters &P, Genome &mapGen, SoloFeature &soloFeat) {
+void BAMbinSortByCoordinate(uint32 iBin, uint binN, uint binS, uint nThreads, string dirBAMsort, Parameters &P, Genome &mapGen, Solo &solo) {
 
     if (binS==0) return; //nothing to do for empty bins
     //allocate arrays
@@ -66,7 +66,10 @@ void BAMbinSortByCoordinate(uint32 iBin, uint binN, uint binS, uint nThreads, st
     for (uint ia=0;ia<binN;ia++) {
         char* bam0=bamIn+startPos[ia*3+2];
         uint32 size0=*((uint32*) bam0)+sizeof(uint32);
-        soloFeat.addBAMtags(bam0,size0,bam1);
+        
+        if (solo.pSolo.samAttrYes)
+            solo.soloFeat[solo.pSolo.samAttrFeature]->addBAMtags(bam0,size0,bam1);
+        
         bgzf_write(bgzfBin, bam0, size0);
     };
 

bin/Linux_x86_64/STARlong

@@ -3,10 +3,10 @@
 #include "IncludeDefine.h"
 #include "Parameters.h"
 #include "Genome.h"
-#include "SoloFeature.h"
+#include "Solo.h"
 
 #include SAMTOOLS_BGZF_H
 
-void BAMbinSortByCoordinate(uint32 iBin, uint binN, uint binS, uint nThreads, string dirBAMsort, Parameters &P, Genome &mapGen, SoloFeature &soloFeat);
+void BAMbinSortByCoordinate(uint32 iBin, uint binN, uint binS, uint nThreads, string dirBAMsort, Parameters &P, Genome &mapGen, Solo &solo);
 
 #endif

bin/Linux_x86_64_static/STAR

@@ -2,7 +2,7 @@
 #include "ErrorWarning.h"
 #include "BAMfunctions.h"
 
-void BAMbinSortUnmapped(uint32 iBin, uint nThreads, string dirBAMsort, Parameters &P, Genome &mapGen, SoloFeature &soloFeat) {
+void BAMbinSortUnmapped(uint32 iBin, uint nThreads, string dirBAMsort, Parameters &P, Genome &mapGen, Solo &solo) {
 
     BGZF *bgzfBin;
     bgzfBin=bgzf_open((dirBAMsort+"/b"+to_string((uint) iBin)).c_str(),("w"+to_string((long long) P.outBAMcompression)).c_str());
@@ -53,7 +53,9 @@ void BAMbinSortUnmapped(uint32 iBin, uint nThreads, string dirBAMsort, Parameter
             //add extra tags to the BAM record
             char* bam0=bamIn[it];
             uint32 size0=bamSize.at(it);
-            soloFeat.addBAMtags(bam0,size0,bam1);
+            
+            if (solo.pSolo.samAttrYes)
+                solo.soloFeat[solo.pSolo.samAttrFeature]->addBAMtags(bam0,size0,bam1);            
 
             bgzf_write(bgzfBin, bam0, size0);
             bamInStream[it].read(bamIn[it],sizeof(int32));//read record size

bin/Linux_x86_64_static/STARlong

@@ -3,10 +3,10 @@
 #include "IncludeDefine.h"
 #include "Parameters.h"
 #include "Genome.h"
-#include "SoloFeature.h"
+#include "Solo.h"
 
 #include SAMTOOLS_BGZF_H
 
-void BAMbinSortUnmapped(uint32 iBin, uint nThreads, string dirBAMsort, Parameters &P, Genome &mapGen, SoloFeature &soloFeat);
+void BAMbinSortUnmapped(uint32 iBin, uint nThreads, string dirBAMsort, Parameters &P, Genome &mapGen, Solo &solo);
 
 #endif

source/BAMbinSortByCoordinate.cpp

@@ -13,8 +13,11 @@ void ParametersSolo::initialize(Parameters *pPin)
 {
     pP=pPin;
 
+    yes = true;
     if (typeStr=="None") {
-        type=0;
+        type = 0;
+        yes = false;
+        samAttrYes = false;
         //solo SAM attributes not allowed
         if (pP->outSAMattrPresent.CR || pP->outSAMattrPresent.CY || pP->outSAMattrPresent.UR || pP->outSAMattrPresent.UY  || pP->outSAMattrPresent.CB  || pP->outSAMattrPresent.UB) {
             ostringstream errOut;

source/BAMbinSortByCoordinate.h

@@ -10,11 +10,10 @@ class ParametersSolo {
     //chemistry, library etc
     string typeStr;
     int type;
+    bool yes;
     string strandStr;
     int32 strand;   
 
-
-    
     //simple barcodes
     uint32 cbS,cbL; //cell barcode start,length
     uint32 umiS,umiL; //umi start,length

source/BAMbinSortUnmapped.cpp

@@ -95,7 +95,10 @@ void ReadAlign::copyRead(ReadAlign &r) {//copy read information only
     readLengthPairOriginal=r.readLengthPairOriginal;
     outFilterMismatchNmaxTotal=r.outFilterMismatchNmaxTotal;
     readName=r.readName;
-
+    iReadAll=r.iReadAll;
+    readFilter=r.readFilter;
+    readFilesIndex=r.readFilesIndex;
+    
     for (uint ii=0;ii<=2;ii++)
         memcpy(Read1[ii],r.Read1[ii],Lread);//need to copy since it will be changed
     Qual1=r.Qual1;

source/BAMbinSortUnmapped.h

@@ -383,7 +383,7 @@ int main(int argInN, char* argIn[]) {
         RAchunk[0]->chunkFilesCat(P.inOut->outSAM, P.outFileTmp + "/Aligned.out.sam.chunk", g_threadChunks.chunkOutN);
     };
 
-    bamSortByCoordinate(P, RAchunk, mainGenome, *soloMain.soloFeat[soloMain.pSolo.samAttrFeature]);
+    bamSortByCoordinate(P, RAchunk, mainGenome, soloMain);
 
     //wiggle output
     if (P.outWigFlags.yes) {

source/ParametersSolo.cpp

@@ -4,9 +4,6 @@
 
 void SoloFeature::addBAMtags(char *&bam0, uint32 &size0, char *bam1)
 {//add extra tags to the BAM record
-
-    if (!pSolo.samAttrYes)
-            return;
 
     uint64 iread = * ((uint64*) (bam0+size0));
     iread = iread >> 32; //iRead was encoded in the upper 32 bitsls

source/ParametersSolo.h

@@ -129,10 +129,12 @@ bool SoloReadBarcode::convertCheckUMI()
 {//check UMIs, return if bad UMIs
     if (convertNuclStrToInt32(umiSeq,umiB)!=-1) {//convert and check for Ns
         stats.V[stats.nNinUMI]++;//UMIs are not allowed to have Ns
+        cbMatch=-23;
         return false;
     };
     if (umiB==homoPolymer[0] || umiB==homoPolymer[1] || umiB==homoPolymer[2] || umiB==homoPolymer[3]) {
         stats.V[stats.nUMIhomopolymer]++;
+        cbMatch=-24;        
         return false;
     };
     return true;
@@ -167,6 +169,8 @@ void SoloReadBarcode::getCBandUMI(string &readNameExtra)
 
         cbSeq="";
         cbQual="";
+        umiSeq="";
+        umiQual="";
 
         uint32 adapterStart=0;
         if (pSolo.adapterYes) {
@@ -183,19 +187,21 @@ void SoloReadBarcode::getCBandUMI(string &readNameExtra)
             return;
         };
 
-        if (!convertCheckUMI())
-            return;        
-        
         bool cbMatchGood=true;
+        if (!convertCheckUMI())
+            cbMatchGood=false;//CB matching will not be done, just extract the sequences
+
         cbMatchInd={0};
         for (auto &cb : pSolo.cbV) {//cycle over multiple barcodes
 
             string cbSeq1, cbQual1;
-            if (!cb.extractBarcode(bSeq, bQual, adapterStart, cbSeq1, cbQual1) 
-                || cbSeq1.size() < cb.minLen || cbSeq1.size() >= cb.wl.size() || cb.wl[cbSeq1.size()].size()==0) {
-                //no match possible for this barcode, or no match for previous barcodes
-                cbMatchGood=false;
-                cbMatch=-11;
+            if ( !cb.extractBarcode(bSeq, bQual, adapterStart, cbSeq1, cbQual1) 
+                 || cbSeq1.size() < cb.minLen || cbSeq1.size() >= cb.wl.size() || cb.wl[cbSeq1.size()].size()==0 ) {
+                //barcode cannot be extracted
+                if (cbMatchGood) {
+                    cbMatch=-11;
+                    cbMatchGood=false;
+                };
             };
             cbSeq  += cbSeq1 + "_";
             cbQual += cbQual1 + "_";

source/ReadAlign_waspMap.cpp

@@ -1 +1 @@
-#define STAR_VERSION "2.7.2x_0708_soloBamTags"
+#define STAR_VERSION "2.7.2x_0719_soloComplexBarcodes"

source/STAR.cpp

@@ -5,7 +5,7 @@
 #include "ErrorWarning.h"
 #include "bam_cat.h"
 
-void bamSortByCoordinate (Parameters &P, ReadAlignChunk **RAchunk, Genome &mainGenome, SoloFeature &soloFeat) {
+void bamSortByCoordinate (Parameters &P, ReadAlignChunk **RAchunk, Genome &mainGenome, Solo &solo) {
     if (P.outBAMcoord) {//sort BAM if needed
         *P.inOut->logStdOut << timeMonthDayTime() << " ..... started sorting BAM\n" <<flush;
         P.inOut->logMain << timeMonthDayTime() << " ..... started sorting BAM\n" <<flush;
@@ -54,7 +54,7 @@ void bamSortByCoordinate (Parameters &P, ReadAlignChunk **RAchunk, Genome &mainG
                 if (binS==0) continue; //empty bin
 
                 if (ibin == nBins-1) {//last bin for unmapped reads
-                    BAMbinSortUnmapped(ibin,P.runThreadN,P.outBAMsortTmpDir, P, mainGenome, soloFeat);
+                    BAMbinSortUnmapped(ibin,P.runThreadN,P.outBAMsortTmpDir, P, mainGenome, solo);
                 } else {
                     uint newMem=binS+binN*24;
                     bool boolWait=true;
@@ -66,7 +66,7 @@ void bamSortByCoordinate (Parameters &P, ReadAlignChunk **RAchunk, Genome &mainG
                         };
                         sleep(0.1);
                     };
-                    BAMbinSortByCoordinate(ibin,binN,binS,P.runThreadN,P.outBAMsortTmpDir, P, mainGenome, soloFeat);
+                    BAMbinSortByCoordinate(ibin,binN,binS,P.runThreadN,P.outBAMsortTmpDir, P, mainGenome, solo);
                     #pragma omp critical
                     totalMem-=newMem;//"release" RAM
                 };

source/SoloFeature_addBAMtags.cpp

@@ -4,8 +4,8 @@
 #include "Parameters.h"
 #include "ReadAlignChunk.h"
 #include "Genome.h"
-#include "SoloFeature.h"
+#include "Solo.h"
 
-void bamSortByCoordinate(Parameters &P, ReadAlignChunk **RAchunk, Genome &mainGenome, SoloFeature &soloFeat) ;
+void bamSortByCoordinate(Parameters &P, ReadAlignChunk **RAchunk, Genome &mainGenome, Solo &solo) ;
 
 #endif