+ fixed minor bugs
authorFabio <fabio@congo.fml.local>
Tue, 14 Oct 2008 16:38:29 +0000 (18:38 +0200)
committerFabio <fabio@congo.fml.local>
Tue, 14 Oct 2008 16:38:29 +0000 (18:38 +0200)
qpalma/Lookup.py
tests/test_qpalma.py

index 2c7fd35..d9ea7be 100644 (file)
@@ -48,30 +48,28 @@ class LookupTable:
          self.prefetch_seq_and_scores(chrId)
 
  
-   def get_seq_and_scores(self,chromo,strand,_start,_stop):
+   def get_seq_and_scores(self,chromo,strand,start,stop):
       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
-
-      assert start <= stop
+         startP   = total_size - stop
+         stopP    = total_size - start
+         _stop     = stopP
+         _start    = startP
 
       chromo_idx = self.chromo_map[chromo]
 
-      currentSequence         = self.genomicSequences[chromo_idx][start:stop]
       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 == '-':
-         currentAcceptorScores   = self.acceptorScoresNeg[chromo_idx][_start:_stop]
-         currentDonorScores      = self.donorScoresNeg[chromo_idx][_start:_stop]
-
-      if strand == '-':
+         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]
 
       return currentSequence, currentAcceptorScores, currentDonorScores
 
@@ -83,7 +81,6 @@ class LookupTable:
       """
       assert chromo in self.chromo_list
 
-      
       genomicSeq_start  = 0
       genomicSeq_stop   = self.seqInfo.getFragmentSize(chromo)
 
@@ -103,16 +100,5 @@ class LookupTable:
       strand = '-'
       genomicSeq,currentAcc,currentDon = self.seqInfo.get_seq_and_scores(chromo,strand,genomicSeq_start,genomicSeq_stop)
 
-      genomicSeq = genomicSeq.lower()
-
-      newCurrentAcc = [-inf] 
-      newCurrentAcc.extend(currentAcc[:-1])
-      currentAcc = newCurrentAcc
-
-      newCurrentDon = [-inf] 
-      newCurrentDon.extend(currentDon[:-1])
-      #pdb.set_trace()
-      currentDon = newCurrentDon
-
       self.acceptorScoresNeg[chromo_idx]   = currentAcc
       self.donorScoresNeg[chromo_idx]      = currentDon
index e585a75..1bb101b 100644 (file)
@@ -318,12 +318,20 @@ def check_reverse_strand_calculation(id,b,e,seqInfo):
 
 def check_example(chromo,strand,b,e,seqInfo,lt1):
    dna,acc,don = seqInfo.get_seq_and_scores(1,strand,b,e)
-   #_dna,_acc,_don = seqInfo.get_seq_and_scores(1,strand,b,e)
-   _dna,_acc,_don= lt1.get_seq_and_scores(1,'-',b,e)
+
+   if lt1 != None:
+      _dna,_acc,_don= lt1.get_seq_and_scores(1,strand,b,e)
+   else:
+      _dna,_acc,_don = seqInfo.get_seq_and_scores(1,strand,b,e)
 
    print 'Current interval: (%d,%d), current strand: %s'%(b,e,strand)
    print 'Results for dna,acc,don:'
    print '%s %s %s'%(str(dna==_dna),str(acc==_acc),str(don==_don))
+   print 'size is %d' % len(dna)
+   if dna != _dna:
+      print dna[:20]
+      print _dna[:20]
+
    print [p for p,e in enumerate(acc) if e != -inf][:10]
    print [p for p,e in enumerate(_acc) if e != -inf][:10]
    print [p for p,e in enumerate(don) if e != -inf][:10]
@@ -340,7 +348,7 @@ def simple_check(settings,seqInfo,lt1):
    total_size = seqInfo.getFragmentSize(chromo)
 
    for strand in ['+','-']:
-      for (b,e) in [(206,874),(545,874),(999,1234),(1000,total_size-1000),(3,total_size-3)]:
+      for (b,e) in [(206,874),(545,874),(999,1234),(1000,total_size-1000),(3,total_size-3),(0,total_size)]:
          check_example(chromo,strand,b,e,seqInfo,lt1)
 
    for (b,e) in intervals:
@@ -363,7 +371,7 @@ def checks():
    accessWrapper = DataAccessWrapper(settings)
    seqInfo = SeqSpliceInfo(accessWrapper,settings['allowed_fragments'])
 
-   lt1 = None
+   #lt1 = None
    lt1 = LookupTable(settings)
 
    print 'Checking with toy data...'
@@ -385,10 +393,10 @@ def checks():
    seqInfo = SeqSpliceInfo(accessWrapper,settings['allowed_fragments'])
 
    lt1 = None
-   lt1 = LookupTable(settings)
+   #lt1 = LookupTable(settings)
 
-   print 'Checking with real data...'
-   simple_check(settings,seqInfo,lt1)
+   #print 'Checking with real data...'
+   #simple_check(settings,seqInfo,lt1)
 
 
 def run():