<commit_before>#(c) 2016 by Authors
#This file is a part of ABruijn program.
#Released under the BSD license (see LICENSE file)

"""
Quick and dirty alignment consensus
"""

import logging
from collections import defaultdict
from itertools import izip
import multiprocessing
import signal

from abruijn.alignment import shift_gaps, SynchronizedReader
import abruijn.config as config

logger = logging.getLogger()

class Profile:
    __slots__ = ("insertions", "matches", "nucl")

    def __init__(self):
        self.insertions = defaultdict(str)
        self.matches = defaultdict(int)
        self.nucl = "-"
        #self.length = length
        #self.insertions = [defaultdict(str) for _ in xrange(length)]
        #self.matches = [defaultdict(int) for _ in xrange(length)]
        #self.nucl = ["-" for _ in xrange(length)]


def _thread_worker(blasr_reader, contigs_info, results_queue,
                   error_queue):
    try:
        blasr_reader.init_reading()

        while not blasr_reader.is_eof():
            ctg_id, ctg_aln = blasr_reader.get_chunk()
            if ctg_id is None:
                break

            profile = _contig_profile(ctg_aln, contigs_info[ctg_id].length)
            sequence = _flattern_profile(profile)
            results_queue.put((ctg_id, sequence))

    except Exception as e:
        error_queue.put(e)


def get_consensus(alignment_path, contigs_info, num_proc):
    """
    Main function
    """
    blasr_reader = SynchronizedReader(alignment_path)
    manager = multiprocessing.Manager()
    results_queue = manager.Queue()
    error_queue = manager.Queue()

    #making sure the main process catches SIGINT
    orig_sigint = signal.signal(signal.SIGINT, signal.SIG_IGN)
    threads = []
    for _ in xrange(num_proc):
        threads.append(multiprocessing.Process(target=_thread_worker,
                                               args=(blasr_reader, contigs_info,
                                                     results_queue, error_queue)))
    signal.signal(signal.SIGINT, orig_sigint)

    for t in threads:
        t.start()
    try:
        for t in threads:
            t.join()
    except KeyboardInterrupt:
        for t in threads:
            t.terminate()

    if not error_queue.empty():
        raise error_queue.get()

    out_fasta = {}
    while not results_queue.empty():
        ctg_id, ctg_seq = results_queue.get()
        out_fasta[ctg_id] = ctg_seq

    return out_fasta


def _contig_profile(alignment, genome_len):
    """
    Computes alignment profile
    """
    MIN_ALIGNMENT = config.vals["min_alignment_length"]

    #profile = Profile(genome_len)
    profile = [Profile() for _ in xrange(genome_len)]
    for aln in alignment:
        if aln.err_rate > 0.5 or aln.trg_end - aln.trg_start < MIN_ALIGNMENT:
            continue

        qry_seq = shift_gaps(aln.trg_seq, aln.qry_seq)
        trg_seq = shift_gaps(qry_seq, aln.trg_seq)

        trg_pos = aln.trg_start
        for trg_nuc, qry_nuc in izip(trg_seq, qry_seq):
            if trg_nuc == "-":
                trg_pos -= 1
            if trg_pos >= genome_len:
                trg_pos -= genome_len
<commit_msg>do not output empty contig<commit_after>#(c) 2016 by Authors
#This file is a part of ABruijn program.
#Released under the BSD license (see LICENSE file)

"""
Quick and dirty alignment consensus
"""

import logging
from collections import defaultdict
from itertools import izip
import multiprocessing
import signal

from abruijn.alignment import shift_gaps, SynchronizedReader
import abruijn.config as config

logger = logging.getLogger()

class Profile:
    __slots__ = ("insertions", "matches", "nucl")

    def __init__(self):
        self.insertions = defaultdict(str)
        self.matches = defaultdict(int)
        self.nucl = "-"
        #self.length = length
        #self.insertions = [defaultdict(str) for _ in xrange(length)]
        #self.matches = [defaultdict(int) for _ in xrange(length)]
        #self.nucl = ["-" for _ in xrange(length)]


def _thread_worker(blasr_reader, contigs_info, results_queue,
                   error_queue):
    try:
        blasr_reader.init_reading()

        while not blasr_reader.is_eof():
            ctg_id, ctg_aln = blasr_reader.get_chunk()
            if ctg_id is None:
                break

            profile = _contig_profile(ctg_aln, contigs_info[ctg_id].length)
            sequence = _flattern_profile(profile)
            results_queue.put((ctg_id, sequence))

    except Exception as e:
        error_queue.put(e)


def get_consensus(alignment_path, contigs_info, num_proc):
    """
    Main function
    """
    blasr_reader = SynchronizedReader(alignment_path)
    manager = multiprocessing.Manager()
    results_queue = manager.Queue()
    error_queue = manager.Queue()

    #making sure the main process catches SIGINT
    orig_sigint = signal.signal(signal.SIGINT, signal.SIG_IGN)
    threads = []
    for _ in xrange(num_proc):
        threads.append(multiprocessing.Process(target=_thread_worker,
                                               args=(blasr_reader, contigs_info,
                                                     results_queue, error_queue)))
    signal.signal(signal.SIGINT, orig_sigint)

    for t in threads:
        t.start()
    try:
        for t in threads:
            t.join()
    except KeyboardInterrupt:
        for t in threads:
            t.terminate()

    if not error_queue.empty():
        raise error_queue.get()

    out_fasta = {}
    while not results_queue.empty():
        ctg_id, ctg_seq = results_queue.get()
        if len(ctg_seq) > 0:
            out_fasta[ctg_id] = ctg_seq

    return out_fasta


def _contig_profile(alignment, genome_len):
    """
    Computes alignment profile
    """
    MIN_ALIGNMENT = config.vals["min_alignment_length"]

    #profile = Profile(genome_len)
    profile = [Profile() for _ in xrange(genome_len)]
    for aln in alignment:
        if aln.err_rate > 0.5 or aln.trg_end - aln.trg_start < MIN_ALIGNMENT:
            continue

        qry_seq = shift_gaps(aln.trg_seq, aln.qry_seq)
        trg_seq = shift_gaps(qry_seq, aln.trg_seq)

        trg_pos = aln.trg_start
        for trg_nuc, qry_nuc in izip(trg_seq, qry_seq):
            if trg_nuc == "-":
                trg_pos -= 1
            if trg_pos >= genome_len:
                trg_pos -= genome_len

{'fraction_non_alphanumeric': 0.08999686422075887, 'fraction_numerical': 0.0023518344308560675, 'mean_word_length': 2.3804981452040277, 'pattern_counts': {'":': 2, '<': 5, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1315111', 'n_tokens_mistral': 2070, 'n_tokens_neox': 1999, 'n_words': 463}