--- a
+++ b/indrops.py
@@ -0,0 +1,1786 @@
+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
+    
+    @property
+    def is_sorted(self):
+        if not hasattr(self, '_is_sorted'):
+            self._is_sorted = os.path.exists(self.sorted_gzipped_fastq_filename) and os.path.exists(self.sorted_gzipped_fastq_index_filename)
+        return self._is_sorted
+
+    @property
+    def part_barcode_counts(self):
+        if not hasattr(self, '_part_barcode_counts'):
+            with open(self.barcode_counts_pickle_filename, 'r') as f:
+                self._part_barcode_counts = pickle.load(f)
+        return self._part_barcode_counts
+
+    @property
+    def sorted_index(self):
+        if not hasattr(self, '_sorted_index'):
+            with open(self.sorted_gzipped_fastq_index_filename, 'r') as f:
+                self._sorted_index = pickle.load(f)
+        return self._sorted_index
+
+    def contains_library_in_query(self, query_libraries):
+        return self.library_name in query_libraries
+
+    def sort_reads_by_barcode(self, abundant_barcodes={}):
+        sorted_barcodes = [j for j,v in sorted(abundant_barcodes.items(), key=lambda i:-i[1][1])]
+        sorted_barcodes = [j for j in sorted_barcodes if j in self.part_barcode_counts]
+
+        barcode_buffers = {}
+        barcode_gzippers = {}
+        for bc in sorted_barcodes + ['ignored']:
+            barcode_buffers[bc] = BytesIO()
+            barcode_gzippers[bc] = gzip.GzipFile(fileobj=barcode_buffers[bc], mode='wb')
+
+        total_processed_reads = 0
+        total_ignored_reads = 0
+        bcs_with_data = set()
+        bcs_with_tmp_data = set()
+        barcode_tmp_filename = lambda bc: '%s.%s.tmp.gz' % (self.sorted_gzipped_fastq_filename, bc)
+
+
+        total_reads = sum(self.part_barcode_counts.values())
+        print_to_stderr('Sorting %d reads from %d barcodes above absolute minimum threshold.' % (total_reads, len(abundant_barcodes)))
+        with open(self.filtered_fastq_filename, 'r') as input_fastq:
+            for name, seq, qual in from_fastq(input_fastq):
+                total_processed_reads += 1
+                bc = name.split(':')[0]
+
+                if total_processed_reads%1000000 == 0:
+                    print_to_stderr('Read in %.02f percent of all reads (%d)' % (100.*total_processed_reads/total_reads, total_processed_reads))
+                
+                if bc in abundant_barcodes:
+                    barcode_gzippers[bc].write(to_fastq(name, seq, qual))
+                    bcs_with_data.add(bc)
+                else:
+                    total_ignored_reads += 1
+                    barcode_gzippers['ignored'].write(to_fastq(name, seq, qual))
+                    bcs_with_data.add('ignored')
+
+
+        sorted_output_index = {}
+        with open(self.sorted_gzipped_fastq_filename, 'wb') as sorted_output:
+            for original_bc in sorted_barcodes + ['ignored']:
+                if original_bc != 'ignored':
+                    new_bc_name = abundant_barcodes[original_bc][0]
+                    barcode_reads_count = self.part_barcode_counts[original_bc]
+                else:
+                    new_bc_name = 'ignored'
+                    barcode_reads_count = total_ignored_reads
+
+                start_pos = sorted_output.tell()
+                barcode_gzippers[original_bc].close()
+                if original_bc in bcs_with_data:
+                    barcode_buffers[original_bc].seek(0)
+                    shutil.copyfileobj(barcode_buffers[original_bc], sorted_output)
+                barcode_buffers[original_bc].close()
+                end_pos = sorted_output.tell()
+
+                if end_pos > start_pos:
+                    sorted_output_index[new_bc_name] = (original_bc, start_pos, end_pos, end_pos-start_pos, barcode_reads_count)
+
+        with open(self.sorted_gzipped_fastq_index_filename, 'w') as f:
+            pickle.dump(sorted_output_index, f)      
+
+    def get_reads_for_barcode(self, barcode):
+        if barcode not in self.sorted_index:
+            raise StopIteration
+
+        original_barcode, start_byte_offset, end_byte_offset, byte_length, barcode_reads = self.sorted_index[barcode]
+
+        with open(self.sorted_gzipped_fastq_filename, 'rb') as sorted_output:
+            sorted_output.seek(start_byte_offset)
+            byte_buffer = BytesIO(sorted_output.read(byte_length))
+            ungzipper = gzip.GzipFile(fileobj=byte_buffer, mode='rb')
+            while True:
+                yield next(ungzipper)
+
+    @contextmanager
+    def trimmomatic_and_low_complexity_filter_process(self):
+        """
+        We start 3 processes that are connected with Unix pipes.
+
+        Process 1 - Trimmomatic. Doesn't support stdin/stdout, so we instead use named pipes (FIFOs). It reads from FIFO1, and writes to FIFO2. 
+        Process 2 - In line complexity filter, a python script. It reads from FIFO2 (Trimmomatic output) and writes to the ouput file. 
+        Process 3 - Indexer that counts the number of reads for every barcode. This reads from stdin, writes the reads to stdout and writes the index as a pickle to stderr.
+
+        When these are done, we start another process to count the results on the FastQ file.
+        """
+        filtered_dir = os.path.dirname(self.filtered_fastq_filename) #We will use the same directory for creating temporary FIFOs, assuming we have write access.
+        
+        self.filtering_statistics_counter = defaultdict(int)
+        with FIFO(dir=filtered_dir) as fifo2, open(self.filtered_fastq_filename, 'w') as filtered_fastq_file, open(self.filtered_fastq_filename+'.counts.pickle', 'w') as filtered_index_file:
+            
+            low_complexity_filter_cmd = [self.project.paths.python, self.project.paths.trim_polyA_and_filter_low_complexity_reads_py,
+                '-input', fifo2.filename, 
+                '--min-post-trim-length', self.project.parameters['trimmomatic_arguments']['MINLEN'],
+                '--max-low-complexity-fraction', str(self.project.parameters['low_complexity_filter_arguments']['max_low_complexity_fraction']),
+                ]
+            counter_cmd = [self.project.paths.python,  self.project.paths.count_barcode_distribution_py]
+
+            p2 = subprocess.Popen(low_complexity_filter_cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
+            p3 = subprocess.Popen(counter_cmd, stdin=p2.stdout, stdout=filtered_fastq_file, stderr=filtered_index_file)
+
+            with FIFO(dir=filtered_dir) as fifo1:
+
+                trimmomatic_cmd = [self.project.paths.java, '-Xmx500m', '-jar', self.project.paths.trimmomatic_jar,
+                        'SE', '-threads', "1", '-phred33', fifo1.filename, fifo2.filename]
+                for arg in self.project.parameters['trimmomatic_arguments']['argument_order']:
+                    val = self.project.parameters['trimmomatic_arguments'][arg]
+                    trimmomatic_cmd.append('%s:%s' % (arg, val))
+
+                p1 = subprocess.Popen(trimmomatic_cmd, stderr=subprocess.PIPE)
+
+                fifo1_filehandle = open(fifo1.filename, 'w')
+                yield fifo1_filehandle
+                fifo1_filehandle.close()
+                trimmomatic_stderr = p1.stderr.read().splitlines()
+                if trimmomatic_stderr[2] != 'TrimmomaticSE: Completed successfully':
+                    raise Exception('Trimmomatic did not complete succesfully on %s' % filtered_filename)
+                trimmomatic_metrics = trimmomatic_stderr[1].split() 
+                # ['Input', 'Reads:', #READS, 'Surviving:', #SURVIVING, (%SURVIVING), 'Dropped:', #DROPPED, (%DROPPED)]
+                trimmomatic_metrics = {'input' : trimmomatic_metrics[2], 'output': trimmomatic_metrics[4], 'dropped': trimmomatic_metrics[7]}
+                p1.wait()
+
+            complexity_filter_metrics = pickle.load(p2.stderr)
+            p2.wait()
+            p3.wait()
+
+
+        filtering_metrics = {
+            'read_structure' : dict(self.filtering_statistics_counter),
+            'trimmomatic' : trimmomatic_metrics,
+            'complexity_filter': complexity_filter_metrics,
+        }
+        with open(self.filtering_metrics_filename, 'w') as f:
+            yaml.dump(dict(filtering_metrics), f, default_flow_style=False)
+
+
+class V1V2Filtering(LibrarySequencingPart):
+
+    def __init__(self, bioread_filename=None, metaread_filename=None, *args, **kwargs):
+
+        self.bioread_filename = bioread_filename
+        self.metaread_filename = metaread_filename
+        LibrarySequencingPart.__init__(self, *args, **kwargs)
+
+
+    def filter_and_count_reads(self):
+        """
+        Input the two raw FastQ files
+        Output: 
+            - A single fastQ file that uses the read name to store the barcoding information
+            - A pickle of the number of reads originating from each barcode 
+        """
+        # Relevant paths
+        r1_filename, r2_filename = self.metaread_filename, self.bioread_filename
+
+        #Get barcode neighborhoods
+        bc1s = self.project.gel_barcode1_revcomp_list_neighborhood
+        bc2s = self.project.gel_barcode2_revcomp_list_neighborhood 
+
+
+        # This starts a Trimmomatic process, a low complexity filter process, and will 
+        # upon closing, start the barcode distribution counting process.
+        last_ping = time.time()
+        ping_every_n_reads = 1000000
+        ping_header = "{0:>12}{1:>16}{2:>12}{3:>10}{4:>10}{5:>10}{6:>10}{7:>10}{8:>10}{9:>10}"
+        ping_header = ping_header.format("Total Reads", "", "Valid Reads", "W1 in R2", "Empty", "No W1", "No polyT", "No BC1", "No BC2", "No UMI")
+        ping_template = "{total:12d}    {rate:5.1f} sec/M {Valid:12.1%}{W1_in_R2:10.1%}{empty_read:10.1%}{No_W1:10.1%}{No_polyT:10.1%}{BC1:10.1%}{BC2:10.1%}{Umi_error:10.1%}"
+        def print_ping_to_log(last_ping):
+            sec_per_mil = (time.time()-last_ping)/(ping_every_n_reads/10**6) if last_ping else 0.0
+            total = self.filtering_statistics_counter['Total']
+            if total > 0:
+                ping_format_data = {k: float(self.filtering_statistics_counter[k])/total for k in ['Valid', 'W1_in_R2', 'empty_read',  'No_W1', 'No_polyT', 'BC1', 'BC2', 'Umi_error']}
+                print_to_stderr(ping_template.format(total=total, rate=sec_per_mil, **ping_format_data))
+
+
+        with self.trimmomatic_and_low_complexity_filter_process() as trim_process:
+            #Iterate over the weaved reads
+            for r_name, r1_seq, r1_qual, r2_seq, r2_qual in self._weave_fastqs(r1_filename, r2_filename):
+                    
+                # Check if they should be kept
+                keep, result = self._process_reads(r1_seq, r2_seq, valid_bc1s=bc1s, valid_bc2s=bc2s)
+
+                # Write the the reads worth keeping
+                if keep:
+                    bc, umi = result
+                    trim_process.write(to_fastq_lines(bc, umi, r2_seq, r2_qual, r_name))
+                    self.filtering_statistics_counter['Valid'] += 1
+                else:
+                    self.filtering_statistics_counter[result] += 1
+
+                # Track speed per M reads
+                self.filtering_statistics_counter['Total'] += 1
+                if self.filtering_statistics_counter['Total']%(10*ping_every_n_reads) == 1:
+                    print_to_stderr(ping_header)
+
+                if self.filtering_statistics_counter['Total']%ping_every_n_reads == 0:
+                    print_ping_to_log(last_ping)
+                    last_ping = time.time()
+
+            print_ping_to_log(False)
+
+        print_to_stderr(self.filtering_statistics_counter)
+
+    def _weave_fastqs(self, r1_fastq, r2_fastq):
+        """
+        Merge 2 FastQ files by returning paired reads for each.
+        Returns only R1_seq, R2_seq and R2_qual.
+        """
+
+        is_gz_compressed = False
+        is_bz_compressed = False
+        if r1_fastq.split('.')[-1] == 'gz' and r2_fastq.split('.')[-1] == 'gz':
+            is_gz_compressed = True
+            
+        #Added bz2 support VS
+        if r1_fastq.split('.')[-1] == 'bz2' and r2_fastq.split('.')[-1] == 'bz2':
+            is_bz_compressed = True
+
+        # Decompress Gzips using subprocesses because python gzip is incredibly slow.
+        if is_gz_compressed:    
+            r1_gunzip = subprocess.Popen("gzip --stdout -d %s" % (r1_fastq), shell=True, stdout=subprocess.PIPE)
+            r1_stream = r1_gunzip.stdout
+            r2_gunzip = subprocess.Popen("gzip --stdout -d %s" % (r2_fastq), shell=True, stdout=subprocess.PIPE)
+            r2_stream = r2_gunzip.stdout
+        elif is_bz_compressed:
+            r1_bunzip = subprocess.Popen("bzcat %s" % (r1_fastq), shell=True, stdout=subprocess.PIPE)
+            r1_stream = r1_bunzip.stdout
+            r2_bunzip = subprocess.Popen("bzcat %s" % (r2_fastq), shell=True, stdout=subprocess.PIPE)
+            r2_stream = r2_bunzip.stdout
+        else:
+            r1_stream = open(r1_fastq, 'r')
+            r2_stream = open(r2_fastq, 'r')
+
+        while True:
+            #Read 4 lines from each FastQ
+            name = next(r1_stream).rstrip()[1:].split()[0] #Read name
+            r1_seq = next(r1_stream).rstrip() #Read seq
+            next(r1_stream) #+ line
+            r1_qual = next(r1_stream).rstrip() #Read qual
+            
+            next(r2_stream) #Read name
+            r2_seq = next(r2_stream).rstrip() #Read seq
+            next(r2_stream) #+ line
+            r2_qual = next(r2_stream).rstrip() #Read qual
+            
+            # changed to allow for empty reads (caused by adapter trimming)
+            if name:
+                yield name, r1_seq, r1_qual, r2_seq, r2_qual
+            else:
+            # if not r1_seq or not r2_seq:
+                break
+
+        r1_stream.close()
+        r2_stream.close()
+
+    def _process_reads(self, name, read, valid_bc1s={}, valid_bc2s={}):
+        """
+        Returns either:
+            True, (barcode, umi)
+                (if read passes filter)
+            False, name of filter that failed
+                (for stats collection)
+        
+        R1 anatomy: BBBBBBBB[BBB]WWWWWWWWWWWWWWWWWWWWWWCCCCCCCCUUUUUUTTTTTTTTTT______________
+            B = Barcode1, can be 8, 9, 10 or 11 bases long.
+            W = 'W1' sequence, specified below
+            C = Barcode2, always 8 bases
+            U = UMI, always 6 bases
+            T = Beginning of polyT tail.
+            _ = Either sequencing survives across the polyT tail, or signal starts dropping off
+                (and start being anything, likely with poor quality)
+        """
+
+        minimal_polyT_len_on_R1 = 7
+        hamming_threshold_for_W1_matching = 3
+
+        w1 = "GAGTGATTGCTTGTGACGCCTT"
+        rev_w1 = "AAGGCGTCACAAGCAATCACTC" #Hard-code so we don't recompute on every one of millions of calls
+        # If R2 contains rev_W1, this is almost certainly empty library
+        if rev_w1 in read:
+            return False, 'W1_in_R2'
+
+        # # With reads sufficiently long, we will often see a PolyA sequence in R2. 
+        # if polyA in read:
+        #     return False, 'PolyA_in_R2'
+
+        # Check for polyT signal at 3' end.
+        # 44 is the length of BC1+W1+BC2+UMI, given the longest PolyT
+        #BC1: 8-11 bases
+        #W1 : 22 bases
+        #BC2: 8 bases
+        #UMI: 6 bases
+
+        # check for empty reads (due to adapter trimming)
+        if not read:
+            return False, 'empty_read'
+        
+        #Check for W1 adapter
+        #Allow for up to hamming_threshold errors
+        if w1 in name:
+            w1_pos = name.find(w1)
+            if not 7 < w1_pos < 12:
+                return False, 'No_W1'
+        else:
+            #Try to find W1 adapter at start positions 8-11
+            #by checking hamming distance to W1.
+            for w1_pos in range(8, 12):
+                if string_hamming_distance(w1, name[w1_pos:w1_pos+22]) <= hamming_threshold_for_W1_matching:
+                    break
+            else:
+                return False, 'No_W1'
+                
+        bc2_pos=w1_pos+22
+        umi_pos=bc2_pos+8
+        polyTpos=umi_pos+6
+        expected_poly_t = name[polyTpos:polyTpos+minimal_polyT_len_on_R1]
+        if string_hamming_distance(expected_poly_t, 'T'*minimal_polyT_len_on_R1) > 3:
+                 return False, 'No_polyT'
+            
+        bc1 = str(name[:w1_pos])
+        bc2 = str(name[bc2_pos:umi_pos])
+        umi = str(name[umi_pos:umi_pos+6])
+        
+        #Validate barcode (and try to correct when there is no ambiguity)
+        if valid_bc1s and valid_bc2s:
+            # Check if BC1 and BC2 can be mapped to expected barcodes
+            if bc1 in valid_bc1s:
+                # BC1 might be a neighboring BC, rather than a valid BC itself. 
+                bc1 = valid_bc1s[bc1]
+            else:
+                return False, 'BC1'
+            if bc2 in valid_bc2s:
+                bc2 = valid_bc2s[bc2]
+            else:
+                return False, 'BC2'
+            if 'N' in umi:
+                return False, 'UMI_error'
+        bc = '%s-%s'%(bc1, bc2)
+        return True, (bc, umi)
+
+class V3Demultiplexer():
+
+    def __init__(self, library_indices, project=None, part_filename="", input_filename="", run_name="", part_name="", run_version_details="v3"):
+
+        self.run_version_details = run_version_details
+        self.input_filename = input_filename
+        self.project = project
+        self.run_name = run_name
+        self.part_name = part_name
+        self.libraries = {}
+        for lib in library_indices:
+            lib_index = lib['library_index']
+            lib_name = lib['library_name']
+            library_part_filename = part_filename.format(library_name=lib_name, library_index=lib_index)
+            self.libraries[lib_index] = LibrarySequencingPart(filtered_fastq_filename=library_part_filename, project=project, run_name=run_name, library_name=lib_name, part_name=part_name)
+
+    def _weave_fastqs(self, fastqs):
+        last_extension = [fn.split('.')[-1] for fn in fastqs]
+        if all(ext == 'gz' for ext in last_extension):
+            processes = [subprocess.Popen("gzip --stdout -d %s" % (fn), shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) for fn in fastqs]
+            streams = [r.stdout for r in processes]
+        elif all(ext == 'bz2' for ext in last_extension):
+            processes = [subprocess.Popen("bzcat %s" % (fn), shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) for fn in fastqs]
+            streams = [r.stdout for r in processes]
+        elif all(ext == 'fastq' for ext in last_extension):
+            streams = [open(fn, 'r') for fn in fastqs]
+        else:
+            raise("ERROR: Different files are compressed differently. Check input.")
+
+        while True:
+            names = [next(s)[:-1].split()[0] for s in streams]
+            seqs = [next(s)[:-1] for s in streams]
+            blanks = [next(s)[:-1]  for s in streams]
+            quals = [next(s)[:-1]  for s in streams]
+
+            assert all(name==names[0] for name in names)
+            yield names[0], seqs, quals
+
+        for s in streams:
+            s.close()
+
+
+    def _process_reads(self, name, seqs, quals, valid_bc1s={}, valid_bc2s={}, valid_libs={}):
+        """
+        Returns either:
+            True, (barcode, umi)
+                (if read passes filter)
+            False, name of filter that failed
+                (for stats collection)
+        """
+
+        r1, r2, r3, r4 = seqs
+        if self.run_version_details=='v3-miseq':
+            r2 = rev_comp(r2)
+            r4 = rev_comp(r4)
+
+        if r3 in valid_libs:
+            lib_index = valid_libs[r3]
+        else:
+            return False, r3, 'Invalid_library_index'
+
+        if r2 in valid_bc1s:
+            bc1 = valid_bc1s[r2]
+        else:
+            return False, lib_index, 'Invalid_BC1'
+
+        orig_bc2 = r4[:8]
+        umi = r4[8:8+6]
+        polyA = r4[8+6:]
+
+        if orig_bc2 in valid_bc2s:
+            bc2 = valid_bc2s[orig_bc2]
+        else:
+            return False, lib_index, 'Invalid_BC2'
+
+        if 'N' in umi:
+            return False, lib_index, 'UMI_contains_N'
+
+        final_bc = '%s-%s' % (bc1, bc2)
+        return True, lib_index, (final_bc, umi)
+
+
+    def filter_and_count_reads(self):
+        # Prepare error corrected index sets
+        self.sequence_to_index_mapping = {}
+        libs = self.libraries.keys()
+        self.sequence_to_index_mapping = dict(zip(libs, libs))
+        index_neighborhoods = [set(seq_neighborhood(lib, 1)) for lib in libs]
+        for lib, clibs in zip(libs, index_neighborhoods):
+            # Quick check that error-correction maps to a single index
+            for clib in clibs:
+                if sum(clib in hood for hood in index_neighborhoods)==1:
+                    self.sequence_to_index_mapping[clib] = lib
+
+        # Prepare error corrected barcode sets
+        error_corrected_barcodes = self.project.gel_barcode2_list_neighborhood
+        error_corrected_rev_compl_barcodes = self.project.gel_barcode2_revcomp_list_neighborhood
+
+        # Open up our context managers
+        manager_order = [] #It's imperative to exit managers the opposite order than we open them!
+        trim_processes = {}
+        trim_processes_managers = {}
+
+        for lib in self.libraries.keys():
+            manager_order.append(lib)
+            trim_processes_managers[lib] = self.libraries[lib].trimmomatic_and_low_complexity_filter_process()
+            trim_processes[lib] = trim_processes_managers[lib].__enter__()
+
+        overall_filtering_statistics = defaultdict(int)
+
+        # Paths for the 4 expected FastQs
+        input_fastqs = []
+        for r in ['R1', 'R2', 'R3', 'R4']:
+            input_fastqs.append(self.input_filename.format(read=r))
+
+        last_ping = time.time()
+        ping_every_n_reads = 1000000
+        ping_header = "{0:>12}{1:>16}{2:>12}{3:>10}{4:>10}{5:>10}{6:>10}   |" + ''.join("{%d:>12.10}"%i for i in range(7,7+len(manager_order)))
+        ping_header = ping_header.format("Total Reads", "", "Valid Reads", "No index", "No BC1", "No BC2", "No UMI", *[self.libraries[k].library_name for k in manager_order])
+        ping_template = "{total:12d}    {rate:5.1f} sec/M {Valid:12.1%}{Invalid_library_index:10.1%}{Invalid_BC1:10.1%}{Invalid_BC2:10.1%}{UMI_contains_N:10.1%}   |{"+":>12.1%}{".join(manager_order)+":>12.1%}"
+        
+        def print_ping_to_log(last_ping):
+            sec_per_mil = (time.time() - last_ping)/(float(ping_every_n_reads)/10**6) if last_ping else 0
+            total = overall_filtering_statistics['Total']
+            ping_format_data = {k: float(overall_filtering_statistics[k])/total for k in ['Valid', 'Invalid_library_index', 'Invalid_BC1',  'Invalid_BC2', 'UMI_contains_N']}
+            if overall_filtering_statistics['Valid'] > 0:
+                ping_format_data.update({k: float(self.libraries[k].filtering_statistics_counter['Valid'])/overall_filtering_statistics['Valid'] for k in manager_order})
+            print_to_stderr(ping_template.format(total=total, rate=sec_per_mil, **ping_format_data))
+
+        common__ = defaultdict(int)
+        print_to_stderr('Filtering %s, file %s' % (self.run_name, self.input_filename))
+        for r_name, seqs, quals in self._weave_fastqs(input_fastqs):
+
+            # Python 3 compatibility in mind!
+            seqs = [s.decode('utf-8') for s in seqs]
+
+            keep, lib_index, result = self._process_reads(r_name, seqs, quals,
+                                                    error_corrected_barcodes, error_corrected_rev_compl_barcodes, 
+                                                    self.sequence_to_index_mapping)
+            common__[seqs[1]] += 1
+            if keep:
+                bc, umi = result
+                bio_read = seqs[0]
+                bio_qual = quals[0]
+                trim_processes[lib_index].write(to_fastq_lines(bc, umi, bio_read, bio_qual, r_name[1:]))
+                self.libraries[lib_index].filtering_statistics_counter['Valid'] += 1
+                self.libraries[lib_index].filtering_statistics_counter['Total'] += 1
+                overall_filtering_statistics['Valid'] += 1
+
+            else:
+                if result != 'Invalid_library_index':
+                    self.libraries[lib_index].filtering_statistics_counter[result] += 1
+                    self.libraries[lib_index].filtering_statistics_counter['Total'] += 1
+                overall_filtering_statistics[result] += 1
+
+            # Track speed per M reads
+            overall_filtering_statistics['Total'] += 1
+
+            if overall_filtering_statistics['Total']%(ping_every_n_reads*10)==1:
+                print_to_stderr(ping_header)
+            
+            if overall_filtering_statistics['Total']%ping_every_n_reads == 0:
+                print_ping_to_log(last_ping)
+                last_ping = time.time()
+                
+        print_ping_to_log(False)
+        # Close up the context managers
+        for lib in manager_order[::-1]:
+            trim_processes_managers[lib].__exit__(None, None, None)
+
+    def contains_library_in_query(self, query_libraries):
+        for lib in self.libraries.values():
+            if lib.contains_library_in_query(query_libraries):
+                return True
+        return False
+
+
+
+
+
+if __name__=="__main__":
+
+    import sys, argparse
+    parser = argparse.ArgumentParser()
+
+    parser.add_argument('project', type=argparse.FileType('r'), help='Project YAML File.')
+    parser.add_argument('-l', '--libraries', type=str, help='[all] Library name(s) to work on. If blank, will iterate over all libraries in project.', nargs='?', default='')
+    parser.add_argument('-r', '--runs', type=str, help='[all] Run name(s) to work on. If blank, will iterate over all runs in project.', nargs='?', default='')
+    parser.add_argument('command', type=str, choices=['info', 'filter', 'identify_abundant_barcodes', 'sort', 'quantify', 'aggregate', 'build_index', 'get_reads', 'output_barcode_fastq'])
+    parser.add_argument('--total-workers', type=int, help='[all] Total workers that are working together. This takes precedence over barcodes-per-worker.', default=1)
+    parser.add_argument('--worker-index', type=int, help='[all] Index of current worker (the first worker should have index 0).', default=0)
+    parser.add_argument('--min-reads', type=int, help='[quantify] Minimum number of reads for barcode to be processed', nargs='?', default=750)
+    parser.add_argument('--max-reads', type=int, help='[quantify] Maximum number of reads for barcode to be processed', nargs='?', default=100000000)
+    parser.add_argument('--min-counts', type=int, help='[aggregate] Minimun number of UMIFM counts for barcode to be aggregated', nargs='?', default=0)
+    parser.add_argument('--analysis-prefix', type=str, help='[quantify/aggregate/convert_bam/merge_bam] Prefix for analysis files.', nargs='?', default='')
+    parser.add_argument('--no-bam', help='[quantify] Do not output alignments to bam file.', action='store_true')
+    parser.add_argument('--genome-fasta-gz', help='[build_index] Path to gzipped soft-masked genomic FASTA file.')
+    parser.add_argument('--ensembl-gtf-gz', help='[build_index] Path to gzipped ENSEMBL GTF file. ')
+    parser.add_argument('--mode', help='[build_index] Stringency mode for transcriptome build. [strict|all_ensembl]', default='strict')
+    parser.add_argument('--override-yaml', help="[all] Dictionnary to update project YAML with.. [You don't need this.]", nargs='?', default='')
+
+    args = parser.parse_args()
+    project = IndropsProject(args.project)
+    if args.override_yaml:
+        override = eval(args.override_yaml)
+        if 'paths' in override:
+            project.yaml['paths'].update(override['paths'])
+        if 'parameters' in override:
+            for k,v in override['parameters'].items():
+                project.yaml['parameters'][k].update(v)
+        if hasattr(project, '_paths'):
+            del project._paths
+        if hasattr(project, '_parameters'):
+            del project._parameters
+
+    target_libraries = []
+    if args.libraries:
+        for lib in args.libraries.split(','):
+            assert lib in project.libraries
+            if lib not in target_libraries:
+                target_libraries.append(lib)
+    else:
+        target_libraries = project.libraries.keys()
+    lib_query = set(target_libraries)
+
+    target_runs = []
+    if args.runs:
+        for run in args.runs.split(','):
+            assert run in project.runs
+            target_runs.append(run)
+    else:
+        target_runs = project.runs.keys()
+
+    target_library_parts = []
+    for lib in target_libraries:
+        for pi, part in enumerate(project.libraries[lib].parts):
+            if part.run_name in target_runs:
+                target_library_parts.append((lib, pi))
+
+    if args.command == 'info':
+        print_to_stderr('Project Name: ' + project.name)
+        target_run_parts = []
+        for run in target_runs:
+            target_run_parts += [part for part in project.runs[run] if part.contains_library_in_query(lib_query)]
+        print_to_stderr('Total library parts in search query: ' + str(len(target_run_parts)))
+
+    elif args.command == 'filter':
+        target_run_parts = []
+        for run in target_runs:
+            target_run_parts += [part for part in project.runs[run] if part.contains_library_in_query(lib_query)]
+
+        for part in worker_filter(target_run_parts, args.worker_index, args.total_workers):
+            print_to_stderr('Filtering run "%s", library "%s", part "%s"' % (part.run_name, part.library_name if hasattr(part, 'library_name') else 'N/A', part.part_name))
+            part.filter_and_count_reads()
+
+    elif args.command == 'identify_abundant_barcodes':
+        for library in worker_filter(target_libraries, args.worker_index, args.total_workers):
+            project.libraries[library].identify_abundant_barcodes()
+
+    elif args.command == 'sort':
+        for library, part_index in worker_filter(target_library_parts, args.worker_index, args.total_workers):
+            print_to_stderr('Sorting %s, part "%s"' % (library, project.libraries[library].parts[part_index].filtered_fastq_filename))
+            project.libraries[library].sort_reads_by_barcode(index=part_index)
+
+    elif args.command == 'quantify':
+        for library in target_libraries:
+            project.libraries[library].quantify_expression(worker_index=args.worker_index, total_workers=args.total_workers,
+                    min_reads=args.min_reads, max_reads=args.max_reads, min_counts=args.min_counts,
+                    analysis_prefix=args.analysis_prefix,
+                    no_bam=args.no_bam, run_filter=target_runs)
+
+            for part in project.libraries[library].parts:
+                if hasattr(part, '_sorted_index'):
+                    del part._sorted_index
+
+    elif args.command == 'aggregate':
+        for library in target_libraries:
+            project.libraries[library].aggregate_counts(analysis_prefix=args.analysis_prefix)
+
+    elif args.command == 'build_index':
+        project.build_transcriptome(args.genome_fasta_gz, args.ensembl_gtf_gz, mode=args.mode)
+
+    elif args.command == 'get_reads':
+        for library in target_libraries:
+            sorted_barcode_names = project.libraries[library].sorted_barcode_names(min_reads=args.min_reads, max_reads=args.max_reads)
+            for bc in sorted_barcode_names:
+                for line in project.libraries[library].get_reads_for_barcode(bc, run_filter=target_runs):
+                    sys.stdout.write(line)
+
+            for part in project.libraries[library].parts:
+                if hasattr(part, '_sorted_index'):
+                    del part._sorted_index
+
+    elif args.command == 'output_barcode_fastq':
+        for library in target_libraries:
+            project.libraries[library].output_barcode_fastq(worker_index=args.worker_index, total_workers=args.total_workers,
+                    min_reads=args.min_reads, max_reads=args.max_reads, analysis_prefix=args.analysis_prefix, run_filter=target_runs)
\ No newline at end of file