+ added checks to see whether reads lie in the 250nuc are from the start/stop of...
authorfabio <fabio@e1793c9e-67f9-0310-80fc-b846ff1f7b36>
Fri, 20 Jun 2008 13:32:31 +0000 (13:32 +0000)
committerfabio <fabio@e1793c9e-67f9-0310-80fc-b846ff1f7b36>
Fri, 20 Jun 2008 13:32:31 +0000 (13:32 +0000)
sequence as there are no splice scores present

git-svn-id: http://svn.tuebingen.mpg.de/ag-raetsch/projects/QPalma@9695 e1793c9e-67f9-0310-80fc-b846ff1f7b36

qpalma/sequence_utils.py

index 46435e1..f9043ec 100644 (file)
@@ -1,11 +1,16 @@
 #!/usr/bin/env python 
 # -*- coding: utf-8 -*- 
 
+#
+# This module holds all functions for queries on the dna flat files and the
+# splice score files.
+# Furthermore it contains some utility functions such as reverse complement and
+# unbracketing est strings
+#
+
 import os
 import pdb
-import random
 import re
-import sys
 import subprocess
 
 from numpy.matlib import inf
@@ -13,6 +18,7 @@ from numpy.matlib import inf
 from Genefinding import *
 from genome_utils import load_genomic
 
+
 def get_flatfile_size(filename):
    cmd =  'wc -c %s | cut -f1 -d \' \'' % filename
    obj = subprocess.Popen(cmd,shell=True,stdout=subprocess.PIPE,stderr=subprocess.PIPE)
@@ -104,10 +110,11 @@ def create_bracket_seq(dna_seq,read_seq):
    read: aac
    is written in bracket notation: aa[ac]
    """
-   assert len(dna_seq) == len(read_seq)
+   assert len(dna_seq) == len(read_seq),pdb.set_trace()
    return "".join(map(lambda x,y: ['[%s%s]'%(x,y),x][x==y],dna_seq,read_seq))
 
 
+
 def getSpliceScores(chr,strand,intervalBegin,intervalEnd,total_size=0):
    """
    Now we want to use interval_query to get the predicted splice scores trained
@@ -135,12 +142,17 @@ def getSpliceScores(chr,strand,intervalBegin,intervalEnd,total_size=0):
    return acc, don
 
 
-def get_seq_and_scores(chr,strand,genomicSeq_start,genomicSeq_stop,dna_flat_files):
+def get_seq_and_scores(chr,strand,genomicSeq_start,genomicSeq_stop,dna_flat_files,only_seq=False):
    """
    This function expects an interval, chromosome and strand information and
    returns then the genomic sequence of this interval and the associated scores.
    """
 
+   # because splice site scores are predicted using a window of the sequence the
+   # first and the last NO_SCORE_WINDOW_SIZE nucleotides of a genomic sequence
+   # do not have any score predictions
+   NO_SCORE_WINDOW_SIZE = 250
+
    assert genomicSeq_start < genomicSeq_stop
 
    chrom         = 'chr%d' % chr
@@ -154,47 +166,48 @@ def get_seq_and_scores(chr,strand,genomicSeq_start,genomicSeq_stop,dna_flat_file
    no_base = re.compile('[^acgt]')
    genomicSeq = no_base.sub('n',genomicSeq)
 
-   if strand == '+':
-      intervalBegin  = genomicSeq_start-100
-      intervalEnd    = genomicSeq_stop+100
+   if only_seq:
+      return genomicSeq
 
-      currentAcc, currentDon = getSpliceScores(chr,strand,intervalBegin,intervalEnd)
+   fn = 'chr%d.dna.flat' % chr
+   filename = os.path.join(dna_flat_files,fn)
+   total_size = get_flatfile_size(filename)
 
-      currentAcc = currentAcc[100:-98]
-      currentAcc = currentAcc[1:]
-      currentDon = currentDon[100:-100]
+   intervalBegin = genomicSeq_start-100
+   intervalEnd    = genomicSeq_stop+100
 
-      length = len(genomicSeq)
-      currentAcc = currentAcc[:length]
+   # shorten the intervalEnd to the maximal size of the genomic sequence if it
+   # would exceed this
+   if intervalEnd > total_size:
+      intervalEnd = total_size-1
 
+   # we have assertions below checking whether every gt or ag has a splice score
+   # not equal to -inf if the gt or ag are within the first NO_SCORE_WINDOW_SIZE
+   # nucleotide then they have no scores at all
+   check_assertions = True
+   if intervalEnd > total_size - NO_SCORE_WINDOW_SIZE:
+      check_assertions = False
+      
+   if intervalBegin < NO_SCORE_WINDOW_SIZE:
+      check_assertions = False
+
+   currentAcc, currentDon = getSpliceScores(chr,strand,intervalBegin,intervalEnd,total_size)
+
+   currentAcc = currentAcc[100:-98]
+   currentAcc = currentAcc[1:]
+   currentDon = currentDon[100:-100]
+
+   length = len(genomicSeq)
+
+   if strand == '+':
+      currentAcc = currentAcc[:length]
       currentDon = currentDon+[-inf]*(length-len(currentDon))
 
       ag_tuple_pos = [p for p,e in enumerate(genomicSeq) if p>1 and genomicSeq[p-1]=='a' and genomicSeq[p]=='g' ]
       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')]
 
-      assert ag_tuple_pos == [p for p,e in enumerate(currentAcc) if e != -inf and p > 1], pdb.set_trace()
-      assert gt_tuple_pos == [p for p,e in enumerate(currentDon) if e != -inf and p > 0], pdb.set_trace()
-      assert len(currentAcc) == len(currentDon)
-
-      return genomicSeq, currentAcc, currentDon
-
    # build reverse complement if on negative strand
    if strand == '-':
-      fn = 'chr%d.dna.flat' % chr
-      filename = os.path.join(dna_flat_files,fn)
-      end = get_flatfile_size(filename)
-
-      intervalBegin = genomicSeq_start-100
-      intervalEnd    = genomicSeq_stop+100
-
-      total_size = end
-      currentAcc, currentDon = getSpliceScores(chr,strand,intervalBegin,intervalEnd,total_size)
-
-      currentAcc = currentAcc[100:-98]
-      currentAcc = currentAcc[1:]
-      currentDon = currentDon[100:-100]
-
-      length = len(genomicSeq)
       currentAcc = currentAcc[:length]
       currentAcc = currentAcc+[-inf]*(length-len(currentAcc))
 
@@ -203,11 +216,12 @@ def get_seq_and_scores(chr,strand,genomicSeq_start,genomicSeq_stop,dna_flat_file
       ag_tuple_pos = [p for p,e in enumerate(genomicSeq) if p>1 and p<len(genomicSeq) and genomicSeq[p-1]=='a' and genomicSeq[p]=='g']
       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')]
 
+   if check_assertions:
       assert ag_tuple_pos == [p for p,e in enumerate(currentAcc) if e != -inf and p > 1], pdb.set_trace()
       assert gt_tuple_pos == [p for p,e in enumerate(currentDon) if e != -inf and p > 0], pdb.set_trace()
       assert len(currentAcc) == len(currentDon), pdb.set_trace()
 
-      return genomicSeq, currentAcc, currentDon
+   return genomicSeq, currentAcc, currentDon
 
 
 def checks():
@@ -222,6 +236,5 @@ def checks():
          print p_seq == reverse_complement(n_seq)
 
          
-
 if __name__ == '__main__':
    checks()