#!/usr/bin/python
# -*- coding:utf-8 -*-
'''Substitution mutation, such as A-->G.


by Wubin Qu <quwubin@gmail.com>,
Copyright @ 2010, All Rights Reserved.
'''

Author  = 'Wubin Qu  <quwubin@gmail.com> CZlab, BIRM, China'
Date    = 'Nov-05-2010 15:57:41'
License = 'GPL v3'
Version = '1.0'

import sys
import os
from optparse import OptionParser
from BifScripts import FastaFormatParser as FFP


def get_opt():
    '''Handle options'''
    usage = '''Usage: %prog [options]

    ./point_mutation.py -i infile.fa -p 263 -m A2T
    
    '''

    version = '%prog Version: ' + '%s [%s]' % (Version, Date)
    parser = OptionParser(usage=usage, version=version)
    parser.add_option('-i', '--infile', dest='infile', help='Input file name. [String]')
    parser.add_option('-o', '--outfile', dest='outfile', help='Output file name. [String]')
    parser.add_option('-p', '--position', dest='position', help='0-index position of the point mutation position. [Integer]')
    parser.add_option('-m', '--mutation', dest='mutation', help='Mutation detailes, legal input is, e.g. A2T. [Integer]')
    [options, args] = parser.parse_args()

    if len(args) > 1:
        parser.error('Incorrect argument, add" "-h" for help.')

    if not options.infile:
        parser.error('Input fasta-format file needed, add" "-h" for help.')

    if not options.position or int(options.position) < 0:
        parser.error('Point mutation position needed, add" "-h" for help.')

    if not options.mutation:
        parser.error('Mutation details needed, add" "-h" for help.')

    return options

def main ():
    '''Main'''
    options = get_opt()
    records = FFP.parse(open(options.infile))
    if len(records) > 1:
        print 'Currently, only one sequence support, %s has multiple sequences.' % options.infile
        exit()

    seq = records[0]['seq']
    id = records[0]['id']
    desc = records[0]['desc']
    position = int(options.position)
    if len(seq) < position + 1:
        print 'Error: Mutation position is longer than the sequence'
        exit()
    mutation_from, mutation_to = options.mutation.split('2')
    if seq[position] != mutation_from:
        print 'Error: nucleotide base in %s is not %s, actually is: %s' % (position, mutation_from, seq[position])
        exit()

    mutated_seq = seq[:position] + mutation_to + seq[position + 1 : ]
    output = '>%s %s%s%s' % (id, desc, os.linesep, mutated_seq)
    try:
        fh = open(options.outfile, 'w')
        fh.write(output)
        fh.close()
    except:
        print output

if __name__ == '__main__':
    main()