<filename>sam2aln.py
#! /usr/bin/env python

"""
!!!!!! ADAPTED FROM MiCall for one-time use in mixed-hcv project

Shipyard-style MiSeq pipeline, step 3
Takes SAM in CSV format as input.  Assumes that these data have been
through remap.py (but not strictly necessary!).  Merges paired-end reads
and outputs aligned sequences (with insertions and deletions, minus soft
clips).

Dependencies:
    bowtie2-build
    bowtie2-align
    settings.py
"""

import argparse
import collections
from csv import DictWriter
import itertools
import multiprocessing.forking
import multiprocessing.pool
import os
import re
import sys
from threading import Semaphore

#from micall.settings import max_prop_N, read_mapping_cutoff, sam2aln_q_cutoffs

from big_counter import BigCounter

max_prop_N = 0.2
read_mapping_cutoff = 0
sam2aln_q_cutoffs = [15]
nthreads = 4
READ_MERGE_CHUNK_SIZE = 100
READ_MERGE_BUFFER_SIZE = READ_MERGE_CHUNK_SIZE * nthreads * 10


def parseArgs():
    parser = argparse.ArgumentParser(
        description='Conversion of SAM data into aligned format.')
    parser.add_argument('samfile',
                        type=argparse.FileType('rU'),
                        help='<input> SAM output of bowtie2')
    parser.add_argument('aligned_csv',
                        type=argparse.FileType('w'),
                        help='<output> CSV containing cleaned and merged reads')
    parser.add_argument('insert_csv',
                        type=argparse.FileType('w'),
                        help='<output> CSV containing insertions relative to sample consensus')
    parser.add_argument('failed_csv',
                        type=argparse.FileType('w'),
                        help='<output> CSV containing reads that failed to merge')
    #parser.add_argument('-p', type=int, default=None, help='(optional) number of threads')
    
    return parser.parse_args()


cigar_re = re.compile('[0-9]+[MIDNSHPX=]')  # CIGAR token
gpfx = re.compile('^[-]+')  # length of gap prefix


# From https://github.com/pyinstaller/pyinstaller/wiki/Recipe-Multiprocessing
class _Popen(multiprocessing.forking.Popen):
    def __init__(self, *args, **kw):
        if hasattr(sys, 'frozen'):
            # We have to set original _MEIPASS2 value from sys._MEIPASS
            # to get --onefile mode working.
            # Last character is stripped in C-loader. We have to add
            # '/' or '\\' at the end.
            os.putenv('_MEIPASS2', sys._MEIPASS)  # @UndefinedVariable
        try:
            super(_Popen, self).__init__(*args, **kw)
        finally:
            if hasattr(sys, 'frozen'):
                # On some platforms (e.g. AIX) 'os.unsetenv()' is not
                # available. In those cases we cannot delete the variable
                # but only set it to the empty string. The bootloader
                # can handle this case.
                if hasattr(os, 'unsetenv'):
                    os.unsetenv('_MEIPASS2')
                else:
                    os.putenv('_MEIPASS2', '')


class Process(multiprocessing.Process):
    _Popen = _Popen

class Pool(multiprocessing.pool.Pool):
    Process = Process

def apply_cigar (cigar, seq, qual):
    """
    Use CIGAR string (Compact Idiosyncratic Gapped Alignment Report) in SAM data
    to apply soft clips, insertions, and deletions to the read sequence.
    Any insertions relative to the sample consensus sequence are discarded to
    enforce a strict pairwise alignment, and returned separately in a
    dict object.
    """
    newseq = ''
    newqual = ''
    insertions = {}
    tokens = cigar_re.findall(cigar)
    if len(tokens) == 0:
        return None, None, None
    # Account for removing soft clipped bases on left
    shift = 0
    if tokens[0].endswith('S'):
        shift = int(tokens[0][:-1])
    left = 0
    for token in tokens:
        length = int(token[:-1])
        # Matching sequence: carry it over
        if token[-1] == 'M':
            newseq += seq[left:(left+length)]
            newqual += qual[left:(left+length)]
            left += length
        # Deletion relative to reference: pad with gaps
        elif token[-1] == 'D':
            newseq += '-'*length
            newqual += ' '*length  # Assign fake placeholder score (Q=-1)
        # Insertion relative to reference: skip it (excise it)
        elif token[-1] == 'I':
            insertions[left] = (seq[left:(left+length)], qual[left:(left+length)])
            left += length
            continue
        # Soft clipping leaves the sequence in the SAM - so we should skip it
        elif token[-1] == 'S':
            left += length
            continue
        else:
            raise RuntimeError(
                "Unable to handle CIGAR token: {} - quitting".format(token))

    return shift, newseq, newqual, insertions


def merge_pairs (seq1, seq2, qual1, qual2, q_cutoff=10, minimum_q_delta=5):
    """
    Combine paired-end reads into a single sequence by managing discordant
    base calls on the basis of quality scores.
    """
    mseq = ''
    # force second read to be longest of the two
    if len(seq1) > len(seq2):
        seq1, seq2 = seq2, seq1
        qual1, qual2 = qual2, qual1

    for i, c2 in enumerate(seq2):
        q2 = ord(qual2[i])-33
        if i < len(seq1):
            c1 = seq1[i]
            q1 = ord(qual1[i])-33
            if c1 == '-' and c2 == '-':
                mseq += '-'
                continue
            if c1 == c2:  # Reads agree and at least one has sufficient confidence
                if q1 > q_cutoff or q2 > q_cutoff:
                    mseq += c1
                else:
                    mseq += 'N'  # neither base is confident
            else:
                if abs(q2 - q1) >= minimum_q_delta:
                    if q1 > max(q2, q_cutoff):
                        mseq += c1
                    elif q2 > max(q1, q_cutoff):
                        mseq += c2
                    else:
                        mseq += 'N'
                else:
                    mseq += 'N'  # cannot resolve between discordant bases
        else:
            # past end of read 1
            if c2 == '-' and q2 == 0:
                mseq += 'n'  # interval between reads
                continue
            mseq += c2 if (q2 > q_cutoff) else 'N'
    return mseq


def len_gap_prefix(s):
    hits = gpfx.findall(s)
    if hits:
        return len(hits[0])
    return 0

        
def is_first_read(flag):
    """
    Interpret bitwise flag from SAM field.
    Returns True or False indicating whether the read is the first read in a pair.
    """
    IS_FIRST_SEGMENT = 0x40
    return (int(flag) & IS_FIRST_SEGMENT) != 0


def matchmaker(samfile, semaphore=None):
    """
    An iterator that returns pairs of reads sharing a common qname from a remap CSV.
    Note that unpaired reads will be left in the cached_rows dictionary and
    discarded.
    :param samfile: open file handle to SAM output generated by bowtie2
    :param semaphore: acquired before each pair is generated, unless it is None
    :return: yields pairs of rows from DictReader corresponding to paired reads
    """
    #reader = DictReader(samfile)
    labels = ['qname', 'flag', 'rname', 'pos', 'mapq', 'cigar', 'rnext', 'pnext',
              'tlen', 'seq', 'qual']
    cached_rows = {}
    for line in samfile:
        if line.startswith('@'):
            continue  # skip header line
        if 'HCV' not in line:
            continue  # skip reads that mapped to another reference

        items = line.strip('\n').split('\t')
        row = dict(zip(labels, items[:11]))
        qname = row['qname']
        old_row = cached_rows.pop(qname, None)
        if old_row is None:
            cached_rows[qname] = row
        else:
            if semaphore is not None:
                semaphore.acquire()
            # current row should be the second read of the pair
            yield old_row, row


def parse_sam(rows):
    """ Merge two matched reads into a single aligned read.
    
    Also report insertions and failed merges.
    @param rows: tuple holding a pair of matched rows - forward and reverse reads
    @return: (refname, merged_seqs, insert_list, failed_list) where
        merged_seqs is {qcut: seq} the merged sequence for each cutoff level
        insert_list is [{'qname': query_name,
                         'fwd_rev': 'F' or 'R',
                         'refname': refname,
                         'pos': pos,
                         'insert': insertion_sequence,
                         'qual': insertion_quality_sequence}] insertions
        relative to the reference sequence.
        failed_list is [{'qname': query_name,
                         'qcut': qcut,
                         'seq1': seq1,
                         'qual1': qual1,
                         'seq2': seq2,
                         'qual2': qual2,
                         'prop_N': proportion_of_Ns,
                         'mseq': merged_sequence}] sequences that failed to
        merge.
    """
    row1, row2 = rows
    mseqs = {}
    failed_list = []
    insert_list = []
    rname = row1['rname']
    qname = row1['qname']
    cigar1 = row1['cigar']
    cigar2 = row2['cigar']

    # filtering criteria
    reason = None
    if cigar1 == '*':
        reason = 'R1 unmapped'
    if int(row1['mapq']) < read_mapping_cutoff:
        reason = 'R1 low mapq'

    if cigar2 == '*':
        reason = 'R2 unmapped'
    if int(row2['mapq']) < read_mapping_cutoff:
        reason = 'R2 low mapq'

    genotype1, genotype2 = None, None
    try:
        genotype1 = row1['rname'].split('-')[1][0]
        genotype2 = row2['rname'].split('-')[1][0]
    except:
        reason = 'discordant map'
        pass

    if genotype1 != genotype2:
        reason = 'map conflict'

    if reason:
        failed_list.append({'qname': qname,
                            'rname1': row1['rname'],
                            'rname2': row2['rname'],
                            'reason': reason})
    else:
        pos1 = int(row1['pos'])-1  # convert 1-index to 0-index
        _, seq1, qual1, inserts = apply_cigar(cigar1, row1['seq'], row1['qual'])
    
        # report insertions relative to sample consensus
        for left, (iseq, iqual) in inserts.iteritems():
            insert_list.append({'qname': qname,
                                'fwd_rev': 'F' if is_first_read(row1['flag']) else 'R',
                                'refname': rname,
                                'pos': pos1+left,
                                'insert': iseq,
                                'qual': iqual})
    
        seq1 = '-'*pos1 + seq1  # pad sequence on left
        qual1 = '!'*pos1 + qual1  # assign lowest quality to gap prefix so it does not override mate
    
    
        # now process the mate
        pos2 = int(row2['pos'])-1  # convert 1-index to 0-index
        _, seq2, qual2, inserts = apply_cigar(cigar2, row2['seq'], row2['qual'])
        for left, (iseq, iqual) in inserts.iteritems():
            insert_list.append({'qname': qname,
                                'fwd_rev': 'F' if is_first_read(row2['flag']) else 'R',
                                'refname': rname,
                                'pos': pos2+left,
                                'insert': iseq,
                                'qual': iqual})
        seq2 = '-'*pos2 + seq2
        qual2 = '!'*pos2 + qual2
    
        # merge reads
        for qcut in sam2aln_q_cutoffs:
            mseq = merge_pairs(seq1, seq2, qual1, qual2, qcut)
            prop_N = mseq.count('N') / float(len(mseq.strip('-')))
            if prop_N > max_prop_N:
                # fail read pair
                failed_list.append({'qname': qname,
                                    'reason': 'merge failure'})
                continue
            mseqs[qcut] = mseq

    return rname, mseqs, insert_list, failed_list

def parse_sam_in_threads(remap_csv, nthreads, semaphore):
    """ Call parse_sam() in multiple processes.
    
    Launch a multiprocessing pool, walk through the iterator, and then be sure
    to close the pool at the end.
    """
    pool = Pool(processes=nthreads)
    try:
        reads = pool.imap(parse_sam,
                          iterable=matchmaker(remap_csv, semaphore),
                          chunksize=100)
        for read in reads:
            yield read
    finally:
        pool.close()
        pool.join()


class CounterFactory:
    def __init__(self, aligned_file):
        self.aligned_file_name = getattr(aligned_file, 'name', None)

    def create_counter(self):
        if self.aligned_file_name is None:
            return collections.Counter()
        dirname = os.path.dirname(self.aligned_file_name)
        file_prefix = os.path.join(os.path.abspath(dirname),
                                   'merged_seq_counts')
        return BigCounter(file_prefix=file_prefix)


def sam2aln(remap_csv, aligned_csv, insert_csv, failed_csv, nthreads=None):
    # prepare outputs
    insert_fields =  ['qname', 'fwd_rev', 'refname', 'pos', 'insert', 'qual']
    insert_writer = DictWriter(insert_csv, insert_fields, lineterminator=os.linesep)
    insert_writer.writeheader()

    failed_fields =  ['qname', 'rname1', 'rname2', 'reason']
    failed_writer = DictWriter(failed_csv, failed_fields, lineterminator=os.linesep)
    failed_writer.writeheader()

    counter_factory = CounterFactory(aligned_csv)
    # We use a funky compound key, so that the BigCounter can write the keys
    # out to a CSV file. Counting separately for each region used hundreds
    # of counters.
    # aligned = {qcut: {region + '\t' + mseq: count}}
    aligned = collections.defaultdict(counter_factory.create_counter)
    if nthreads:
        semaphore = Semaphore(READ_MERGE_BUFFER_SIZE)
        iter = parse_sam_in_threads(remap_csv, nthreads, semaphore)
    else:
        semaphore = None
        iter = itertools.imap(parse_sam, matchmaker(remap_csv))

    for rname, mseqs, insert_list, failed_list in iter:
        if semaphore is not None:
            semaphore.release()
        for qcut, mseq in mseqs.iteritems():
            # collect identical merged sequences
            mseq_counter = aligned[qcut]
            mseq_counter[rname + '\t' + mseq] += 1

        # write out inserts to CSV
        insert_writer.writerows(insert_list)

        # write out failed read mergers to CSV
        failed_writer.writerows(failed_list)

    # write out merged sequences to file
    aligned_fields = ['refname', 'qcut', 'rank', 'count', 'offset', 'seq']
    aligned_writer = DictWriter(aligned_csv, aligned_fields, lineterminator=os.linesep)
    aligned_writer.writeheader()
    prev_rname = None
    rank = 0
    # for rname, region in sorted(aligned.iteritems()):
    for qcut, mseq_counter in aligned.iteritems():
        for key, count in mseq_counter.iteritems():
            rname, seq = key.split('\t')
            if rname != prev_rname:
                rank = 0
                prev_rname = rname
            else:
                rank += 1
            aligned_writer.writerow(dict(refname=rname,
                                         qcut=qcut,
                                         rank=rank,
                                         count=count,
                                         offset=len_gap_prefix(seq),
                                         seq=seq.strip('-')))


def main():
    args = parseArgs()
    sam2aln(remap_csv=args.samfile,
            aligned_csv=args.aligned_csv,
            insert_csv=args.insert_csv,
            failed_csv=args.failed_csv,
            nthreads=nthreads)

if __name__ == '__main__':
    main()
