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,12 +1039,13 @@ void Parameters::inputParameters (int argInN, char* argIn[]) {//input parameters
     pCh.out.junctions=false;
     pCh.out.samOld=false;
     pCh.out.bamHardClip=true;//default
+    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 (pCh.out.type.at(0)=="SeparateSAMold") {
+            } else if (type1=="SeparateSAMold") {
                 pCh.out.samOld=true;
-        } else if (pCh.out.type.at(0)=="Junctions") {
+            } else if (type1=="Junctions") {
                 pCh.out.junctions=true;
             } else if (type1=="HardClip") {
                 pCh.out.bamHardClip=true;
@@ -1106,7 +1107,7 @@ void Parameters::inputParameters (int argInN, char* argIn[]) {//input parameters
                 exitWithError(errOut.str(),std::cerr, inOut->logMain, EXIT_CODE_PARAMETER, *this);
             };
         };
-
+    };
 
     alignEndsType.ext[0][0]=false;
     alignEndsType.ext[0][1]=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,6 +61,8 @@ void ReadAlign::waspMap() {
             };
 
             waspRA->mapOneRead();
+            waspRA->multMapSelect();
+            waspRA->mappedFilter();
                     
             if (waspRA->unmapType!=-1) {
                 waspType=4;

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"