basic XS inference for ref-overlapping alns

gpertea · gpertea · commit ba97372eb5cf · 2018-09-04T12:24:35.000-04:00
diff --git a/gclib/GBam.cpp b/gclib/GBam.cpp
@@ -388,9 +388,8 @@ switch (cop) {
 
  char GBamRecord::spliceStrand() { // '+', '-' from the XS tag, or 0 if no XS tag
    char c=tag_char("XS");
-   if (c) return c;
-     else return '.';
-   }
+   return ((c=='+' || c=='-') ? c : '.');
+ }
 
  char* GBamRecord::sequence() { //user must free this after use
    char *s = (char*)bam1_seq(b);
diff --git a/rlink.h b/rlink.h
@@ -428,7 +428,8 @@ struct GReadAlnData {
 	//GPVec< GVec<RC_ExonOvl> > g_exonovls; //>5bp overlaps with guide exons, for each read "exon"
 	GReadAlnData(GBamRecord* bamrec=NULL, char nstrand=0, int num_hits=0,
 			int hit_idx=0, TAlnInfo* tif=NULL):brec(bamrec), strand(nstrand),
-					nh(num_hits), hi(hit_idx), juncs(true), tinfo(tif) { } //, g_exonovls(true)
+					nh(num_hits), hi(hit_idx), juncs(true), tinfo(tif) { }
+	void clear() { delete tinfo; }
 	~GReadAlnData() { delete tinfo; }
 };
 
@@ -526,18 +527,13 @@ struct BundleData {
 	  			rc_tdata, rc_edata, rc_idata);
 	  }
  }
- /* after reference annotation was loaded
- void rc_finalize_refs() {
-     if (rc_data==NULL) return;
-     //rc_data->setupCov();
-	}
-	Not needed here, we update the coverage span as each transcript is added
- */
+
  void keepGuide(GffObj* scaff, GPVec<RC_TData>* rc_tdata=NULL,
 		 GPVec<RC_Feature>* rc_edata=NULL, GPVec<RC_Feature>* rc_idata=NULL);
 
- //bool evalReadAln(GBamRecord& brec, char& strand, int nh); //, int hi);
  bool evalReadAln(GReadAlnData& alndata, char& strand);
+   //update read alndata with info about reference overlaps
+   //also updates the transcription strand if needed/possible
 
  void Clear() {
 	keepguides.Clear();
diff --git a/stringtie.cpp b/stringtie.cpp
@@ -406,7 +406,6 @@ if (tstackSize<DEF_TSTACK_SIZE) defStackSize=DEF_TSTACK_SIZE;
 	 bool chr_changed=false;
 	 int pos=0;
 	 const char* refseqName=NULL;
-	 //char strand=0;
 	 char xstrand=0;
 	 int nh=1;
 	 int hi=0;
@@ -441,7 +440,8 @@ if (tstackSize<DEF_TSTACK_SIZE) defStackSize=DEF_TSTACK_SIZE;
 		 }
 		 if (xstrand=='.' && brec->exons.Count()>1) {
 			 no_xs++;
-			 continue; //skip spliced alignments lacking XS tag (e.g. HISAT alignments)
+			 //continue; //skip spliced alignments lacking XS tag (e.g. HISAT alignments) ?
+			 //this is changed in StringTie 2 !
 		 }
 		 if (refseqName==NULL) GError("Error: cannot retrieve target seq name from BAM record!\n");
 		 pos=brec->start; //BAM is 0 based, but GBamRecord makes it 1-based
@@ -461,10 +461,12 @@ if (tstackSize<DEF_TSTACK_SIZE) defStackSize=DEF_TSTACK_SIZE;
 			 if (alncounts.Count()<=gseq_id) {
 				 alncounts.Resize(gseq_id+1, 0);
 			 }
-			 else if (alncounts[gseq_id]>0) GError(ERR_BAM_SORT);
+			 else if (alncounts[gseq_id]>0)
+				 GError(ERR_BAM_SORT);
 			 prev_pos=0;
 		 }
-		 if (pos<prev_pos) GError(ERR_BAM_SORT);
+		 if (pos<prev_pos)
+			 GError(ERR_BAM_SORT);
 		 prev_pos=pos;
 		 if (skipGseq) continue;
 		 alncounts[gseq_id]++;
@@ -657,19 +659,18 @@ if (tstackSize<DEF_TSTACK_SIZE) defStackSize=DEF_TSTACK_SIZE;
 				 }
 			 } while (cend_changed);
 		 }
-	 } //adjusted currentend and checked for overlapping reference transcripts
+	 } //adjusted currentend ; selected overlapping reference transcripts
 	 GReadAlnData alndata(brec, 0, nh, hi, tinfo);
      bool ovlpguide=bundle->evalReadAln(alndata, xstrand);
-     if(!eonly || ovlpguide) { // in eonly case consider read only if it overlaps guide
-    	 //check for overlaps with ref transcripts which may set xstrand
+	 //overlaps with ref transcripts may set xstrand if needed
+     if (!eonly || ovlpguide) {
+    	 // in -e mode: a read alignment is discarded if it doesn't overlap a reference transcript!
     	 if (xstrand=='+') alndata.strand=1;
     	 else if (xstrand=='-') alndata.strand=-1;
-    	 //GMessage("%s\t%c\t%d\thi=%d\n",brec->name(), xstrand, alndata.strand,hi);
-    	 //countFragment(*bundle, *brec, hi,nh); // we count this in build_graphs to only include mapped fragments that we consider correctly mapped
-    	 //fprintf(stderr,"fragno=%d fraglen=%lu\n",bundle->num_fragments,bundle->frag_len);if(bundle->num_fragments==100) exit(0);
-    	 //if (!ballgown || ref_overlap)
-    	 processRead(currentstart, currentend, *bundle, hashread, alndata);
-			  // *brec, strand, nh, hi);
+    	 //this is changed in StringTie 2 ! (improved transcription strand guessing)
+    	 bool XSmissing=(alndata.strand==0 && brec->exons.Count()>1);
+    	 if (!XSmissing)
+    	   processRead(currentstart, currentend, *bundle, hashread, alndata);
      }
  } //for each read alignment
 
