CHANGES.md

+STAR 2.6.0c 2018/05/10
+======================
+
+* Fixed bugs in merging mates (--peOverlap*) and WASP filtering algorithms. Please see CHANGES and RELEASEnotes from 2.6.0a.
+
+
 STAR 2.6.0b 2018/05/02
 ======================
 
-* Fixed bugs introduced in 2.6.0a. Please see CHANGES and RELEASEnotes from 2.5.0a.
+* Fixed bugs introduced in 2.6.0a. Please see CHANGES and RELEASEnotes from 2.6.0a.
 
 
 STAR 2.6.0a 2018/04/23

debian/changelog

+rna-star (2.6.0c+dfsg-1) unstable; urgency=medium
+
+  * New upstream release.
+
+ -- Sascha Steinbiss <satta@debian.org>  Wed, 16 May 2018 08:39:24 +0200
+
 rna-star (2.6.0b+dfsg-1) unstable; urgency=medium
 
   * New upstream release.

source/Parameters.cpp

@@ -1039,75 +1039,76 @@ void Parameters::inputParameters (int argInN, char* argIn[]) {//input parameters
     pCh.out.junctions=false;
     pCh.out.samOld=false;
     pCh.out.bamHardClip=true;//default
-    for (const auto& type1 : pCh.out.type) {
-        if (type1=="WithinBAM") {
-            pCh.out.bam=true;
-        } else if (pCh.out.type.at(0)=="SeparateSAMold") {
-            pCh.out.samOld=true;
-        } else if (pCh.out.type.at(0)=="Junctions") {
-            pCh.out.junctions=true;
-        } else if (type1=="HardClip") {
-            pCh.out.bamHardClip=true;
-        } else if (type1=="SoftClip") {
-            pCh.out.bamHardClip=false;     
-        } else {
-            ostringstream errOut;
-            errOut <<"EXITING because of FATAL INPUT ERROR: unknown/unimplemented value for --chimOutType: "<<type1 <<"\n";
-            errOut <<"SOLUTION: re-run STAR with --chimOutType Junctions , SeparateSAMold  , WithinBAM , HardClip \n";
-            exitWithError(errOut.str(), std::cerr, inOut->logMain, EXIT_CODE_PARAMETER, *this);
+    if (pCh.segmentMin>0) {//only if chimeric detection is activated
+        for (const auto& type1 : pCh.out.type) {
+            if (type1=="WithinBAM") {
+                pCh.out.bam=true;
+            } else if (type1=="SeparateSAMold") {
+                pCh.out.samOld=true;
+            } else if (type1=="Junctions") {
+                pCh.out.junctions=true;
+            } else if (type1=="HardClip") {
+                pCh.out.bamHardClip=true;
+            } else if (type1=="SoftClip") {
+                pCh.out.bamHardClip=false;     
+            } else {
+                ostringstream errOut;
+                errOut <<"EXITING because of FATAL INPUT ERROR: unknown/unimplemented value for --chimOutType: "<<type1 <<"\n";
+                errOut <<"SOLUTION: re-run STAR with --chimOutType Junctions , SeparateSAMold  , WithinBAM , HardClip \n";
+                exitWithError(errOut.str(), std::cerr, inOut->logMain, EXIT_CODE_PARAMETER, *this);
+            };
         };
-    };
 
-    if (pCh.out.bam && !outBAMunsorted && !outBAMcoord) {
-            ostringstream errOut;
-            errOut <<"EXITING because of fatal PARAMETERS error: --chimOutType WithinBAM requires BAM output\n";
-            errOut <<"SOLUTION: re-run with --outSAMtype BAM Unsorted/SortedByCoordinate\n";
-            exitWithError(errOut.str(), std::cerr, inOut->logMain, EXIT_CODE_PARAMETER, *this);
-    };
-    
-    if (pCh.multimapNmax>0 && (pCh.out.bam || pCh.out.samOld)) {       
-            ostringstream errOut;
-            errOut <<"EXITING because of fatal PARAMETERS error: --chimMultimapNmax > 0 (new chimeric detection) presently only works with --chimOutType Junctions\n";
-            errOut <<"SOLUTION: re-run with --chimOutType Junctions\n";
-            exitWithError(errOut.str(), std::cerr, inOut->logMain, EXIT_CODE_PARAMETER, *this);
-    };
-    
-    if (peOverlap.NbasesMin > 0) {
-        if (pCh.multimapNmax == 0 && (pCh.out.junctions || pCh.out.samOld)) {
+        if (pCh.out.bam && !outBAMunsorted && !outBAMcoord) {
                 ostringstream errOut;
-                errOut <<"EXITING because of fatal PARAMETERS error: --chimMultimapNmax 0 (default old chimeric detection) and --peOverlapNbasesMin > 0 (merging ovelrapping mates) presently only works with --chimOutType WithinBAM\n";
-                errOut <<"SOLUTION: re-run with --chimOutType WithinBAM\n";
+                errOut <<"EXITING because of fatal PARAMETERS error: --chimOutType WithinBAM requires BAM output\n";
+                errOut <<"SOLUTION: re-run with --outSAMtype BAM Unsorted/SortedByCoordinate\n";
                 exitWithError(errOut.str(), std::cerr, inOut->logMain, EXIT_CODE_PARAMETER, *this);
         };
-    };
-        
-    
-    
-    if (pCh.out.bam && !outSAMattrPresent.NM) {
-       outSAMattrOrder.push_back(ATTR_NM);
-       inOut->logMain << "WARNING --chimOutType=WithinBAM, therefore STAR will output NM attribute" <<endl;
-    };
-   
-    pCh.filter.genomicN=false;
-    for (uint ii=0; ii<pCh.filter.stringIn.size(); ii++)
-    {
-        if (pCh.filter.stringIn.at(ii)=="banGenomicN")
-        {
-            pCh.filter.genomicN=true;
-        }
-        else if (pCh.filter.stringIn.at(ii)=="None")
-        {//nothing to do
-        }
-        else
+
+        if (pCh.multimapNmax>0 && (pCh.out.bam || pCh.out.samOld)) {       
+                ostringstream errOut;
+                errOut <<"EXITING because of fatal PARAMETERS error: --chimMultimapNmax > 0 (new chimeric detection) presently only works with --chimOutType Junctions\n";
+                errOut <<"SOLUTION: re-run with --chimOutType Junctions\n";
+                exitWithError(errOut.str(), std::cerr, inOut->logMain, EXIT_CODE_PARAMETER, *this);
+        };
+
+        if (peOverlap.NbasesMin > 0) {
+            if (pCh.multimapNmax == 0 && (pCh.out.junctions || pCh.out.samOld)) {
+                    ostringstream errOut;
+                    errOut <<"EXITING because of fatal PARAMETERS error: --chimMultimapNmax 0 (default old chimeric detection) and --peOverlapNbasesMin > 0 (merging ovelrapping mates) presently only works with --chimOutType WithinBAM\n";
+                    errOut <<"SOLUTION: re-run with --chimOutType WithinBAM\n";
+                    exitWithError(errOut.str(), std::cerr, inOut->logMain, EXIT_CODE_PARAMETER, *this);
+            };
+        };
+
+
+
+        if (pCh.out.bam && !outSAMattrPresent.NM) {
+           outSAMattrOrder.push_back(ATTR_NM);
+           inOut->logMain << "WARNING --chimOutType=WithinBAM, therefore STAR will output NM attribute" <<endl;
+        };
+
+        pCh.filter.genomicN=false;
+        for (uint ii=0; ii<pCh.filter.stringIn.size(); ii++)
         {
-            ostringstream errOut;
-            errOut << "EXITING because of fatal PARAMETERS error: unrecognized value of --chimFilter="<<pCh.filter.stringIn.at(ii)<<"\n";
-            errOut << "SOLUTION: use allowed values: banGenomicN || None";
-            exitWithError(errOut.str(),std::cerr, inOut->logMain, EXIT_CODE_PARAMETER, *this);
+            if (pCh.filter.stringIn.at(ii)=="banGenomicN")
+            {
+                pCh.filter.genomicN=true;
+            }
+            else if (pCh.filter.stringIn.at(ii)=="None")
+            {//nothing to do
+            }
+            else
+            {
+                ostringstream errOut;
+                errOut << "EXITING because of fatal PARAMETERS error: unrecognized value of --chimFilter="<<pCh.filter.stringIn.at(ii)<<"\n";
+                errOut << "SOLUTION: use allowed values: banGenomicN || None";
+                exitWithError(errOut.str(),std::cerr, inOut->logMain, EXIT_CODE_PARAMETER, *this);
+            };
         };
     };
 
-
     alignEndsType.ext[0][0]=false;
     alignEndsType.ext[0][1]=false;
     alignEndsType.ext[1][0]=false;

source/ReadAlign_peOverlapMergeMap.cpp

@@ -3,10 +3,9 @@
 
 void ReadAlign::peOverlapMergeMap() {
     
-
     if (!P.peOverlap.yes || P.readNmates!=2 ) {//no peOverlap
-        peOv.yes=false;
         return;
+        peOv.yes=false;
     };
 
     //debug
@@ -17,10 +16,10 @@ void ReadAlign::peOverlapMergeMap() {
     peMergeRA->copyRead(*this);
     peMergeRA->peMergeMates();
     peOv=peMergeRA->peOv;
+    peOv.yes=false;
 
     if (peOv.nOv==0) {//check if mates can be merged, if not - return
         //cout <<"\n-1\n";
-        peOv.yes=false;
         return;
     };
 
@@ -36,7 +35,6 @@ void ReadAlign::peOverlapMergeMap() {
         //    cout <<P.genomeNumToNT[peMergeRA->Read1[0][ii]];
         //};
         //cout << "\n";
-
         return;
     };
     
@@ -71,10 +69,8 @@ void ReadAlign::peOverlapMergeMap() {
 
     //P.alignSplicedMateMapLminOverLmate=P_alignSplicedMateMapLminOverLmate;
     
-    if (peScore<=trBest->maxScore || chimRecord) {
+    if (peScore<=trBest->maxScore || chimRecord) {//otherwise peOv.yes=false
         peOv.yes=true;
-    } else {
-        peOv.yes=false;
     };
     
     return;
@@ -82,8 +78,8 @@ void ReadAlign::peOverlapMergeMap() {
 
 void ReadAlign::peMergeMates() {
     
-    uint s1=localSearch(Read1[0],readLength[0],Read1[0]+readLength[0]+1,readLength[1],P.peOverlap.MMp);
-    uint s0=localSearch(Read1[0]+readLength[0]+1,readLength[1],Read1[0],readLength[0],P.peOverlap.MMp);
+    uint s1=localSearchNisMM(Read1[0],readLength[0],Read1[0]+readLength[0]+1,readLength[1],P.peOverlap.MMp);
+    uint s0=localSearchNisMM(Read1[0]+readLength[0]+1,readLength[1],Read1[0],readLength[0],P.peOverlap.MMp);
 
     uint o1=min(readLength[1],readLength[0]-s1);
     uint o0=min(readLength[0],readLength[1]-s0);
@@ -203,12 +199,12 @@ void Transcript::peOverlapSEtoPE(uint* mateStart, Transcript &t) {//convert alig
 
             exons[nExons][EX_iFrag]=(imate==0 ? t.Str : 1-t.Str);
             exons[nExons][EX_sjA]=t.exons[iex][EX_sjA];
-            if (nExons>0 && iex>0) {
-                canonSJ[nExons-1]=t.canonSJ[iex-1];
-                sjAnnot[nExons-1]=t.sjAnnot[iex-1];
-                sjStr[nExons-1]=t.sjStr[iex-1];
-                shiftSJ[nExons-1][0]=t.shiftSJ[iex-1][0];
-                shiftSJ[nExons-1][1]=t.shiftSJ[iex-1][1];
+            if (iex<t.nExons-1) {
+                canonSJ[nExons]=t.canonSJ[iex];
+                sjAnnot[nExons]=t.sjAnnot[iex];
+                sjStr[nExons]=t.sjStr[iex];
+                shiftSJ[nExons][0]=t.shiftSJ[iex][0];
+                shiftSJ[nExons][1]=t.shiftSJ[iex][1];
             };
             //record these exons for mate2
             if (t.exons[iex][EX_R]>=mSta[imate]) {//exon left is inside the mate

source/ReadAlign_waspMap.cpp

@@ -61,7 +61,9 @@ void ReadAlign::waspMap() {
             };
 
             waspRA->mapOneRead();
-
+            waspRA->multMapSelect();
+            waspRA->mappedFilter();
+                    
             if (waspRA->unmapType!=-1) {
                 waspType=4;
                 return;

source/STAR.cpp

@@ -48,6 +48,18 @@ void usage() {
 
 
 int main(int argInN, char* argIn[]) {
+    
+    //debug
+//     uint nn=4000000000;
+//     uint* dummy=new uint[nn];
+//     std::cout <<"Allocated"<<endl;
+//     for (uint ii=0;ii<nn;ii++) {
+//         dummy[ii]=14829735431805718528LLU;;
+//     };
+//     std::cout <<"Filled "<<dummy[0]<<endl;
+//     delete[] dummy;     
+//     std::cout <<"Deleted"<<endl;
+    
     // If no argument is given, or the first argument is either '-h' or '--help', run usage()
     if((argInN == 1) || (argInN == 2 && (strcmp("-h",argIn[1]) == 0 || strcmp ("--help",argIn[1]) == 0 ))) {
         usage();

source/SequenceFuns.cpp

@@ -201,6 +201,32 @@ uint localSearch(const char *x, uint nx, const char *y, uint ny, double pMM){
     return ixBest;
 };
 
+uint localSearchNisMM(const char *x, uint nx, const char *y, uint ny, double pMM){
+    //find the best alignment of two short sequences x and y
+    //pMM is the maximum percentage of mismatches
+    //Ns in x OR y are considered mismatches
+    uint nMatch=0, nMM=0, nMatchBest=0, nMMbest=0, ixBest=nx;
+    for (uint ix=0;ix<nx;ix++) {
+        nMatch=0; nMM=0;
+        for (uint iy=0;iy<min(ny,nx-ix);iy++) {
+            if (x[ix+iy]==y[iy] && y[iy]<4) {
+                nMatch++;
+            } else {
+                nMM++;
+            };
+        };
+
+        if ( ( nMatch>nMatchBest || (nMatch==nMatchBest && nMM<nMMbest) ) && double(nMM)/double(nMatch)<=pMM) {
+            ixBest=ix;
+            nMatchBest=nMatch;
+            nMMbest=nMM;
+        };
+    };
+    return ixBest;
+};
+
+
+
 uint qualitySplit(char* r, char* q, uint L, char Qsplit, uint maxNsplit, uint  minLsplit, uint** splitR) {
     //splits the read r[L] by quality scores q[L], outputs in splitR - split coordinate/length - per base
     //returns number of good split regions

source/SequenceFuns.h

@@ -21,6 +21,8 @@ char convertNt01234(const char R0);//transform sequence  from ACGT into 0-1-2-3-
 
 uint chrFind(uint, uint, uint*); // find chromosome from global locus
 uint localSearch(const char*, uint, const char*, uint, double); //local search to clip adapter
+uint localSearchNisMM(const char *x, uint nx, const char *y, uint ny, double pMM);
+
 uint qualitySplit(char*, char*, uint, char, uint, uint, uint**);
 
 #endif

source/VERSION

-#define STAR_VERSION "STAR_2.6.0b"
+#define STAR_VERSION "STAR_2.6.0c"