+ fixed alignment index calculation for negative strand
[qpalma.git] / qpalma / sequence_utils.py
index f9043ec..806413f 100644 (file)
@@ -1,6 +1,14 @@
-#!/usr/bin/env python 
-# -*- coding: utf-8 -*- 
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# Written (W) 2008 Fabio De Bona
+# Copyright (C) 2008 Max-Planck-Society
 
+#
+# This file contains the main interface to the QPalma pipeline.
+#
 #
 # This module holds all functions for queries on the dna flat files and the
 # splice score files.
 #
 
 import os
+import os.path
 import pdb
 import re
 import subprocess
+import numpy
 
 from numpy.matlib import inf
 
-from Genefinding import *
-from genome_utils import load_genomic
+from Genefinding import doIntervalQuery,new_intp,intp_assign
+
+from QPalmaDP import createIntArrayFromList
+
+
+def load_genomic(chromosome, strand, start, stop, genome, one_based=True):
+   """
+   This function stems from Cheng Soon Ongs genome_utils package.
+   """
+
+   if (type(start)==numpy.ndarray) or (type(stop)==numpy.ndarray):
+      assert(len(start) == len(stop))
+      assert((start[1:]-stop[:-1]>0).all())
+      if strand == '+':
+          idx = xrange(len(start))
+      else:
+          idx = xrange(len(start)-1,-1,-1)
+
+      seq = ''.join([load_genomic(chromosome, strand, start[ix], stop[ix], genome)\
+                     for ix in idx])
+      return seq
+  
+      #fname = '/fml/ag-raetsch/share/databases/genomes/' + genome + '/' + chromosome + '.dna.flat'
+      fname = genome + '/' + chromosome + '.dna.flat'
+      f=file(fname)
+      if one_based:
+          f.seek(start-1)
+          str=f.read(stop-start+1)
+      else:
+          f.seek(start)
+          str=f.read(stop-start)
+  
+      if strand=='-':
+          return reverse_complement(str)
+      elif strand=='+':
+          return str
+      else:
+          print 'strand must be + or -'
+          raise KeyError
+
+
+extended_alphabet    = ['-','a','c','g','t','n','[',']']
+alphabet             = ['-','a','c','g','t','n']
 
 
-def get_flatfile_size(filename):
+class DataAccessionException:
+   pass
+
+
+def get_flat_file_size(filename):
    cmd =  'wc -c %s | cut -f1 -d \' \'' % filename
    obj = subprocess.Popen(cmd,shell=True,stdout=subprocess.PIPE,stderr=subprocess.PIPE)
    out,err = obj.communicate()
    
    if err != '':
       print 'Error occurred while trying to obtain file size'
-   return int(out)
+      raise DataAccessionException
+
+   # subtract one because of wc -c output
+   return int(out)-1
 
 
 def reverse_complement(seq):
@@ -43,7 +101,7 @@ def reverse_complement(seq):
    """
 
    bpos = seq.find('[')
-   rc = lambda x: {'a':'t','c':'g','g':'c','t':'a'}[x]
+   rc = lambda x: {'a':'t','c':'g','g':'c','t':'a','n':'n'}[x]
 
    # check first whether seq contains no brackets at all
    if bpos == -1:
@@ -51,7 +109,7 @@ def reverse_complement(seq):
       ret_val.reverse()
       ret_val = "".join(ret_val)
    else:
-      brc = lambda x: {'a':'t','c':'g','g':'c','t':'a','[':'[',']':']'}[x]
+      brc = lambda x: {'a':'t','c':'g','g':'c','t':'a','n':'n','[':'[',']':']'}[x]
 
       # first_part is the part of the seq up to the first occurrence of a
       # bracket
@@ -69,12 +127,12 @@ def reverse_complement(seq):
    return ret_val
 
 
-def unbracket_est(est):
+def unbracket_seq(seq):
    """
    This function takes a read in bracket notation and restores the read sequence from it.
    I.e.
 
-   est = cgt[ac][tg]aa
+   seq = cgt[ac][tg]aa
 
    leads to
 
@@ -84,21 +142,44 @@ def unbracket_est(est):
    entry is the base from the dna.
    """
 
-   new_est = ''
+   new_seq = ''
+   e = 0
+
+   while True:
+      if e >= len(seq):
+         break
+
+      if seq[e] == '[':
+         new_seq += seq[e+2]
+         e += 4
+      else:
+         new_seq += seq[e]
+         e += 1
+
+   return "".join(new_seq).lower()
+
+
+def reconstruct_dna_seq(seq):
+   """
+
+   """
+
+   new_seq = ''
    e = 0
 
    while True:
-      if e >= len(est):
+      if e >= len(seq):
          break
 
-      if est[e] == '[':
-         new_est += est[e+2]
+      if seq[e] == '[':
+         new_seq += seq[e+1]
          e += 4
       else:
-         new_est += est[e]
+         new_seq += seq[e]
          e += 1
 
-   return "".join(new_est).lower()
+   return "".join(new_seq).lower()
+
 
 
 def create_bracket_seq(dna_seq,read_seq):
@@ -115,126 +196,366 @@ def create_bracket_seq(dna_seq,read_seq):
 
 
 
-def getSpliceScores(chr,strand,intervalBegin,intervalEnd,total_size=0):
+def my_load_genomic(fname, strand, start, stop, one_based=True):
+
+   if (type(start)==numpy.ndarray) or (type(stop)==numpy.ndarray):
+      assert(len(start) == len(stop))
+      assert((start[1:]-stop[:-1]>0).all())
+      if strand == '+':
+         idx = xrange(len(start))
+      else:
+         idx = xrange(len(start)-1,-1,-1)
+
+      seq = ''.join([load_genomic(chromosome, strand, start[ix], stop[ix], genome)\
+      for ix in idx])
+      return seq
+
+   try:
+      f=file(fname)
+      if one_based:
+         f.seek(start-1)
+         str=f.read(stop-start+1)
+      else:
+         f.seek(start)
+         str=f.read(stop-start)
+
+      if strand=='-':
+         return reverse_complement(str)
+      elif strand=='+':
+         return str
+      else:
+         print 'strand must be + or -'
+         raise KeyError
+
+   except:
+      print fname,strand,start,stop
+
+
+class DataAccessWrapper:
    """
-   Now we want to use interval_query to get the predicted splice scores trained
-   on the TAIR sequence and annotation.
+   The purpose of this class is to wrap the genomic/splice site score data
+   access.
    """
 
-   size = intervalEnd-intervalBegin
-   assert size > 1, 'Error (getSpliceScores): interval size is less than 2!'
+   def __init__(self,settings):
+      self.genomic_dir        = settings['genome_dir'] 
+      self.acc_score_dir      = settings['acceptor_scores_loc'] 
+      self.don_score_dir      = settings['donor_scores_loc'] 
+      self.genomic_fmt        = settings['genome_file_fmt'] 
+      self.sscore_fmt         = settings['splice_score_file_fmt'] 
 
-   acc = size*[0.0]
-   don = size*[0.0]
 
-   interval_matrix = createIntArrayFromList([intervalBegin,intervalEnd])
-   pos_size = new_intp()
-   intp_assign(pos_size,1)
+   def get_genomic_fragment_fn(self,id):
+      return os.path.join(self.genomic_dir,self.genomic_fmt%id)
 
-   # fetch acceptor scores
-   sscore_filename = '/fml/ag-raetsch/home/fabio/tmp/interval_query_files/acc/contig_%d%s'
-   acc = doIntervalQuery(chr,strand,intervalBegin,intervalEnd,sscore_filename,total_size)
 
-   # fetch donor scores
-   sscore_filename = '/fml/ag-raetsch/home/fabio/tmp/interval_query_files/don/contig_%d%s'
-   don = doIntervalQuery(chr,strand,intervalBegin,intervalEnd,sscore_filename,total_size)
+   def get_splice_site_scores_fn(self,id,strand):
+      acc_fn = os.path.join(self.acc_score_dir,self.sscore_fmt%(id,strand))
+      don_fn = os.path.join(self.don_score_dir,self.sscore_fmt%(id,strand))
+      return acc_fn,don_fn
 
-   return acc, don
 
 
-def get_seq_and_scores(chr,strand,genomicSeq_start,genomicSeq_stop,dna_flat_files,only_seq=False):
+class SeqSpliceInfo():
    """
-   This function expects an interval, chromosome and strand information and
-   returns then the genomic sequence of this interval and the associated scores.
+   
    """
+   
+   def __init__(self,accessWrapper,fragment_list):
+      self.accessWrapper = accessWrapper
+      
+      self.chromo_sizes = {}
+
+      # take care that it also works say if we get a chromo_list [1,3,5]
+      # those are mapped then to 0,1,2
+      self.fragment_list = fragment_list
+      self.chromo_map = {}
+      for idx,elem in enumerate(fragment_list):
+         self.chromo_map[elem] = idx
+   
+      print "Fetching fragment sizes..."
+      for id in fragment_list:
+         filename = self.accessWrapper.get_genomic_fragment_fn(id)
+         total_size = get_flat_file_size(filename)
+         self.chromo_sizes[self.chromo_map[id]] = total_size
 
-   # 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
+      print "done"
 
-   assert genomicSeq_start < genomicSeq_stop
 
-   chrom         = 'chr%d' % chr
-   genomicSeq = load_genomic(chrom,strand,genomicSeq_start-1,genomicSeq_stop,dna_flat_files,one_based=False)
-   genomicSeq = genomicSeq.lower()
+   def getFragmentSize(self,id):
+      return self.chromo_sizes[self.chromo_map[id]]
 
-   # check the obtained dna sequence
-   assert genomicSeq != '', 'load_genomic returned empty sequence!'
-   
-   # all entries other than a c g t are set to n
-   no_base = re.compile('[^acgt]')
-   genomicSeq = no_base.sub('n',genomicSeq)
 
-   if only_seq:
-      return genomicSeq
+   def getSpliceScores(self,id,strand,intervalBegin,intervalEnd):
+      """
+      Now we want to use interval_query to get the predicted splice scores trained
+      on the TAIR sequence and annotation.
+      """
+
+      size = intervalEnd-intervalBegin
+      assert size > 1, 'Error (getSpliceScores): interval size is less than 2!'
+
+      acc = size*[0.0]
+      don = size*[0.0]
+
+      interval_matrix = createIntArrayFromList([intervalBegin,intervalEnd])
+      pos_size = new_intp()
+      intp_assign(pos_size,1)
+
+      total_size = self.getFragmentSize(id)
+
+      acc_fn, don_fn = self.accessWrapper.get_splice_site_scores_fn(id,strand)
+
+      # fetch acceptor scores
+      acc = doIntervalQuery(id,strand,intervalBegin,intervalEnd,acc_fn,total_size)
+
+      # fetch donor scores
+      don = doIntervalQuery(id,strand,intervalBegin,intervalEnd,don_fn,total_size)
 
-   fn = 'chr%d.dna.flat' % chr
-   filename = os.path.join(dna_flat_files,fn)
-   total_size = get_flatfile_size(filename)
+      return acc, don
 
-   intervalBegin = genomicSeq_start-100
-   intervalEnd    = genomicSeq_stop+100
 
-   # shorten the intervalEnd to the maximal size of the genomic sequence if it
-   # would exceed this
-   if intervalEnd > total_size:
-      intervalEnd = total_size-1
+   def fetch_sequence(self,chromo,strand,genomicSeq_start,genomicSeq_stop):
+      filename = self.accessWrapper.get_genomic_fragment_fn(chromo)
+      genomicSeq = my_load_genomic(filename,'+',genomicSeq_start,genomicSeq_stop,one_based=False)
+      genomicSeq = genomicSeq.strip().lower()
 
-   # 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
+      # check the obtained dna sequence
+      assert genomicSeq != '', 'load_genomic returned empty sequence!'
       
-   if intervalBegin < NO_SCORE_WINDOW_SIZE:
-      check_assertions = False
+      # all entries other than a c g t are set to n
+      no_base = re.compile('[^acgt]')
+      genomicSeq = no_base.sub('n',genomicSeq)
 
-   currentAcc, currentDon = getSpliceScores(chr,strand,intervalBegin,intervalEnd,total_size)
+      return genomicSeq
 
-   currentAcc = currentAcc[100:-98]
-   currentAcc = currentAcc[1:]
-   currentDon = currentDon[100:-100]
 
-   length = len(genomicSeq)
+   def get_seq_and_scores(self,id,strand,genomicSeq_start,genomicSeq_stop,only_seq=False,perform_checks=True):
+      """
+      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 = 205
 
-   if strand == '+':
-      currentAcc = currentAcc[:length]
-      currentDon = currentDon+[-inf]*(length-len(currentDon))
+      chromo = self.chromo_map[id]
 
-      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')]
+      total_size = self.chromo_sizes[chromo]
 
-   # build reverse complement if on negative strand
-   if strand == '-':
-      currentAcc = currentAcc[:length]
-      currentAcc = currentAcc+[-inf]*(length-len(currentAcc))
+      print 'genomicSeq: ' + str(genomicSeq_start),str(genomicSeq_stop)
+
+      assert genomicSeq_start < genomicSeq_stop, pdb.set_trace()
+      assert genomicSeq_start >= 0
+      assert genomicSeq_stop <= total_size
+
+      genomicSeq = self.fetch_sequence(id,strand,genomicSeq_start,genomicSeq_stop)
+
+      if strand == '-':
+         genomicSeq = reverse_complement(genomicSeq)
+
+      if only_seq:
+         return genomicSeq
 
+      # shorten the intervalEnd to the maximal size of the genomic sequence if it would exceed this size
+      lookup_offset_begin  = 10
+      lookup_offset_end    = 10
+
+      intervalBegin  = max(genomicSeq_start-lookup_offset_begin,0)
+      if intervalBegin == 0:
+         lookup_offset_begin = 0
+         intervalBegin = genomicSeq_start
+         
+      intervalEnd    = min(genomicSeq_stop+lookup_offset_end,total_size-1)
+      if intervalEnd == total_size-1:
+         lookup_offset_end = 0
+         intervalEnd = genomicSeq_stop
+
+      currentAcc, currentDon = self.getSpliceScores(id,strand,intervalBegin,intervalEnd)
+
+      # remove the offset positions
+      if strand == '+':
+         currentAcc = currentAcc[lookup_offset_begin+2:]
+         currentDon = currentDon[lookup_offset_begin+1:]
+      elif strand == '-':
+         currentAcc = currentAcc[lookup_offset_begin:]
+         currentDon = [-inf]+currentDon
+         currentDon = currentDon[(lookup_offset_begin):]
+      else:
+         assert False
+
+      currentAcc = currentAcc[:len(genomicSeq)]
+      currentDon = currentDon[:len(genomicSeq)]
+
+      length = len(genomicSeq)
+      currentAcc = currentAcc+[-inf]*(length-len(currentAcc))
       currentDon = currentDon+[-inf]*(length-len(currentDon))
+      currentDon[-1] = -inf
+      
+      assert len(genomicSeq) == len(currentAcc) == len(currentDon), pdb.set_trace()
+
+      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']
+
+      for pos in ag_tuple_pos:
+         if pos+intervalBegin < NO_SCORE_WINDOW_SIZE and currentAcc[pos] == -inf:
+            currentAcc[pos] = 1e-6
+
+         if pos+intervalBegin > total_size - NO_SCORE_WINDOW_SIZE and currentAcc[pos] == -inf:
+            currentAcc[pos] = 1e-6
+
+      for pos in gt_tuple_pos:
+         if pos+intervalBegin < NO_SCORE_WINDOW_SIZE and currentDon[pos] == -inf:
+            currentDon[pos] = 1e-6
+
+         if pos+intervalBegin > total_size - NO_SCORE_WINDOW_SIZE and currentDon[pos] == -inf:
+            currentDon[pos] = 1e-6
+
+      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]
+
+      check_window_size = 30
+
+      if perform_checks and not ag_tuple_pos == acc_pos:
+         print 'ACC: Chromo/strand: %d/%s' % (id,strand)
+         print ag_tuple_pos[:check_window_size]
+         print acc_pos[:check_window_size]
+         print 
+         print ag_tuple_pos[-check_window_size:]
+         print acc_pos[-check_window_size:]
+         pdb.set_trace()
+
+      if perform_checks and not gt_tuple_pos == don_pos:
+         print 'DON: Chromo/strand: %d/%s' % (id,strand)
+         print gt_tuple_pos[:check_window_size]
+         print don_pos[:check_window_size]
+         print 
+         print gt_tuple_pos[-check_window_size:]
+         print don_pos[-check_window_size:]
+         pdb.set_trace()
+
+
+      return genomicSeq, currentAcc, currentDon
+
+
+
 
-      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
+##################################
+##################################
+##################################
+##################################
 
 
-def checks():
-   start  = 1001
-   stop   = start+1234
-   dna_flat_files    =  '/fml/ag-raetsch/share/projects/genomes/A_thaliana_best/genome/'
+   def _get_seq_and_scores(self,id,strand,genomicSeq_start,genomicSeq_stop,only_seq=False,perform_checks=True):
+      """
+      This function expects an interval, chromosome and strand information and
+      returns then the genomic sequence of this interval and the associated scores.
+      """
 
-   for chromo in range(1,6):
-         p_seq,_p_acc,_p_don = get_seq_and_scores(chromo,'+',start,stop,dna_flat_files)
-         n_seq,n_acc,n_don = get_seq_and_scores(chromo,'-',start,stop,dna_flat_files)
+      chromo = self.chromo_map[id]
 
-         print p_seq == reverse_complement(n_seq)
+      total_size = self.chromo_sizes[chromo]
 
+      if strand == '-':
+         s_start  = total_size - genomicSeq_stop
+         s_stop   = total_size - genomicSeq_start
+         genomicSeq_start = s_start
+         genomicSeq_stop  = s_stop
+
+
+      assert genomicSeq_start < genomicSeq_stop, pdb.set_trace()
+
+      genomicSeq = self.fetch_sequence(id,strand,genomicSeq_start,genomicSeq_stop)
+
+      if strand == '-':
+         genomicSeq = reverse_complement(genomicSeq)
+
+      if only_seq:
+         return genomicSeq
+
+      currentAcc,currentDon = self.get_scores(id,strand,genomicSeq_start,genomicSeq_stop,genomicSeq,total_size,perform_checks)
+
+      return genomicSeq, currentAcc, currentDon
+
+
+   def get_scores(self,id,strand,genomicSeq_start,genomicSeq_stop,genomicSeq,total_size,perform_checks):
+      # 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 = 205
+
+      seq_size = len(genomicSeq)
+      # shorten the intervalEnd to the maximal size of the genomic sequence if it would exceed this size
+      lookup_offset_begin  = 10
+      lookup_offset_end    = 10
+
+      intervalBegin  = max(genomicSeq_start-lookup_offset_begin,0)
+      if intervalBegin == 0:
+         lookup_offset_begin = 0
          
-if __name__ == '__main__':
-   checks()
+      intervalEnd    = min(genomicSeq_stop+lookup_offset_end,total_size-1)
+      if intervalEnd == total_size:
+         lookup_offset_end = 0
+
+      print 'interval is: %d %d' % (intervalBegin,intervalEnd)
+
+      # splice score indices are 1-based
+      currentAcc, currentDon = self.getSpliceScores(id,strand,intervalBegin+1,intervalEnd+1)
+
+      print 'original sizes: %d %d' % (len(currentAcc),len(currentDon))
+
+      # remove the offset positions
+      currentAcc = currentAcc[(lookup_offset_begin+1):(lookup_offset_begin+1+seq_size)]
+      currentDon = currentDon[lookup_offset_begin:lookup_offset_begin+seq_size]
+
+      print 'all sizes'
+      print len(genomicSeq)
+      print len(currentAcc)
+      print len(currentDon)
+
+      gt_tuple_pos = [p for p,e in enumerate(genomicSeq) if 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 e=='g' and genomicSeq[p-1]=='a' ]
+
+      for pos in ag_tuple_pos:
+         if pos+genomicSeq_start < NO_SCORE_WINDOW_SIZE and currentAcc[pos] == -inf:
+            currentAcc[pos] = 1e-6
+
+         if pos+genomicSeq_start > total_size - NO_SCORE_WINDOW_SIZE and currentAcc[pos] == -inf:
+            currentAcc[pos] = 1e-6
+
+      for pos in gt_tuple_pos:
+         if pos+genomicSeq_start < NO_SCORE_WINDOW_SIZE and currentDon[pos] == -inf:
+            currentDon[pos] = 1e-6
+
+         if pos+genomicSeq_start > total_size - NO_SCORE_WINDOW_SIZE and currentDon[pos] == -inf:
+            currentDon[pos] = 1e-6
+
+      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]
+
+      check_window_size = 30
+
+      if perform_checks and not ag_tuple_pos == acc_pos:
+         print 'ACC: Chromo/strand: %d/%s' % (id,strand)
+         print ag_tuple_pos[:check_window_size]
+         print acc_pos[:check_window_size]
+         print 
+         print ag_tuple_pos[-check_window_size:]
+         print acc_pos[-check_window_size:]
+         pdb.set_trace()
+
+      if perform_checks and not gt_tuple_pos == don_pos:
+         print 'DON: Chromo/strand: %d/%s' % (id,strand)
+         print gt_tuple_pos[:check_window_size]
+         print don_pos[:check_window_size]
+         print ''
+         print gt_tuple_pos[-check_window_size:]
+         print don_pos[-check_window_size:]
+         pdb.set_trace()
+
+      assert len(genomicSeq) == len(currentAcc) == len(currentDon), pdb.set_trace()
+
+      return currentAcc,currentDon