"""

  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 3 of the License, or
  (at your option) any later version.
  
  Written (W) 2095-2010 Andre Kahles
  Copyright (C) 2009-2010 by Friedrich Miescher Laboratory, Tuebingen, Germany

  This script filters any given alignment by given criteria.
  
  For detailed usage information type:

    python filter_alignment.py 
"""

import sys
import re
import subprocess

class Feature(object):
    """Is an intron feature object"""


    def __init__(self, max_mm=80, feature_string=''):
        if feature_string == '':
            self.alignment_support = 0
            self.submission_support = 1
            self.mm_ex = dict()
            self.max_mm = max_mm + 1
        else:
            self.alignment_support = int(feature_string[0])
            self.submission_support = int(feature_string[1])
            self.mm_ex = dict()
            self.max_mm = max_mm + 1
            for _sl in feature_string[2:]:
                (key, value) = _sl.split(':')
                self.mm_ex[key] = int(value)

    def merge_features(self, feature_string):
        """Merges information in feature_string into current feature object"""
        self.alignment_support += int(feature_string[0])
        self.submission_support += int(feature_string[1])
        for _sl in feature_string[2:]:
            (key, value) = _sl.split(':')
            try:
                self.mm_ex[key] += int(value)
            except KeyError:
                self.mm_ex[key] = int(value)

    def add_mm_ex(self, ex, mm):
        """Adds mm ex information"""

        self.alignment_support += 1
        try:
            self.mm_ex[(ex*self.max_mm) + mm] += 1
        except KeyError:
            self.mm_ex[(ex*self.max_mm) + mm] = 1

    def get_feature_string(self):
        """Returns string with mm ex elements."""
        _line = (str(self.alignment_support) + '\t' + str(self.submission_support) + '\t')
        for key in self.mm_ex:
            _line += (str(key) + ':' + str(self.mm_ex[key]) + '\t')

        return _line[:-1]


    def get_submission_support(self):
        """Returns  submission support"""

        return int(self.submission_support)


def parse_options(argv):

    """Parses options from the command line """

    from optparse import OptionParser, OptionGroup

    parser = OptionParser()
    io_opts = OptionGroup(parser, 'I/O OPTIONS')
    io_opts.add_option('-a', '--alignment', dest='align', metavar='FILE', help='alignment file in sam format (- for stdin)', default='-')
    io_opts.add_option('-o', '--outfile', dest='outfile', metavar='FILE', help='outfile - default is tagged infile or stdout for stdin processing', default='-')
    io_opts.add_option('-b', '--bam_input', dest='bam_input', action='store_true', help='input has BAM format - does not work for STDIN', default=False)
    io_opts.add_option('-s', '--samtools', dest='samtools', metavar='PATH', help='if SAMtools is not in your PATH, provide the right path here (only neccessary for BAM input)', default='')
    filter_crit = OptionGroup(parser, 'FILTER CRITERIA')
    filter_crit.add_option('-R', '--ignore_multireads', dest='multireads', metavar='FILE', help='file containing the multireads to ignore', default='-')
    filter_crit.add_option('-M', '--max_intron_len', dest='max_intron_len', metavar='INT', type='int', help='maximal intron length [100000000]', default='100000000')
    filter_crit.add_option('-e', '--min_exon_len', dest='min_exon_len', metavar='INT', type='int', help='minimal exon length [0]', default=0)
    filter_crit.add_option('-X', '--max_mismatches', dest='max_mismatches', metavar='INT', type='int', help='maximum number of allowed mismathes [10000]', default=10000)
    filter_crit.add_option('-V', '--variants', dest='variants', action='store_true', help='use variant tags XM and XG for mismatch filtering [off]', default=False)
    filter_crit.add_option('-i', '--intron_features', dest='intron_features', metavar='FILE', help='intron features file - only spliced reads present in this file are kept', default='-')
    filter_crit.add_option('-c', '--clip_filter', dest='clip_filter', action='store_true', help='filters clipped reads', default=False)
    filter_hand = OptionGroup(parser, 'FILTER HANDLING')
    filter_hand.add_option('-d', '--del_worse', dest='del_worse', action='store_true', \
        help='all alignments of the same read with more mismatches are deleted as well. File must be sorted by read ID!', default=False)
    filter_hand.add_option('-n', '--no_suboptimal', dest='no_suboptimal', action='store_true', \
        help='if suboptimal spliced alignments occur and a batter alignment is available, keep only the best alignment. File must be sorted by read id and flag', default=False)
    filter_hand.add_option('-w', '--window', dest='window', metavar='INT', type='int', help='size of overlap-window to count suboptimal alignments as one stratum [2]', default=2)
    others = OptionGroup(parser, 'OTHERS')
    others.add_option('-v', '--verbose', dest='verbose', action='store_true', help='verbosity [off]', default=False)
    parser.add_option_group(io_opts)
    parser.add_option_group(filter_crit)
    parser.add_option_group(filter_hand)
    parser.add_option_group(others)

    (options, args) = parser.parse_args()
    
    if len(argv) < 2:
        parser.print_help()
        sys.exit(2)

    return options

def filter_lines_by_criteria(curr_lines, multireads, options, filter_counter):
    """Filters the list of lines by given criteria"""

    min_mm = 100000
    _curr_lines = []

    for sl in curr_lines:
        mm = -1 
        if options.del_worse:
            try:
                for opt in sl[11:]:
                    if options.variants and opt[:3] in ['XG:', 'XM:']:
                        if mm == -1:
                            mm = int(opt[5:])
                        else:
                            mm += int(opt[5:])
                    if not options.variants and opt[:3] == 'NM:':
                        mm = int(opt[5:])
                        break
                if mm == -1:
                    print >> sys.stderr, 'No mismatch information available or read string missing in %s' % options.align
                    sys.exit(1)
            except IndexError:
                print >> sys.stderr, 'No mismatch information available or read string missing in %s' % options.align
                sys.exit(1)

        if sl[5].find('N') > -1 and options.intron_features != '-':
            cont_flag = True
            (op, size) = (re.split('[0-9]*', sl[5])[1:], re.split('[^0-9]', sl[5])[:-1])
            size = [int(i) for i in size]
            offset = 0
            for o in range(len(op)):
                if op[o] in ['M', 'D']:
                    offset += size[o]
                if op[o] == 'N':
                    start = int(sl[3]) - 1 + offset
                    stop = start + size[o]
                    if options.features.has_key((sl[2], str(start), str(stop))):
                        cont_flag = False
                        break
            if cont_flag:
                filter_counter += 1
                if mm != -1:
                    min_mm = min(mm, min_mm)
                continue

        if (sl[0], int(sl[1]) & 196) in multireads:
            filter_counter += 1
            if mm != -1:
                min_mm = min(mm, min_mm)
            continue

        if options.max_mismatches < 10000:
            cont_flag = False
            if options.del_worse and mm > options.max_mismatches:
                filter_counter += 1
                if mm != -1:
                    min_mm = min(mm, min_mm)
                continue
            else:
                try:
                    _mm = 0
                    for opt in sl[11:]:
                        if options.variants and opt[:3] in ['XM:', 'XG:']:
                            _mm += int(opt[5:])        
                        elif (not options.variants) and (opt[:3] == 'NM:'):
                            _mm += int(opt[5:])
                        if _mm > options.max_mismatches:
                            cont_flag = True
                            break
                    if cont_flag:
                        filter_counter += 1
                        if mm != -1:
                            min_mm = min(mm, min_mm)