python/computeSpliceAlign.py

   1 #!/usr/bin/env python
   2 # -*- coding: utf-8 -*-
   3
   4 from numpy.matlib import zeros
   5 import pdb
   6
   7 def  computeSpliceAlign(dna, exons):
   8    """
   9    Exonpos: Anfang Exon bis 1 rechts neben Ende Exon
  10    Daraus folgt: SpliceposAnfang ist gleich Ende Exon von Vorgaenger
  11    (ausser dem letzten und ersten)
  12     und Ende Splicestelle ist gleich Anfang Exon vom naechsten, da Splicestelleende notiert 1 rechts neben Splciestellenende
  13    bsp:
  14    cccccA1cccccE1cccccA2 ... AmcccccEmcccccccc
  15    cccccXXcccccA1cccccE1 ... Em-1cccXXcccccccc
  16    """
  17
  18    numberOfExons = (exons.shape)[0] # how many rows ?
  19    exonSizes = [-1]*numberOfExons
  20
  21    assert numberOfExons == 3
  22
  23    for idx in range(numberOfExons):
  24       exonSizes[idx] = exons[idx,1] - exons[idx,0]
  25
  26    sizeMatchmatrix = 6 # -acgtn
  27
  28    # SpliceAlign vector describes alignment:
  29    # 1:donorpos, 3:intron 2:acceptorpos, 0:exon, 4: dangling end
  30    SpliceAlign = []
  31
  32    if exons[0,0] > 0:
  33       SpliceAlign.extend([4]*(exons[0,0]))
  34
  35    for idx in range(numberOfExons):
  36       exonLength = exonSizes[idx]
  37       SpliceAlign.extend([0]*exonLength)
  38
  39       if idx < numberOfExons-1:
  40          intronLength = exons[idx+1,0] - exons[idx,1]
  41          SpliceAlign.extend([1]+[3]*(intronLength-2)+[2])
  42
  43    if len(dna) > exons[2,1]:
  44       SpliceAlign.extend([4]*(len(dna)-exons[2,1]))
  45
  46    assert len(SpliceAlign) == len(dna), pdb.set_trace()
  47
  48    # number of matches: just have to look at the underlying est
  49    # est = dna(find(SpliceAlign==0)) # exon nucleotides
  50    exonPos = [pos for pos,elem in enumerate(SpliceAlign) if elem == 0]
  51    est = [elem for pos,elem in enumerate(dna) if pos in exonPos]
  52
  53    #length_est = sum(exon_sizes) ;
  54    totalESTLength = 0
  55    for elem in exonSizes:
  56       totalESTLength += elem
  57
  58    assert totalESTLength == len(est)
  59
  60    # counts the occurences of a,c,g,t,n in this order
  61    numChar = [0]*sizeMatchmatrix
  62
  63    for elem in est:
  64       if elem == 'a':
  65          numChar[0] += 1
  66       if elem == 'c':
  67          numChar[1] += 1
  68       if elem == 'g':
  69          numChar[2] += 1
  70       if elem == 't':
  71          numChar[3] += 1
  72       if elem == 'n':
  73          numChar[4] += 1
  74
  75    totalNumChar = 0
  76    for idx in range(sizeMatchmatrix):
  77       totalNumChar += numChar[idx]
  78
  79    assert totalNumChar == len(est)
  80
  81    # writing in weight match matrix
  82    # matrix is saved columnwise
  83    trueWeightMatch = zeros((sizeMatchmatrix*sizeMatchmatrix,1)) # Scorematrix fuer Wahrheit
  84    for idx in range(1,sizeMatchmatrix):
  85       trueWeightMatch[(sizeMatchmatrix+1)*idx] = numChar[idx-1]
  86
  87    return SpliceAlign, trueWeightMatch