import os, subprocess

import itertools

import operator

from collections import defaultdict, OrderedDict

import errno

# cPickle is a faster version of pickle that isn't installed in python3

# inserted try statement just in case

try:

   import cPickle as pickle

except:

   import pickle

from io import BytesIO

import numpy as np

import re

import shutil

import gzip

# product: product(A, B) returns the same as ((x,y) for x in A for y in B).

# combination: Return r length subsequences of elements from the input iterable.

from itertools import product, combinations

import time

import yaml

import tempfile

import string

from contextlib import contextmanager

# -----------------------

#

# Helper functions

#

# -----------------------

def string_hamming_distance(str1, str2):

    """

    Fast hamming distance over 2 strings known to be of same length.

    In information theory, the Hamming distance between two strings of equal 

    length is the number of positions at which the corresponding symbols 

    are different.

    eg "karolin" and "kathrin" is 3.

    """

    return sum(itertools.imap(operator.ne, str1, str2))

___tbl = {'A':'T', 'T':'A', 'C':'G', 'G':'C', 'N':'N'}

def rev_comp(seq):

    return ''.join(___tbl[s] for s in seq[::-1])

def to_fastq(name, seq, qual):

    """

    Return string that can be written to fastQ file

    """

    return '@'+name+'\n'+seq+'\n+\n'+qual+'\n'

def to_fastq_lines(bc, umi, seq, qual, read_name=''):

    """

    Return string that can be written to fastQ file

    """

    reformated_name = read_name.replace(':', '_')

    name = '%s:%s:%s' % (bc, umi, reformated_name)

    return to_fastq(name, seq, qual)

def from_fastq(handle):

    while True:

        name = next(handle).rstrip()[1:] #Read name

        seq = next(handle).rstrip() #Read seq

        next(handle) #+ line

        qual = next(handle).rstrip() #Read qual

        if not name or not seq or not qual:

            break

        yield name, seq, qual

def seq_neighborhood(seq, n_subs=1):

    """

    Given a sequence, yield all sequences within n_subs substitutions of 

    that sequence by looping through each combination of base pairs within

    each combination of positions.

    """

    for positions in combinations(range(len(seq)), n_subs):

    # yields all unique combinations of indices for n_subs mutations

        for subs in product(*("ATGCN",)*n_subs):

        # yields all combinations of possible nucleotides for strings of length

        # n_subs

            seq_copy = list(seq)

            for p, s in zip(positions, subs):

                seq_copy[p] = s

            yield ''.join(seq_copy)

def build_barcode_neighborhoods(barcode_file, expect_reverse_complement=True):

    """

    Given a set of barcodes, produce sequences which can unambiguously be

    mapped to these barcodes, within 2 substitutions. If a sequence maps to 

    multiple barcodes, get rid of it. However, if a sequences maps to a bc1 with 

    1change and another with 2changes, keep the 1change mapping.

    """

    # contains all mutants that map uniquely to a barcode

    clean_mapping = dict()

    # contain single or double mutants 

    mapping1 = defaultdict(set)

    mapping2 = defaultdict(set)

    #Build the full neighborhood and iterate through barcodes

    with open(barcode_file, 'rU') as f:

        # iterate through each barcode (rstrip cleans string of whitespace)

        for line in f:

            barcode = line.rstrip()

            if expect_reverse_complement:

                barcode = rev_comp(line.rstrip())

            # each barcode obviously maps to itself uniquely

            clean_mapping[barcode] = barcode

            # for each possible mutated form of a given barcode, either add

            # the origin barcode into the set corresponding to that mutant or 

            # create a new entry for a mutant not already in mapping1

            # eg: barcodes CATG and CCTG would be in the set for mutant CTTG

            # but only barcode CATG could generate mutant CANG

            for n in seq_neighborhood(barcode, 1):

                mapping1[n].add(barcode)

            # same as above but with double mutants

            for n in seq_neighborhood(barcode, 2):

                mapping2[n].add(barcode)   

    # take all single-mutants and find those that could only have come from one

    # specific barcode

    for k, v in mapping1.items():

        if k not in clean_mapping:

            if len(v) == 1:

                clean_mapping[k] = list(v)[0]

    for k, v in mapping2.items():

        if k not in clean_mapping:

            if len(v) == 1:

                clean_mapping[k] = list(v)[0]

    del mapping1

    del mapping2

    return clean_mapping

def check_dir(path):

    """

    Checks if directory already exists or not and creates it if it doesn't

    """

    try:

        os.makedirs(path)

    except OSError as exception:

        if exception.errno != errno.EEXIST:

            raise

def print_to_stderr(msg, newline=True):

    """

    Wrapper to eventually write to stderr

    """

    sys.stderr.write(str(msg))

    if newline:

        sys.stderr.write('\n')

def worker_filter(iterable, worker_index, total_workers):

    return (p for i,p in enumerate(iterable) if (i-worker_index)%total_workers==0)

class FIFO():

    """

    A context manager for a named pipe.

    """

    def __init__(self, filename="", suffix="", prefix="tmp_fifo_dir", dir=None):

        if filename:

            self.filename = filename

        else:

            self.tmpdir = tempfile.mkdtemp(suffix=suffix, prefix=prefix, dir=dir)

            self.filename = os.path.join(self.tmpdir, 'fifo')

    def __enter__(self):

        if os.path.exists(self.filename):

            os.unlink(self.filename)

        os.mkfifo(self.filename)

        return self

    def __exit__(self, type, value, traceback):

        os.remove(self.filename)

        if hasattr(self, 'tmpdir'):

            shutil.rmtree(self.tmpdir)

# -----------------------

#

# Core objects

#

# -----------------------

class IndropsProject():

    def __init__(self, project_yaml_file_handle, read_only=False):

        self.yaml = yaml.load(project_yaml_file_handle)

        self.name = self.yaml['project_name']

        self.project_dir = self.yaml['project_dir']

        self.libraries = OrderedDict()

        self.runs = OrderedDict()

        self.read_only = read_only

        for run in self.yaml['sequencing_runs']:

            """

            After filtering, each sequencing run generates between 1 ... X files with filtered reads.

               X = (N x M)

             - N: The run is often split into several files (a typical NextSeq run is split into L001,

                  L002, L003, L004 which match different lanes, but this can also be done artificially.

             - M: The same run might contain several libraries. The demultiplexing can be handled by the script (or externally).

                  If demultiplexing is done externally, there will be a different .fastq file for each library.

            """

            version = run['version']

            filtered_filename = '{library_name}_{run_name}'

            if run['version'] == 'v3':

                filtered_filename += '_{library_index}'

            # Prepare to iterate over run split into several files

            if 'split_affixes' in run:

                filtered_filename += '_{split_affix}'

                split_affixes = run['split_affixes']

            else:

                split_affixes = ['']

            filtered_filename += '.fastq'

            # Prepare to iterate over libraries

            if 'libraries' in run:

                run_libraries = run['libraries']

            elif 'library_name' in run:

                run_libraries = [{'library_name' : run['library_name'], 'library_prefix':''}]

            else:

                raise Exception('No library name or libraries specified.')

            if run['version']=='v1' or run['version']=='v2':

                for affix in split_affixes:

                    for lib in run_libraries:

                        lib_name = lib['library_name']

                        if lib_name not in self.libraries:

                            self.libraries[lib_name] = IndropsLibrary(name=lib_name, project=self, version=run['version'])

                        else:

                            assert self.libraries[lib_name].version == run['version']

                        if version == 'v1':

                            metaread_filename = os.path.join(run['dir'],run['fastq_path'].format(split_affix=affix, read='R1', library_prefix=lib['library_prefix']))

                            bioread_filename = os.path.join(run['dir'],run['fastq_path'].format(split_affix=affix, read='R2', library_prefix=lib['library_prefix']))

                        elif version == 'v2':

                            metaread_filename  = os.path.join(run['dir'],run['fastq_path'].format(split_affix=affix, read='R2', library_prefix=lib['library_prefix']))

                            bioread_filename = os.path.join(run['dir'],run['fastq_path'].format(split_affix=affix, read='R1', library_prefix=lib['library_prefix']))

                        filtered_part_filename = filtered_filename.format(run_name=run['name'], split_affix=affix, library_name=lib_name)

                        filtered_part_path = os.path.join(self.project_dir, lib_name, 'filtered_parts', filtered_part_filename)

                        part = V1V2Filtering(filtered_fastq_filename=filtered_part_path,

                            project=self, 

                            bioread_filename=bioread_filename,

                            metaread_filename=metaread_filename,

                            run_name=run['name'],

                            library_name=lib_name,

                            part_name=affix)

                        if run['name'] not in self.runs:

                            self.runs[run['name']] = []

                        self.runs[run['name']].append(part)

                        self.libraries[lib_name].parts.append(part)

            elif run['version'] == 'v3' or run['version'] == 'v3-miseq':

                for affix in split_affixes:

                    filtered_part_filename = filtered_filename.format(run_name=run['name'], split_affix=affix,

                        library_name='{library_name}', library_index='{library_index}')

                    part_filename = os.path.join(self.project_dir, '{library_name}', 'filtered_parts', filtered_part_filename)

                    input_filename = os.path.join(run['dir'], run['fastq_path'].format(split_affix=affix, read='{read}'))

                    part = V3Demultiplexer(run['libraries'], project=self, part_filename=part_filename, input_filename=input_filename, run_name=run['name'], part_name=affix,

                        run_version_details=run['version'])

                    if run['name'] not in self.runs:

                        self.runs[run['name']] = []

                    self.runs[run['name']].append(part)

                    for lib in run_libraries:

                        lib_name = lib['library_name']

                        lib_index = lib['library_index']

                        if lib_name not in self.libraries:

                            self.libraries[lib_name] = IndropsLibrary(name=lib_name, project=self, version=run['version'])

                        self.libraries[lib_name].parts.append(part.libraries[lib_index])

    @property

    def paths(self):

        if not hasattr(self, '_paths'):

            script_dir = os.path.dirname(os.path.realpath(__file__))

            #Read defaults

            with open(os.path.join(script_dir, 'default_parameters.yaml'), 'r') as f:

                paths = yaml.load(f)['paths']

            # Update with user provided values

            paths.update(self.yaml['paths'])

            paths['python'] = os.path.join(paths['python_dir'], 'python')

            paths['java'] = os.path.join(paths['java_dir'], 'java')

            paths['bowtie'] = os.path.join(paths['bowtie_dir'], 'bowtie')

            paths['samtools'] = os.path.join(paths['samtools_dir'], 'samtools')

            paths['trimmomatic_jar'] = os.path.join(script_dir, 'bins', 'trimmomatic-0.33.jar')

            paths['rsem_tbam2gbam'] = os.path.join(paths['rsem_dir'], 'rsem-tbam2gbam')

            paths['rsem_prepare_reference'] = os.path.join(paths['rsem_dir'], 'rsem-prepare-reference')

            self._paths = type('Paths_anonymous_object',(object,),paths)()

            self._paths.trim_polyA_and_filter_low_complexity_reads_py = os.path.join(script_dir, 'trim_polyA_and_filter_low_complexity_reads.py')

            self._paths.quantify_umifm_from_alignments_py = os.path.join(script_dir, 'quantify_umifm_from_alignments.py')

            self._paths.count_barcode_distribution_py = os.path.join(script_dir, 'count_barcode_distribution.py')

            self._paths.gel_barcode1_list = os.path.join(script_dir, 'ref/barcode_lists/gel_barcode1_list.txt')

            self._paths.gel_barcode2_list = os.path.join(script_dir, 'ref/barcode_lists/gel_barcode2_list.txt')

        return self._paths

    @property

    def parameters(self):

        if not hasattr(self, '_parameters'):

            #Read defaults

            with open(os.path.join(os.path.dirname(os.path.realpath(__file__)), 'default_parameters.yaml'), 'r') as f:

                self._parameters = yaml.load(f)['parameters']

            # Update with user provided values

            if 'parameters' in self.yaml:

                for k, d in self.yaml['parameters'].items():

                    self._parameters[k].update(d)

        return self._parameters

    @property

    def gel_barcode1_revcomp_list_neighborhood(self):

        if not hasattr(self, '_gel_barcode1_list_neighborhood'):

            self._gel_barcode1_revcomp_list_neighborhood = build_barcode_neighborhoods(self.paths.gel_barcode1_list, True)

        return self._gel_barcode1_revcomp_list_neighborhood

    @property

    def gel_barcode2_revcomp_list_neighborhood(self):

        if not hasattr(self, '_gel_barcode2_revcomp_list_neighborhood'):

            self._gel_barcode2_revcomp_list_neighborhood = build_barcode_neighborhoods(self.paths.gel_barcode2_list, True)

        return self._gel_barcode2_revcomp_list_neighborhood

    @property

    def gel_barcode2_list_neighborhood(self):

        if not hasattr(self, '_gel_barcode2_list_neighborhood'):

            self._gel_barcode2_list_neighborhood = build_barcode_neighborhoods(self.paths.gel_barcode2_list, False)

        return self._gel_barcode2_list_neighborhood

    @property

    def stable_barcode_names(self):

        if not hasattr(self, '_stable_barcode_names'):

            with open(self.paths.gel_barcode1_list) as f:

                rev_bc1s = [rev_comp(line.rstrip()) for line in f]

            with open(self.paths.gel_barcode2_list) as f:

                bc2s = [line.rstrip() for line in f]

                rev_bc2s = [rev_comp(bc2) for bc2 in bc2s]

            # V1, V2 names:

            v1v2_names = {}

            barcode_iter = product(rev_bc1s, rev_bc2s)

            name_iter = product(string.ascii_uppercase, repeat=4)

            for barcode, name in zip(barcode_iter, name_iter):

                v1v2_names['-'.join(barcode)] = 'bc' + ''.join(name)

            # V3 names:

            v3_names = {}

            barcode_iter = product(bc2s, rev_bc2s)

            name_iter = product(string.ascii_uppercase, repeat=4)

            for barcode, name in zip(barcode_iter, name_iter):

                v3_names['-'.join(barcode)] = 'bc' + ''.join(name)

            self._stable_barcode_names = {

                'v1' : v1v2_names,

                'v2' : v1v2_names,

                'v3': v3_names,

                'v3-miseq':v3_names,

            }

        return self._stable_barcode_names

    def project_check_dir(self, path):

        if not self.read_only:

            check_dir(path)

    def filter_gtf(self, gzipped_transcriptome_gtf, gtf_with_genenames_in_transcript_id):

        # A small number of gene are flagged as having two different biotypes.

        gene_biotype_dict = defaultdict(set)

        # Read through GTF file once to get all gene names

        for line in subprocess.Popen(["gzip", "--stdout", "-d", gzipped_transcriptome_gtf], stdout=subprocess.PIPE).stdout:

            # Skip non-gene feature lines.

            if '\tgene\t' not in line:

                continue

            gene_biotype_match = re.search(r'gene_biotype \"(.*?)\";', line)

            gene_name_match = re.search(r'gene_name \"(.*?)\";', line)

            if gene_name_match and gene_biotype_match:

                gene_name = gene_name_match.group(1)

                gene_biotype = gene_biotype_match.group(1)

                # Record biotype.

                gene_biotype_dict[gene_name].add(gene_biotype)

        # Detect read-through genes by name. Name must be a fusion of two other gene names 'G1-G2'.

        readthrough_genes = set()

        for gene in gene_biotype_dict.keys():

            if '-' in gene and len(gene.split('-')) == 2:

                g1, g2 = gene.split('-')

                if g1 in gene_biotype_dict and g2 in gene_biotype_dict:

                    readthrough_genes.add(gene)

        # Detect pseudogenes: genes where all associated biotypes have 'pseudogene' in name

        pseudogenes = set()

        for gene, biotypes in gene_biotype_dict.items():

            if all('pseudogene' in b for b in biotypes):

                pseudogenes.add(gene)

        all_genes = set(gene_biotype_dict.keys())

        valid_genes = all_genes.difference(pseudogenes).difference(readthrough_genes)

        transcripts_counter = defaultdict(int)

        # Go through GTF file again, annotating each transcript_id with the gene and outputting to a new GTF file.

        output_gtf = open(gtf_with_genenames_in_transcript_id, 'w')

        for line in subprocess.Popen(["gzip", "--stdout", "-d", gzipped_transcriptome_gtf], stdout=subprocess.PIPE).stdout:

            # Skip non-transcript feature lines.

            if 'transcript_id' not in line:

                continue

            gene_name_match = re.search(r'gene_name \"(.*?)\";', line)

            if gene_name_match:

                gene_name = gene_name_match.group(1)

                if gene_name in valid_genes:

                    # An unusual edgecase in the GTF for Danio Rerio rel89

                    if ' ' in gene_name:

                        gene_name = gene_name.replace(' ', '_')

                    out_line = re.sub(r'(?<=transcript_id ")(.*?)(?=";)', r'\1|'+gene_name, line)

                    output_gtf.write(out_line)

                    if '\ttranscript\t' in line:

                        transcripts_counter['valid'] += 1

                elif gene_name in pseudogenes and '\ttranscript\t' in line:

                    transcripts_counter['pseudogenes'] += 1

                elif gene_name in readthrough_genes and '\ttranscript\t' in line:

                    transcripts_counter['readthrough_genes'] += 1

        output_gtf.close()

        print_to_stderr('Filtered GTF contains %d transcripts (%d genes)' % (transcripts_counter['valid'], len(valid_genes)))

        print_to_stderr('   - ignored %d transcripts from %d pseudogenes)' % (transcripts_counter['pseudogenes'], len(pseudogenes)))

        print_to_stderr('   - ignored %d read-through transcripts (%d genes)' % (transcripts_counter['readthrough_genes'], len(readthrough_genes)))

    def build_transcriptome(self, gzipped_genome_softmasked_fasta_filename, gzipped_transcriptome_gtf,

            mode='strict'):

        import pyfasta

        index_dir = os.path.dirname(self.paths.bowtie_index)

        self.project_check_dir(index_dir)

        genome_filename = os.path.join(index_dir, '.'.join(gzipped_genome_softmasked_fasta_filename.split('.')[:-1]))

        gtf_filename = os.path.join(index_dir, gzipped_transcriptome_gtf.split('/')[-1])

        gtf_prefix = '.'.join(gtf_filename.split('.')[:-2])

        # gtf_with_genenames_in_transcript_id = gtf_prefix + '.annotated.gtf'

        gtf_with_genenames_in_transcript_id = self.paths.bowtie_index + '.gtf'

        print_to_stderr('Filtering GTF')

        self.filter_gtf(gzipped_transcriptome_gtf, gtf_with_genenames_in_transcript_id)

        # accepted_gene_biotypes_for_NA_transcripts = set(["protein_coding","IG_V_gene","IG_J_gene","TR_J_gene","TR_D_gene","TR_V_gene","IG_C_gene","IG_D_gene","TR_C_gene"])

        # tsl1_or_tsl2_strings = ['transcript_support_level "1"', 'transcript_support_level "1 ', 'transcript_support_level "2"', 'transcript_support_level "2 ']

        # tsl_NA =  'transcript_support_level "NA'

        # def filter_ensembl_transcript(transcript_line):

        #     line_valid_for_output = False

        #     if mode == 'strict':

        #         for string in tsl1_or_tsl2_strings:

        #             if string in line:

        #                 line_valid_for_output = True

        #                 break

        #         if tsl_NA in line:

        #             gene_biotype = re.search(r'gene_biotype \"(.*?)\";', line)

        #             if gene_biotype and gene_biotype.group(1) in accepted_gene_biotypes_for_NA_transcripts:

        #                 line_valid_for_output = True

        #         return line_valid_for_output

        #     elif mode == 'all_ensembl':

        #         line_valid_for_output = True

        #         return line_valid_for_output

        # print_to_stderr('Filtering GTF')

        # output_gtf = open(gtf_with_genenames_in_transcript_id, 'w')

        # for line in subprocess.Popen(["gzip", "--stdout", "-d", gzipped_transcriptome_gtf], stdout=subprocess.PIPE).stdout:

        #     if 'transcript_id' not in line:

        #         continue

        #     if filter_ensembl_transcript(line):

        #         gene_name = re.search(r'gene_name \"(.*?)\";', line)

        #         if gene_name:

        #             gene_name = gene_name.group(1)

        #             out_line = re.sub(r'(?<=transcript_id ")(.*?)(?=";)', r'\1|'+gene_name, line)

        #             output_gtf.write(out_line)

        # output_gtf.close()

        print_to_stderr('Gunzipping Genome')

        p_gzip = subprocess.Popen(["gzip", "-dfc", gzipped_genome_softmasked_fasta_filename], stdout=open(genome_filename, 'wb'))

        if p_gzip.wait() != 0:

            raise Exception(" Error in rsem-prepare reference ")

        p_rsem = subprocess.Popen([self.paths.rsem_prepare_reference, '--bowtie', '--bowtie-path', self.paths.bowtie_dir,

                            '--gtf', gtf_with_genenames_in_transcript_id, 

                            '--polyA', '--polyA-length', '5', genome_filename, self.paths.bowtie_index])

        if p_rsem.wait() != 0:

            raise Exception(" Error in rsem-prepare reference ")

        print_to_stderr('Finding soft masked regions in transcriptome')

        transcripts_fasta = pyfasta.Fasta(self.paths.bowtie_index + '.transcripts.fa')

        soft_mask = {}

        for tx, seq in transcripts_fasta.items():

            seq = str(seq)

            soft_mask[tx] = set((m.start(), m.end()) for m in re.finditer(r'[atcgn]+', seq))

        with open(self.paths.bowtie_index + '.soft_masked_regions.pickle', 'w') as out:

            pickle.dump(soft_mask, out)

class IndropsLibrary():

    def __init__(self, name='', project=None, version=''):

        self.project = project

        self.name = name

        self.parts = []

        self.version = version

        self.paths = {}

        for lib_dir in ['filtered_parts', 'quant_dir']:

            dir_path = os.path.join(self.project.project_dir, self.name, lib_dir)

            self.project.project_check_dir(dir_path)

            self.paths[lib_dir] = dir_path

        self.paths = type('Paths_anonymous_object',(object,),self.paths)()

        self.paths.abundant_barcodes_names_filename = os.path.join(self.project.project_dir, self.name, 'abundant_barcodes.pickle')

        self.paths.filtering_statistics_filename = os.path.join(self.project.project_dir, self.name, self.name+'.filtering_stats.csv')

        self.paths.barcode_abundance_histogram_filename = os.path.join(self.project.project_dir, self.name, self.name+'.barcode_abundance.png')

        self.paths.barcode_abundance_by_barcode_filename = os.path.join(self.project.project_dir, self.name, self.name+'.barcode_abundance_by_barcode.png')

        self.paths.missing_quants_filename = os.path.join(self.project.project_dir, self.name, self.name+'.missing_barcodes.pickle')

    @property

    def barcode_counts(self):

        if not hasattr(self, '_barcode_counts'):

            self._barcode_counts = defaultdict(int)

            for part in self.parts:

                for k, v in part.part_barcode_counts.items():

                    self._barcode_counts[k] += v

        return self._barcode_counts

    @property

    def abundant_barcodes(self):

        if not hasattr(self, '_abundant_barcodes'):

            with open(self.paths.abundant_barcodes_names_filename) as f:

                self._abundant_barcodes = pickle.load(f)

        return self._abundant_barcodes

    def sorted_barcode_names(self, min_reads=0, max_reads=10**10):

        return [name for bc,(name,abun) in sorted(self.abundant_barcodes.items(), key=lambda i:-i[1][1]) if (abun>min_reads) & (abun<max_reads)]

    def identify_abundant_barcodes(self, make_histogram=True, absolute_min_reads=250):

        """

        Identify which barcodes are above the absolute minimal abundance, 

        and make a histogram summarizing the barcode distribution

        """

        keep_barcodes = []

        for k, v in self.barcode_counts.items():

            if v > absolute_min_reads:

                keep_barcodes.append(k)

        abundant_barcodes = {}

        print_to_stderr(" %d barcodes above absolute minimum threshold" % len(keep_barcodes))

        for bc in keep_barcodes:

            abundant_barcodes[bc] = (self.project.stable_barcode_names[self.version][bc], self.barcode_counts[bc])

        self._abundant_barcodes = abundant_barcodes

        with open(self.paths.abundant_barcodes_names_filename, 'w') as f:

            pickle.dump(abundant_barcodes, f)

        # Create table about the filtering process

        with open(self.paths.filtering_statistics_filename, 'w') as filtering_stats:

            header = ['Run', 'Part', 'Input Reads', 'Valid Structure', 'Surviving Trimmomatic', 'Surviving polyA trim and complexity filter']

            if self.version == 'v1' or self.version == 'v2':

                structure_parts = ['W1_in_R2', 'empty_read',  'No_W1', 'No_polyT', 'BC1', 'BC2', 'Umi_error']

                header += ['W1 in R2', 'empty read',  'No W1 in R1', 'No polyT', 'BC1', 'BC2', 'UMI_contains_N']

            elif self.version == 'v3' or self.version == 'v3-miseq':

                structure_parts = ['Invalid_BC1', 'Invalid_BC2', 'UMI_contains_N']

                header += ['Invalid BC1', 'Invalid BC2', 'UMI_contains_N']

            trimmomatic_parts = ['dropped']

            header += ['Dropped by Trimmomatic']

            complexity_filter_parts = ['rejected_because_too_short', 'rejected_because_complexity_too_low']

            header += ['Too short after polyA trim', 'Read complexity too low']

            filtering_stats.write(','.join(header)+'\n')

            for part in self.parts:

                with open(part.filtering_metrics_filename) as f:

                    part_stats = yaml.load(f)

                    line = [part.run_name, part.part_name, part_stats['read_structure']['Total'], part_stats['read_structure']['Valid'], part_stats['trimmomatic']['output'], part_stats['complexity_filter']['output']]

                    line += [part_stats['read_structure'][k] if k in part_stats['read_structure'] else 0 for k in structure_parts]

                    line += [part_stats['trimmomatic'][k] if k in part_stats['trimmomatic'] else 0 for k in trimmomatic_parts]

                    line += [part_stats['complexity_filter'][k] if k in part_stats['complexity_filter'] else 0 for k in complexity_filter_parts]

                    line = [str(l) for l in line]

                    filtering_stats.write(','.join(line)+'\n')

        print_to_stderr("Created Library filtering summary:")

        print_to_stderr("  " + self.paths.filtering_statistics_filename)

        # Make the histogram figure

        if not make_histogram:

            return

        count_freq = defaultdict(int)

        for bc, count in self.barcode_counts.items():

            count_freq[count] += 1

        x = np.array(count_freq.keys())

        y = np.array(count_freq.values())

        w = x*y

        # need to use non-intenactive Agg backend

        import matplotlib

        matplotlib.use('Agg')

        from matplotlib import pyplot as plt

        fig = plt.figure()

        ax = fig.add_subplot(111)

        ax.hist(x, bins=np.logspace(0, 6, 50), weights=w, color='green')

        ax.set_xscale('log')

        ax.set_xlabel('Reads per barcode')

        ax.set_ylabel('#reads coming from bin')

        fig.savefig(self.paths.barcode_abundance_histogram_filename)

        print_to_stderr("Created Barcode Abundance Histogram at:")

        print_to_stderr("  " + self.paths.barcode_abundance_histogram_filename)

        fig = plt.figure()

        ax = fig.add_subplot(111)

        ax.hist(list(self.barcode_counts.values()), bins=np.logspace(2, 6, 50), color='green')

        ax.set_xlim((1, 10**6))

        ax.set_xscale('log')

        ax.set_xlabel('Reads per barcode')

        ax.set_ylabel('# of barcodes')

        fig.savefig(self.paths.barcode_abundance_by_barcode_filename)

        print_to_stderr("Created Barcode Abundance Histogram by barcodes at:")

        print_to_stderr("  " + self.paths.barcode_abundance_by_barcode_filename)

    def sort_reads_by_barcode(self, index=0):

        self.parts[index].sort_reads_by_barcode(self.abundant_barcodes)

    def get_reads_for_barcode(self, barcode, run_filter=[]):

        for part in self.parts:

            if (not run_filter) or (part.run_name in run_filter):

                for line in part.get_reads_for_barcode(barcode):

                    yield line

    def output_barcode_fastq(self, analysis_prefix='', min_reads=750, max_reads=10**10, total_workers=1, worker_index=0, run_filter=[]):

        if analysis_prefix:

            analysis_prefix = analysis_prefix + '.'

        output_dir_path = os.path.join(self.project.project_dir, self.name, 'barcode_fastq')

        self.project.project_check_dir(output_dir_path)

        sorted_barcode_names = self.sorted_barcode_names(min_reads=min_reads, max_reads=max_reads)

        # Identify which barcodes belong to this worker

        barcodes_for_this_worker = []

        i = worker_index

        while i < len(sorted_barcode_names):

            barcodes_for_this_worker.append(sorted_barcode_names[i])

            i += total_workers

        print_to_stderr("""[%s] This worker assigned %d out of %d total barcodes.""" % (self.name, len(barcodes_for_this_worker), len(sorted_barcode_names)))        

        for barcode in barcodes_for_this_worker:

            barcode_fastq_filename = analysis_prefix+'%s.%s.fastq' % (self.name, barcode)

            print_to_stderr("  "+barcode_fastq_filename)

            with open(os.path.join(output_dir_path, barcode_fastq_filename), 'w') as f:

                for line in self.get_reads_for_barcode(barcode, run_filter):

                    f.write(line)

    def quantify_expression(self, analysis_prefix='', max_reads=10**10, min_reads=750, min_counts=0, total_workers=1, worker_index=0, no_bam=False, run_filter=[]):

        if analysis_prefix:

            analysis_prefix = analysis_prefix + '.'

        sorted_barcode_names = self.sorted_barcode_names(min_reads=min_reads, max_reads=max_reads)

        #print_to_stderr("   min_reads: %d sorted_barcode_names counts: %d" % (min_reads, len(sorted_barcode_names)))

        # Identify which barcodes belong to this worker

        barcodes_for_this_worker = []

        i = worker_index

        while i < len(sorted_barcode_names):

            barcodes_for_this_worker.append(sorted_barcode_names[i])

            i += total_workers

        counts_output_filename = os.path.join(self.paths.quant_dir, '%sworker%d_%d.counts.tsv' % (analysis_prefix, worker_index, total_workers))

        ambig_counts_output_filename = os.path.join(self.paths.quant_dir, '%sworker%d_%d.ambig.counts.tsv' % (analysis_prefix, worker_index, total_workers))

        ambig_partners_output_filename = os.path.join(self.paths.quant_dir, '%sworker%d_%d.ambig.partners' % (analysis_prefix, worker_index, total_workers))

        metrics_output_filename = os.path.join(self.paths.quant_dir, '%sworker%d_%d.metrics.tsv' % (analysis_prefix, worker_index, total_workers))

        ignored_for_output_filename = counts_output_filename+'.ignored'

        merged_bam_filename = os.path.join(self.paths.quant_dir, '%sworker%d_%d.bam'% (analysis_prefix, worker_index, total_workers))

        merged_bam_index_filename = merged_bam_filename + '.bai'

        get_barcode_genomic_bam_filename = lambda bc: os.path.join(self.paths.quant_dir, '%s%s.genomic.sorted.bam' % (analysis_prefix, bc))

        # If we wanted BAM output, and the merge BAM and merged BAM index are present, then we are done

        if (not no_bam) and (os.path.isfile(merged_bam_filename) and os.path.isfile(merged_bam_index_filename)):

            print_to_stderr('Indexed, merged BAM file detected for this worker. Done.')

            return 

        # Otherwise, we have to check what we need to quantify

        """

        Function to determine which barcodes this quantification worker might have already quantified.

        This tries to handle interruption during any step of the process.

        The worker is assigned some list of barcodes L. For every barcode:

            - It could have been quantified

                - but have less than min_counts ---> so it got written to `ignored` file.

                - and quantification succeeded, meaning

                    1. there is a line (ending in \n) in the `metrics` file. 

                    2. there is a line (ending in \n) in the `quantification` file.

                    3. there (could) be a line (ending in \n) in the `ambiguous quantification` file.

                    4. there (could) be a line (ending in \n) in the `ambiguous quantification partners` file.

                        [If any line doesn't end in \n, then likely the output of that line was interrupted!]

                    5. (If BAM output is desired) There should be a sorted genomic BAM

                    6. (If BAM output is desired) There should be a sorted genomic BAM index

        """

        succesfully_previously_quantified = set()

        previously_ignored = set()

        header_written = False

        if os.path.isfile(counts_output_filename) and os.path.isfile(metrics_output_filename):

            # Load in list of ignored barcodes

            if os.path.isfile(ignored_for_output_filename):

                with open(ignored_for_output_filename, 'r') as f:

                    previously_ignored = set([line.rstrip().split('\t')[0] for line in f])

            # Load the metrics data into memory

            # (It should be fairly small, this is fast and safe)

            existing_metrics_data = {}

            with open(metrics_output_filename, 'r') as f:

                existing_metrics_data = dict((line.partition('\t')[0], line) for line in f if line[-1]=='\n')

            # Quantification data could be large, read it line by line and output it back for barcodes that have a matching metrics line.

            with open(counts_output_filename, 'r') as in_counts, \

                     open(counts_output_filename+'.tmp', 'w') as tmp_counts, \

                     open(metrics_output_filename+'.tmp', 'w') as tmp_metrics:

                for line in in_counts:

                    # The first worker is reponsible for written the header.

                    # Make sure we carry that over

                    if (not header_written) and (worker_index==0):

                        tmp_counts.write(line)

                        tmp_metrics.write(existing_metrics_data['Barcode'])

                        header_written = True

                        continue

                    # This line has incomplete output, skip it.

                    # (This can only happen with the last line)

                    if line[-1] != '\n':

                        continue

                    barcode = line.partition('\t')[0]

                    # Skip barcode if we don't have existing metrics data

                    if barcode not in existing_metrics_data:

                        continue

                    # Check if we BAM required BAM files exist

                    barcode_genomic_bam_filename = get_barcode_genomic_bam_filename(barcode)

                    bam_files_required_and_present = no_bam or (os.path.isfile(barcode_genomic_bam_filename) and os.path.isfile(barcode_genomic_bam_filename+'.bai'))

                    if not bam_files_required_and_present:

                        continue

                    # This passed all the required checks, write the line to the temporary output files

                    tmp_counts.write(line)

                    tmp_metrics.write(existing_metrics_data[barcode])

                    succesfully_previously_quantified.add(barcode)

            shutil.move(counts_output_filename+'.tmp', counts_output_filename)

            shutil.move(metrics_output_filename+'.tmp', metrics_output_filename)

            # For any 'already quantified' barcode, make sure we also copy over the ambiguity data

            with open(ambig_counts_output_filename, 'r') as in_f, \

                 open(ambig_counts_output_filename+'.tmp', 'w') as tmp_f:

                 f_first_line = (worker_index == 0)

                 for line in in_f:

                    if f_first_line:

                        tmp_f.write(line)

                        f_first_line = False

                        continue

                    if (line.partition('\t')[0] in succesfully_previously_quantified) and (line[-1]=='\n'):

                        tmp_f.write(line)

            shutil.move(ambig_counts_output_filename+'.tmp', ambig_counts_output_filename)

            with open(ambig_partners_output_filename, 'r') as in_f, \

                 open(ambig_partners_output_filename+'.tmp', 'w') as tmp_f:

                 for line in in_f:

                    if (line.partition('\t')[0] in succesfully_previously_quantified) and (line[-1]=='\n'):

                        tmp_f.write(line)

            shutil.move(ambig_partners_output_filename+'.tmp', ambig_partners_output_filename)

        barcodes_to_quantify = [bc for bc in barcodes_for_this_worker if (bc not in succesfully_previously_quantified and bc not in previously_ignored)]

        print_to_stderr("""[%s] This worker assigned %d out of %d total barcodes.""" % (self.name, len(barcodes_for_this_worker), len(sorted_barcode_names)))

        if len(barcodes_for_this_worker)-len(barcodes_to_quantify) > 0:

            print_to_stderr("""    %d previously quantified, %d previously ignored, %d left for this run.""" % (len(succesfully_previously_quantified), len(previously_ignored), len(barcodes_to_quantify)))

        print_to_stderr(('{0:<14.12}'.format('Prefix') if analysis_prefix else '') + '{0:<14.12}{1:<9}'.format("Library", "Barcode"), False)

        print_to_stderr("{0:<8s}{1:<8s}{2:<10s}".format("Reads", "Counts", "Ambigs"))

        for barcode in barcodes_to_quantify:

            self.quantify_expression_for_barcode(barcode,

                counts_output_filename, metrics_output_filename,

                ambig_counts_output_filename, ambig_partners_output_filename,

                no_bam=no_bam, write_header=(not header_written) and (worker_index==0), analysis_prefix=analysis_prefix,

                min_counts = min_counts, run_filter=run_filter)

            header_written = True

        print_to_stderr("Per barcode quantification completed.")

        if no_bam:

            return

        #Gather list of barcodes with output from the metrics file

        genomic_bams = []

        with open(metrics_output_filename, 'r') as f:

            for line in f:

                bc = line.partition('\t')[0]

                if bc == 'Barcode': #This is the line in the header

                    continue

                genomic_bams.append(get_barcode_genomic_bam_filename(bc))

        print_to_stderr("Merging BAM output.")

        try:

            subprocess.check_output([self.project.paths.samtools, 'merge', '-f', merged_bam_filename]+genomic_bams, stderr=subprocess.STDOUT)

        except subprocess.CalledProcessError, err:

            print_to_stderr("   CMD: %s" % str(err.cmd)[:400])

            print_to_stderr("   stdout/stderr:")

            print_to_stderr(err.output)

            raise Exception(" === Error in samtools merge === ")

        print_to_stderr("Indexing merged BAM output.")

        try:

            subprocess.check_output([self.project.paths.samtools, 'index', merged_bam_filename], stderr=subprocess.STDOUT)

        except subprocess.CalledProcessError, err:

            print_to_stderr("   CMD: %s" % str(err.cmd)[:400])

            print_to_stderr("   stdout/stderr:")

            print_to_stderr(err.output)

            raise Exception(" === Error in samtools index === ")

        print(genomic_bams)

        for filename in genomic_bams:

            os.remove(filename)

            os.remove(filename + '.bai')

    def quantify_expression_for_barcode(self, barcode, counts_output_filename, metrics_output_filename,

            ambig_counts_output_filename, ambig_partners_output_filename,

            min_counts=0, analysis_prefix='', no_bam=False, write_header=False, run_filter=[]):

        print_to_stderr(('{0:<14.12}'.format(analysis_prefix) if analysis_prefix else '') + '{0:<14.12}{1:<9}'.format(self.name, barcode), False)

        unaligned_reads_output = os.path.join(self.paths.quant_dir, '%s%s.unaligned.fastq' % (analysis_prefix,barcode))

        aligned_bam = os.path.join(self.paths.quant_dir, '%s%s.aligned.bam' % (analysis_prefix,barcode))

        # Bowtie command

        bowtie_cmd = [self.project.paths.bowtie, self.project.paths.bowtie_index, '-q', '-',

            '-p', '1', '-a', '--best', '--strata', '--chunkmbs', '1000', '--norc', '--sam',

            '-shmem', #should sometimes reduce memory usage...?

            '-m', str(self.project.parameters['bowtie_arguments']['m']),

            '-n', str(self.project.parameters['bowtie_arguments']['n']),

            '-l', str(self.project.parameters['bowtie_arguments']['l']),

            '-e', str(self.project.parameters['bowtie_arguments']['e']),

            ]

        if self.project.parameters['output_arguments']['output_unaligned_reads_to_other_fastq']:

            bowtie_cmd += ['--un', unaligned_reads_output]

        # Quantification command

        script_dir = os.path.dirname(os.path.realpath(__file__))

        quant_cmd = [self.project.paths.python, self.project.paths.quantify_umifm_from_alignments_py,

            '-m', str(self.project.parameters['umi_quantification_arguments']['m']),

            '-u', str(self.project.parameters['umi_quantification_arguments']['u']),

            '-d', str(self.project.parameters['umi_quantification_arguments']['d']),

            '--min_non_polyA', str(self.project.parameters['umi_quantification_arguments']['min_non_polyA']),

            '--library', str(self.name),

            '--barcode', str(barcode),

            '--counts', counts_output_filename,

            '--metrics', metrics_output_filename,

            '--ambigs', ambig_counts_output_filename,

            '--ambig-partners', ambig_partners_output_filename,

            '--min-counts', str(min_counts),

        ]

        if not no_bam:

            quant_cmd += ['--bam', aligned_bam]

        if write_header:

            quant_cmd += ['--write-header']

        if self.project.parameters['umi_quantification_arguments']['split-ambigs']:

            quant_cmd.append('--split-ambig')

        if self.project.parameters['output_arguments']['filter_alignments_to_softmasked_regions']:

            quant_cmd += ['--soft-masked-regions', self.project.paths.bowtie_index + '.soft_masked_regions.pickle']

        # Spawn processes

        p1 = subprocess.Popen(bowtie_cmd, stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE)

        p2 = subprocess.Popen(quant_cmd, stdin=p1.stdout, stderr=subprocess.PIPE)

        for line in self.get_reads_for_barcode(barcode, run_filter=run_filter):

            try:

                p1.stdin.write(line)

            except IOError as e:

                print_to_stderr('\n')

                print_to_stderr(p1.stderr.read())

                raise Exception('\n === Error on piping data to bowtie ===')

        p1.stdin.close()

        if p1.wait() != 0:

            print_to_stderr('\n')

            print_to_stderr(p1.stderr.read())

            raise Exception('\n === Error on bowtie ===')

        if p2.wait() != 0:

            print_to_stderr(p2.stderr.read())

            raise Exception('\n === Error on Quantification Script ===')

        print_to_stderr(p2.stderr.read(), False)

        if no_bam:

            # We are done here

            return False

        if not os.path.isfile(aligned_bam):

            raise Exception("\n === No aligned bam was output for barcode %s ===" % barcode)

        genomic_bam = os.path.join(self.paths.quant_dir, '%s%s.genomic.bam' % (analysis_prefix,barcode))

        sorted_bam = os.path.join(self.paths.quant_dir, '%s%s.genomic.sorted.bam' % (analysis_prefix,barcode))

        try:

            subprocess.check_output([self.project.paths.rsem_tbam2gbam, self.project.paths.bowtie_index, aligned_bam, genomic_bam], stderr=subprocess.STDOUT)

        except subprocess.CalledProcessError, err:

            print_to_stderr("   CMD: %s" % str(err.cmd)[:100])

            print_to_stderr("   stdout/stderr:")

            print_to_stderr(err.output)

            raise Exception(" === Error in rsem-tbam2gbam === ")

        try:

            subprocess.check_output([self.project.paths.samtools, 'sort', '-o', sorted_bam, genomic_bam], stderr=subprocess.STDOUT)

        except subprocess.CalledProcessError, err:

            print_to_stderr("   CMD: %s" % str(err.cmd)[:100])

            print_to_stderr("   stdout/stderr:")

            print_to_stderr(err.output)

            raise Exception(" === Error in samtools sort === ")

        try:

            subprocess.check_output([self.project.paths.samtools, 'index', sorted_bam], stderr=subprocess.STDOUT)

        except subprocess.CalledProcessError, err:

            print_to_stderr("   CMD: %s" % str(err.cmd)[:100])

            print_to_stderr("   stdout/stderr:")

            print_to_stderr(err.output)

            raise Exception(" === Error in samtools index === ")

        os.remove(aligned_bam)

        os.remove(genomic_bam)

        return True

    def aggregate_counts(self, analysis_prefix='', process_ambiguity_data=False):

        if analysis_prefix:

            analysis_prefix = analysis_prefix + '.'

            quant_output_files = [fn[len(analysis_prefix):].split('.')[0] for fn in os.listdir(self.paths.quant_dir) if ('worker' in fn and fn[:len(analysis_prefix)]==analysis_prefix)]

        else:

            quant_output_files = [fn.split('.')[0] for fn in os.listdir(self.paths.quant_dir) if (fn[:6]=='worker')]

        worker_names = [w[6:] for w in quant_output_files]

        worker_indices = set(int(w.split('_')[0]) for w in worker_names)

        total_workers = set(int(w.split('_')[1]) for w in worker_names)

        if len(total_workers) > 1:

            raise Exception("""Quantification for library %s, prefix '%s' was run with different numbers of total_workers.""" % (self.name, analysis_prefix))

        total_workers = list(total_workers)[0]

        missing_workers = []

        for i in range(total_workers):

            if i not in worker_indices:

                missing_workers.append(i)

        if missing_workers:

            missing_workers = ','.join([str(i) for i in sorted(missing_workers)])

            raise Exception("""Output from workers %s (total %d) is missing. """ % (missing_workers, total_workers))

        aggregated_counts_filename = os.path.join(self.project.project_dir, self.name, self.name+'.'+analysis_prefix+'counts.tsv')

        aggregated_quant_metrics_filename = os.path.join(self.project.project_dir, self.name, self.name+'.'+analysis_prefix+'quant_metrics.tsv')

        aggregated_ignored_filename = os.path.join(self.project.project_dir, self.name, self.name+'.'+analysis_prefix+'ignored_barcodes.txt')

        aggregated_bam_output = os.path.join(self.project.project_dir, self.name, self.name+'.'+analysis_prefix+'bam')

        aggregated_ambig_counts_filename = os.path.join(self.project.project_dir, self.name, self.name+'.'+analysis_prefix+'ambig_counts.tsv')

        aggregated_ambig_partners_filename = os.path.join(self.project.project_dir, self.name, self.name+'.'+analysis_prefix+'ambig_partners.tsv')

        agg_counts = open(aggregated_counts_filename, mode='w')

        agg_metrics = open(aggregated_quant_metrics_filename, mode='w')

        agg_ignored = open(aggregated_ignored_filename, mode='w')

        if process_ambiguity_data:

            agg_ambigs = open(aggregated_ambig_counts_filename, mode='w')

            agg_ambig_partners = open(aggregated_ambig_partners_filename, mode='w')

        end_of_counts_header = 0

        end_of_metrics_header = 0

        end_of_ambigs_header = 0

        print_to_stderr('  Concatenating output from all workers.')

        for worker_index in range(total_workers):

            counts_output_filename = os.path.join(self.paths.quant_dir, '%sworker%d_%d.counts.tsv' % (analysis_prefix, worker_index, total_workers))

            ambig_counts_output_filename = os.path.join(self.paths.quant_dir, '%sworker%d_%d.ambig.counts.tsv' % (analysis_prefix, worker_index, total_workers))

            ambig_partners_output_filename = os.path.join(self.paths.quant_dir, '%sworker%d_%d.ambig.partners' % (analysis_prefix, worker_index, total_workers))

            metrics_output_filename = os.path.join(self.paths.quant_dir, '%sworker%d_%d.metrics.tsv' % (analysis_prefix, worker_index, total_workers))

            ignored_for_output_filename = counts_output_filename+'.ignored'

            # Counts

            with open(counts_output_filename, 'r') as f:

                shutil.copyfileobj(f, agg_counts)

            # Metrics

            with open(metrics_output_filename, 'r') as f:

                shutil.copyfileobj(f, agg_metrics)

            # Ignored

            if os.path.isfile(counts_output_filename+'.ignored'):

                with open(counts_output_filename+'.ignored', 'r') as f:

                    shutil.copyfileobj(f, agg_ignored)

            if process_ambiguity_data:

                with open(ambig_counts_output_filename, 'r') as f:

                    shutil.copyfileobj(f, agg_ambigs)

                with open(ambig_partners_output_filename, 'r') as f:

                    shutil.copyfileobj(f, agg_ambig_partners)

        print_to_stderr('  GZIPping concatenated output.')

        agg_counts.close()

        subprocess.Popen(['gzip', '-f', aggregated_counts_filename]).wait()

        agg_metrics.close()

        subprocess.Popen(['gzip', '-f', aggregated_quant_metrics_filename]).wait()

        print_to_stderr('Aggregation completed in %s.gz' % aggregated_counts_filename)

        if process_ambiguity_data:

            agg_ambigs.close()

            subprocess.Popen(['gzip', '-f', aggregated_ambig_counts_filename]).wait()

            agg_ambig_partners.close()

            subprocess.Popen(['gzip', '-f', aggregated_ambig_partners_filename]).wait()

        target_bams = [os.path.join(self.paths.quant_dir, '%sworker%d_%d.bam'% (analysis_prefix, worker_index, total_workers)) for worker_index in range(total_workers)]

        target_bams = [t for t in target_bams if os.path.isfile(t)]

        if target_bams:

            print_to_stderr('  Merging BAM files.')

            p1 = subprocess.Popen([self.project.paths.samtools, 'merge', '-f', aggregated_bam_output]+target_bams, stderr=subprocess.PIPE, stdout=subprocess.PIPE)

            if p1.wait() == 0:

                print_to_stderr('  Indexing merged BAM file.')

                p2 = subprocess.Popen([self.project.paths.samtools, 'index', aggregated_bam_output], stderr=subprocess.PIPE, stdout=subprocess.PIPE)

                if p2.wait() == 0:

                    for filename in target_bams:

                        os.remove(filename)

                        os.remove(filename + '.bai')

                else:

                    print_to_stderr(" === Error in samtools index ===")

                    print_to_stderr(p2.stderr.read())

            else:

                print_to_stderr(" === Error in samtools merge ===")

            print_to_stderr(p1.stderr.read())     

        # print_to_stderr('Deleting per-worker counts files.')

        # for worker_index in range(total_workers):

        #     counts_output_filename = os.path.join(self.paths.quant_dir, '%sworker%d_%d.counts.tsv' % (analysis_prefix, worker_index, total_workers))

        #     os.remove(counts_output_filename)

        #     ambig_counts_output_filename = os.path.join(self.paths.quant_dir, '%sworker%d_%d.ambig.counts.tsv' % (analysis_prefix, worker_index, total_workers))

        #     os.remove(ambig_counts_output_filename)

        #     ambig_partners_output_filename = os.path.join(self.paths.quant_dir, '%sworker%d_%d.ambig.partners' % (analysis_prefix, worker_index, total_workers))

        #     os.remove(ambig_partners_output_filename)

        #     metrics_output_filename = os.path.join(self.paths.quant_dir, '%sworker%d_%d.metrics.tsv' % (analysis_prefix, worker_index, total_workers))

        #     os.remove(metrics_output_filename)

        #     ignored_for_output_filename = counts_output_filename+'.ignored'

        #     os.remove(ignored_for_output_filename)

class LibrarySequencingPart():

    def __init__(self, filtered_fastq_filename=None, project=None, run_name='', library_name='', part_name=''):

        self.project = project

        self.run_name = run_name

        self.part_name = part_name

        self.library_name = library_name

        self.filtered_fastq_filename = filtered_fastq_filename

        self.barcode_counts_pickle_filename = filtered_fastq_filename + '.counts.pickle'

        self.filtering_metrics_filename = '.'.join(filtered_fastq_filename.split('.')[:-1]) + 'metrics.yaml'

        self.sorted_gzipped_fastq_filename = filtered_fastq_filename + '.sorted.fastq.gz'

        self.sorted_gzipped_fastq_index_filename = filtered_fastq_filename + '.sorted.fastq.gz.index.pickle'

    @property

    def is_filtered(self):

        if not hasattr(self, '_is_filtered'):

            self._is_filtered = os.path.exists(self.filtered_fastq_filename) and os.path.exists(self.barcode_counts_pickle_filename)

        return self._is_filtered