#! /usr/bin/python

import sys,seq_convert,getopt

def usage():
  print "seq_mw.py [options] "
  print "          -h|--help == print this help text"
  print "          -b|--brief == print only brief output (total mass)"
  print "          -t 0|--terminal=0 == turn off the addition of one H20-equivalent for termini"
  print "          -a|--average == use natural abundance atomic masses"
  sys.exit(0)


brief_output = 0
average = 0
termini = 1
try:
  opts,args = getopt.getopt(sys.argv[1:],'hbat:',['help','brief','terminal=','average'])
except:
  sys.stderr.write("\n*********************************\n")
  sys.stderr.write("\n      Unknown options %s\n" % str(sys.argv[1:]))
  sys.stderr.write("\n*********************************\n\n")
  usage()

if len(opts) > 0:
  for o,a in opts:
    if o in ('-h', '--help'):
      usage()
    elif o in ('-b', '--brief'):
      brief_output = 1
    elif o in ('-a', '--average'):
      average = 1
    elif o in ('-t', '--terminal'):
      termini = int(a)

res_list = 'GAVLIPFWMSTCYNQDEKRH'
#    G   A   V   L   I   P   F   W   M   S   T   C   Y   N   Q   D   E   K   R   H
c = [2,  3,  5,  6,  6,  5,  9, 11,  5,  3,  4,  3,  9,  4,  5,  4,  5,  6,  6,  6]
n = [1,  1,  1,  1,  1,  1,  1,  2,  1,  1,  1,  1,  1,  2,  2,  1,  1,  2,  4,  3]
o = [1,  1,  1,  1,  1,  1,  1,  1,  1,  2,  2,  1,  2,  2,  2,  3,  3,  1,  1,  1]
s = [0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0]
h = [3,  5,  9, 11, 11,  7,  9, 10,  9,  5,  7,  5,  9,  6,  8,  5,  7, 12, 12,  7]

part_spec_vol = [0.64,0.74,0.86,0.90,0.90,0.76,0.77,0.74, 0.75,0.63,0.70,0.61,0.71,0.60,0.67,0.60,0.66,0.82,0.70,0.67]
avg_atomic_mass_table = {'C': 12.011, 'N': 14.007, 'O': 15.999, 'S': 32.064, 'H': 1.008}
isotopic_mass_table = {'C': 12., 'N': 14.00307, 'O': 15.99491, 'S': 31.97207, 'H': 1.007822}
if average:
  mass_table = avg_atomic_mass_table
else:
  mass_table = isotopic_mass_table

header,seq_in = seq_convert.readseq(sys.stdin,'1')
seq_in = seq_in.upper()

seq_convert.res_per_line = 60
seq_convert.out_format = '1'
seq_convert.print_seq_num = 1
seq_convert.write_out(seq_in,'')
print ''

mass = 0
num_el = 0
num_at = {'C': 0,'N':0,'O':0,'S':0,'H':0}
total_part_spec_vol = 0
aa_count = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]

num_res = len(seq_in)
for aa in seq_in:
  i = res_list.index(aa)
  num_at['C'] += c[i]
  num_at['N'] += n[i]
  num_at['O'] += o[i]
  num_at['S'] += s[i]
  num_at['H'] += h[i]
  aa_count[i] += 1
  total_part_spec_vol += part_spec_vol[i]

# fudge termini
if termini:
  num_at['O'] += 1
  num_at['H'] += 2

mass = (num_at['H']*mass_table['H'] + num_at['C']*mass_table['C'] + num_at['N']*mass_table['N'] 
       + num_at['O']*mass_table['O'] + num_at['S']*mass_table['S'])
num_el = num_at['H'] + num_at['C']*6 + num_at['N']*7 + num_at['O']*8 + num_at['S']*16
num_el += 10

print 'Mass: %10.2f' % mass
print 'Number of Electrons: %8d' % num_el
print 'Total partial specific volume: %6.4f '% (total_part_spec_vol/float(len(seq_in)))

print ''

if not brief_output:
  for k in num_at:
    print "Number of %s: %6d" % (k,num_at[k])

  print ''
  for i in range(len(res_list)):
    print "Number of %s: %5d  %5.1f%%" % (res_list[i],aa_count[i],(aa_count[i]*100/num_res))