+ fixed issue with negative strand in lookup table
[qpalma.git] / qpalma / Lookup.py
index d9ea7be..baa6419 100644 (file)
@@ -52,24 +52,27 @@ class LookupTable:
       assert chromo in self.chromo_list
       assert start <= stop
 
-      if strand == '-':
-         total_size = self.seqInfo.getFragmentSize(chromo)
-         startP   = total_size - stop
-         stopP    = total_size - start
-         _stop     = stopP
-         _start    = startP
-
       chromo_idx = self.chromo_map[chromo]
 
+
       if strand == '+':
          currentSequence         = self.genomicSequences[chromo_idx][start:stop]
          currentAcceptorScores   = self.acceptorScoresPos[chromo_idx][start:stop]
          currentDonorScores      = self.donorScoresPos[chromo_idx][start:stop]
       elif strand == '-':
-         currentSequence         = self.genomicSequences[chromo_idx][_start:_stop]
+         currentSequence         = self.genomicSequences[chromo_idx][start:stop]
          currentSequence = reverse_complement(currentSequence)
-         currentAcceptorScores   = self.acceptorScoresNeg[chromo_idx][start:stop]
-         currentDonorScores      = self.donorScoresNeg[chromo_idx][start:stop]
+
+         total_size = self.seqInfo.getFragmentSize(chromo)
+         _start = total_size - stop
+         _stop = total_size - start
+         currentAcceptorScores   = self.acceptorScoresNeg[chromo_idx][_start:_stop]
+         currentDonorScores      = self.donorScoresNeg[chromo_idx][_start:_stop]
+
+         #gt_tuple_pos = [p for p,e in enumerate(currentSequence) if (p>0 and p<len(currentSequence)-1 and e=='g' ) and (currentSequence[p+1]=='t' or currentSequence[p+1]=='c')]
+         #ag_tuple_pos = [p for p,e in enumerate(currentSequence) if p>1 and currentSequence[p-1]=='a' and currentSequence[p]=='g']
+         #acc_pos = [p for p,e in enumerate(currentAcceptorScores) if e != -inf and p > 1]
+         #don_pos = [p for p,e in enumerate(currentDonorScores) if e != -inf and p > 0]
 
       return currentSequence, currentAcceptorScores, currentDonorScores
 
@@ -86,9 +89,7 @@ class LookupTable:
 
       print 'lt total_size %d' % self.seqInfo.getFragmentSize(chromo)
 
-      strand = '+'
-      genomicSeq,currentAcc,currentDon = self.seqInfo.get_seq_and_scores(chromo,strand,genomicSeq_start,genomicSeq_stop)
-      #print len(currentAcc),len(currentDon)
+      genomicSeq,currentAcc,currentDon = self.seqInfo.get_seq_and_scores(chromo,'+',genomicSeq_start,genomicSeq_stop)
       genomicSeq = genomicSeq.lower()
 
       chromo_idx = self.chromo_map[chromo]
@@ -97,8 +98,14 @@ class LookupTable:
       self.acceptorScoresPos[chromo_idx]   = currentAcc
       self.donorScoresPos[chromo_idx]      = currentDon
 
-      strand = '-'
-      genomicSeq,currentAcc,currentDon = self.seqInfo.get_seq_and_scores(chromo,strand,genomicSeq_start,genomicSeq_stop)
+      dummy,currentAcc,currentDon = self.seqInfo.get_seq_and_scores(chromo,'-',genomicSeq_start,genomicSeq_stop)
+
+      #genomicSeq = reverse_complement(self.genomicSequences[chromo_idx])
+      #gt_tuple_pos = [p for p,e in enumerate(genomicSeq) if (p>0 and p<len(genomicSeq)-1 and e=='g' ) and (genomicSeq[p+1]=='t' or genomicSeq[p+1]=='c')]
+      #ag_tuple_pos = [p for p,e in enumerate(genomicSeq) if p>1 and genomicSeq[p-1]=='a' and genomicSeq[p]=='g']
+      #acc_pos = [p for p,e in enumerate(currentAcc) if e != -inf and p > 1]
+      #don_pos = [p for p,e in enumerate(currentDon) if e != -inf and p > 0]
+      #pdb.set_trace()
 
       self.acceptorScoresNeg[chromo_idx]   = currentAcc
       self.donorScoresNeg[chromo_idx]      = currentDon