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
"""
assert chromo in self.chromo_list
-
genomicSeq_start = 0
genomicSeq_stop = self.seqInfo.getFragmentSize(chromo)
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
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]
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:
accessWrapper = DataAccessWrapper(settings)
seqInfo = SeqSpliceInfo(accessWrapper,settings['allowed_fragments'])
- lt1 = None
+ #lt1 = None
lt1 = LookupTable(settings)
print 'Checking with toy data...'
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():