1 """PmfSimulation - peptide length distribution using GnuPlotPy
   2 """
   3 
   4 import unittest
   5 import Gnuplot
   6 """This module require ImitProFound, enzyme module made by yong27's company
   7 """
   8 import ImitProFound, enzyme
   9 
  10 class PeptideStat:
  11     def __init__(self,aTax='',aStep=20):
  12         self.tax=aTax
  13         self.seqIter=self.getFastaSeq()
  14         self.step=20
  15         self.peptides=[]
  16         self.histo=[0]*50000
  17     def getFastaSeq(self):
  18         return ImitProFound.getFastaSeqIter('fasta_'+self.tax)
  19     def digestEnzymeBy(self,anEnzyme):
  20         while 1:
  21             fastaSeq = self.seqIter.next()
  22             if fastaSeq is None: break
  23             self.peptides.append(ImitProFound.PeptideMap(fastaSeq,anEnzyme))
  24 
  25     def getHistoData(self):
  26         if not self.peptides:
  27             print "you must execute diegstEnzymeBy(anEnzyme)"
  28             sys.exit()
  29         for protein in self.peptides:
  30             for peptide in protein.peptideMasses:
  31                 self.histo[int(peptide/self.step)]+=1
  32         return self.histo
  33     def getHistoDataForGnuplot(self,aMethod):
  34         gpinput=[];i=0
  35         for each in aMethod():
  36             gpinput.append([i,each])
  37             i+=1
  38         return gpinput
  39     def getHistoDataByAaLen(self):
  40         if not self.peptides:
  41             print "you must execute diegstEnzymeBy(anEnzyme)"
  42             sys.exit()
  43         aaLen=[0]*5000
  44         for protein in self.peptides:
  45             for peptide in protein.peptSeq:
  46                 aaLen[len(peptide)]+=1
  47         return aaLen
  48 
  49 
  50 def main():
  51     pt=PeptideStat('Escherichia_coli')
  52     #pt=PeptideStat('All')
  53     pt.digestEnzymeBy(enzyme.Trypsin())
  54 
  55     g=Gnuplot.Gnuplot()
  56     g("set term png")
  57     g("set output 'plot.png'")
  58     g("set data style l")
  59     g("set xrange[0:50]")
  60     #g.plot(pt.getHistoDataForGnuplot(pt.getHistoData))
  61     g.plot(pt.getHistoDataForGnuplot(pt.getHistoDataByAaLen))
  62 
  63 class TestPeptideStat(unittest.TestCase):
  64     def test1(self):
  65         pt=PeptideStat('Escherichia_coli')
  66         pt.digestEnzymeBy(enzyme.Trypsin())
  67         print pt.getHistoData()
  68 
  69 if __name__=='__main__':
  70     #unittest.main(argv=('','-v'))
  71     main()