class FileEntry(MyEntry):
    """
    Creates a labeled Entry field for a file or directory.

    *text* is the Label/error box text.
    *directory* is ``True`` if selecting a directory (instead of a file).
    *extensions* is a list of valid file endings

    """

    def __init__(
        self, text, cfg, key, optional=False, directory=False, extensions=None
    ):
        MyEntry.__init__(self, text, cfg, key, optional)
        self.choose = None
        self.clear = None

        self.directory = directory
        if extensions is not None:
            self.extensions = [x.lower() for x in extensions]
        else:
            self.extensions = None

    def body(self, master, row, columns=DEFAULT_COLUMNS, **kwargs):
        """
        Place the required elements using the grid layout method.

        Returns the number of rows taken by this element.
        """
        label = ttk.Label(master, text=self.text)
        label.grid(row=row, column=0, columnspan=1, sticky="e")
        self.entry = ttk.Entry(master, textvariable=self.value)
        self.entry.grid(row=row, column=1, columnspan=columns - 1, sticky="ew")
        if self.directory:
            self.choose = ttk.Button(
                master,
                text="Choose...",
                command=lambda: self.value.set(tkFileDialog.askdirectory()),
            )
        else:
            self.choose = ttk.Button(
                master,
                text="Choose...",
                command=lambda: self.value.set(tkFileDialog.askopenfilename()),
            )
        self.choose.grid(row=row + 1, column=1, sticky="w")
        if self.optional:
            self.clear = ttk.Button(
                master, text="Clear", command=lambda: self.value.set("")
            )
            self.clear.grid(row=row + 1, column=2, sticky="e")
        return 2

    def validate(self):
        if not self.enabled:
            return True
        elif len(self.value.get()) == 0:
            if not self.optional:
                tkMessageBox.showwarning("", "{} not specified.".format(self.text))
                return False
            else:
                return True
        else:
            if os.path.exists(self.value.get()):
                if self.extensions is not None:
                    if any(
                        self.value.get().lower().endswith(x) for x in self.extensions
                    ):
                        return True
                    else:
                        tkMessageBox.showwarning(
                            "", "Invalid file extension " "for {}.".format(self.text)
                        )
                        return False
                else:  # no extension restriction
                    return True
            else:
                tkMessageBox.showwarning(
                    "", "{} file does not exist." "".format(self.text)
                )
                return False

    def enable(self):
        self.enabled = True
        self.entry.state(["!disabled"])
        self.choose.state(["!disabled"])
        if self.optional:
            self.clear.state(["!disabled"])

    def disable(self):
        self.enabled = False
        self.entry.state(["disabled"])
        self.choose.state(["disabled"])
        if self.optional:
            self.clear.state(["disabled"]) class StringEntry(MyEntry):
    """
    Creates a labeled Entry field for a string.

    *text* is the Label/error box text.
    """

    def __init__(self, text, cfg, key, optional=False):
        MyEntry.__init__(self, text, cfg, key, optional)

    def body(self, master, row, columns=DEFAULT_COLUMNS, **kwargs):
        """
        Place the required elements using the grid layout method.

        Returns the number of rows taken by this element.
        """
        label = ttk.Label(master, text=self.text)
        label.grid(row=row, column=0, columnspan=1, sticky="e")
        self.entry = ttk.Entry(master, textvariable=self.value)
        self.entry.grid(row=row, column=1, columnspan=columns - 1, sticky="ew")
        return 1 DEFAULT_COLUMNS = 3 SEQLIB_LABEL_TEXT = OrderedDict(
    [
        ("BcvSeqLib", "Barcoded Variant"),
        ("BcidSeqLib", "Barcoded Identifier"),
        ("OverlapSeqLib", "Overlap"),
        ("BasicSeqLib", "Basic"),
        ("BarcodeSeqLib", "Barcodes Only"),
        ("IdOnlySeqLib", "Identifiers Only"),
    ]
) class BarcodeSeqLib(SeqLib):
    """
    Class for count data from barcoded sequencing libraries. Designed for
    barcode-only scoring or as a parent class for
    :py:class:`~seqlib.barcodevariant.BcvSeqLib` and
    :py:class:`~seqlib.barcodeid.BcidSeqLib`.
    """

    treeview_class_name = "Barcode SeqLib"

    def __init__(self):
        # Init step handled by VariantSeqLib's init for Barcode-variant
        if type(self).__name__ != "BcvSeqLib":
            SeqLib.__init__(self)
        self.reads = None
        self.revcomp_reads = None
        self.trim_start = None
        self.trim_length = None
        self.barcode_min_count = None
        self.add_label("barcodes")
        self.logger = logging.getLogger("{}.{}".format(__name__, self.__class__))

    def configure(self, cfg):
        """
        Set up the object using the config object *cfg*, usually derived from
        a ``.json`` file.
        """
        SeqLib.configure(self, cfg)
        self.logger = logging.getLogger(
            "{}.{} - {}".format(__name__, self.__class__.__name__, self.name)
        )

        # handle non-FASTQ config options
        try:
            if "min count" in cfg["barcodes"]:
                self.barcode_min_count = int(cfg["barcodes"]["min count"])
            else:
                self.barcode_min_count = 0
        except KeyError as key:
            raise KeyError("Missing required config value {}".format(key), self.name)

        # if counts are specified, copy them later
        # else handle the FASTQ config options and check the files
        if self.counts_file is None:
            self.configure_fastq(cfg)
            try:
                if split_fastq_path(self.reads) is None:
                    raise ValueError(
                        "FASTQ file error: unrecognized file " "extension", self.name
                    )
            except IOError as fqerr:
                raise IOError("FASTQ file error: {}".format(fqerr), self.name)

    def serialize(self):
        """
        Format this object (and its children) as a config object suitable for
        dumping to a config file.
        """
        cfg = SeqLib.serialize(self)

        cfg["barcodes"] = dict()
        if self.barcode_min_count > 0:
            cfg["barcodes"]["min count"] = self.barcode_min_count

        cfg["fastq"] = self.serialize_fastq()

        return cfg

    def configure_fastq(self, cfg):
        """
        Set up the object's FASTQ_ file handling and filtering options.
        """
        try:
            self.reads = cfg["fastq"]["reads"]
            self.revcomp_reads = cfg["fastq"]["reverse"]

            if "start" in cfg["fastq"]:
                self.trim_start = cfg["fastq"]["start"]
            else:
                self.trim_start = 1

            if "length" in cfg["fastq"]:
                self.trim_length = cfg["fastq"]["length"]
            else:
                self.trim_length = sys.maxsize

            self.filters = cfg["fastq"]["filters"]
        except KeyError as key:
            raise KeyError("Missing required config value {}".format(key), self.name)

    def serialize_fastq(self):
        """
        Serialize this object's FASTQ_ file handling and filtering options.
        """
        fastq = {
            "reads": self.reads,
            "reverse": self.revcomp_reads,
            "filters": self.serialize_filters(),
        }
        if self.trim_start > 1:
            fastq["start"] = self.trim_start

        if self.trim_length < sys.maxsize:
            fastq["length"] = self.trim_length

        return fastq

    def counts_from_reads(self):
        """
        Reads the forward or reverse FASTQ_ file (reverse reads are
        reverse-complemented), performs quality-based filtering, and counts
        the barcodes.

        Barcode counts after read-level filtering are stored under
        ``"/raw/barcodes/counts"``.
        """
        df_dict = dict()

        filter_flags = dict()
        for key in self.filters:
            filter_flags[key] = False

        # count all the barcodes
        self.logger.info("Counting barcodes")
        for fqr in read_fastq(self.reads):
            fqr.trim_length(self.trim_length, start=self.trim_start)
            if self.revcomp_reads:
                fqr.revcomp()

            if self.read_quality_filter(fqr):  # passed filtering
                try:
                    df_dict[fqr.sequence.upper()] += 1
                except KeyError:
                    df_dict[fqr.sequence.upper()] = 1

        self.save_counts("barcodes", df_dict, raw=True)
        del df_dict

    def calculate(self):
        """
        Counts the barcodes from the FASTQ file or from the provided counts
        file depending on the config.

        Barcodes that pass the minimum count
        filtering are stored under ``"/main/barcodes/counts"``.

        If ``"/main/barcodes/counts"`` already exists, those will be used
        instead of re-counting.
        """
        if self.check_store("/main/barcodes/counts"):
            return

        # no raw counts present
        if not self.check_store("/raw/barcodes/counts"):
            if self.counts_file is not None:
                self.counts_from_file(self.counts_file)
            else:
                self.counts_from_reads()

        if len(self.labels) == 1:  # only barcodes
            self.save_filtered_counts("barcodes", "count >= self.barcode_min_count")
            self.save_filter_stats() class BcvSeqLib(VariantSeqLib, BarcodeSeqLib):
    """
    Class for counting variant data from barcoded sequencing libraries. 
    Creating a :py:class:`BcvSeqLib` requires a valid *config* 
    object with an ``'barcodes'`` entry and information about the wild type 
    sequence.

    The ``barcode_map`` keyword argument can be used to pass an existing 
    :py:class:`~seqlib.barcodemap.BarcodeMap`. Ensuring this is the 
    right :py:class:`~seqlib.barcodemap.BarcodeMap` is the responsibility 
    of the caller.
    """

    treeview_class_name = "Barcoded Variant SeqLib"

    def __init__(self):
        VariantSeqLib.__init__(self)
        BarcodeSeqLib.__init__(self)
        self.barcode_map = None
        self.logger = logging.getLogger("{}.{}".format(__name__, self.__class__))

    def configure(self, cfg, barcode_map=None):
        """
        Set up the object using the config object *cfg*, usually derived from 
        a ``.json`` file.
        """
        VariantSeqLib.configure(self, cfg)
        BarcodeSeqLib.configure(self, cfg)
        self.logger = logging.getLogger(
            "{}.{} - {}".format(__name__, self.__class__.__name__, self.name)
        )
        try:
            if barcode_map is not None:
                if barcode_map.filename == cfg["barcodes"]["map file"]:
                    self.barcode_map = barcode_map
                else:
                    raise ValueError(
                        "Attempted to assign non-matching barcode map [{}]".format(
                            self.name
                        )
                    )
            else:
                self.barcode_map = BarcodeMap(
                    cfg["barcodes"]["map file"], is_variant=True
                )
        except KeyError as key:
            raise KeyError(
                "Missing required config value {key} [{name}]".format(
                    key=key, name=self.name
                )
            )

    def serialize(self):
        """
        Format this object (and its children) as a config object suitable for dumping to a config file.
        """
        cfg = VariantSeqLib.serialize(self)
        cfg.update(BarcodeSeqLib.serialize(self))

        if self.barcode_map is not None:  # required for creating new objects in GUI
            cfg["barcodes"]["map file"] = self.barcode_map.filename

        return cfg

    def calculate(self):
        """
        Counts the barcodes using :py:meth:`BarcodeSeqLib.count` and combines them into 
        variant counts using the :py:class:`BarcodeMap`.
        """
        if not self.check_store("/main/variants/counts"):
            BarcodeSeqLib.calculate(self)  # count the barcodes
            df_dict = dict()
            barcode_variants = dict()

            self.logger.info("Converting barcodes to variants")
            # store mapped barcodes
            self.save_filtered_counts(
                "barcodes",
                "index in self.barcode_map.keys() & count >= self.barcode_min_count",
            )

            # count variants associated with the barcodes
            max_mut_barcodes = 0
            max_mut_variants = 0
            for bc, count in self.store["/main/barcodes/counts"].iterrows():
                count = count["count"]
                variant = self.barcode_map[bc]
                mutations = self.count_variant(variant)
                if mutations is None:  # variant has too many mutations
                    max_mut_barcodes += 1
                    max_mut_variants += count
                    if self.report_filtered:
                        self.report_filtered_variant(variant, count)
                else:
                    try:
                        df_dict[mutations] += count
                    except KeyError:
                        df_dict[mutations] = count
                    barcode_variants[bc] = mutations

            # save counts, filtering based on the min count
            self.save_counts(
                "variants",
                {k: v for k, v in df_dict.iteritems() if v >= self.variant_min_count},
                raw=False,
            )
            del df_dict

            # write the active subset of the BarcodeMap to the store
            barcodes = barcode_variants.keys()
            barcode_variants = pd.DataFrame(
                {"value": [barcode_variants[bc] for bc in barcodes]}, index=barcodes
            )
            del barcodes
            barcode_variants.sort_values("value", inplace=True)
            self.store.put(
                "/raw/barcodemap",
                barcode_variants,
                data_columns=barcode_variants.columns,
                format="table",
            )
            del barcode_variants

            if self.aligner is not None:
                self.logger.info("Aligned {} variants".format(self.aligner.calls))
                self.aligner_cache = None
            # self.report_filter_stats()
            self.logger.info(
                "Removed {} unique barcodes ({} total variants) "
                "with excess mutations".format(max_mut_barcodes, max_mut_variants)
            )
            self.save_filter_stats()

        self.count_synonymous()

    def make_plots(self):
        """
        Make plots for :py:class:`~seqlib.seqlib.BcvSeqLib` objects.

        Creates plot of the number of barcodes mapping to each variant.
        """
        if self.plots_requested:
            SeqLib.make_plots(self)
            # open the PDF file
            pdf = PdfPages(os.path.join(self.plot_dir, "barcodes_per_variant.pdf"))
            barcodemap_plot(self, pdf)
            pdf.close() class BcidSeqLib(BarcodeSeqLib):
    """
    Class for counting data from barcoded sequencing libraries with non-variant
    identifiers. 
    Creating a :py:class:`BcidSeqLib` requires a valid *config* 
    object with an ``'barcodes'`` entry and information.

    The ``barcode_map`` keyword argument can be used to pass an existing 
    :py:class:`~seqlib.barcodemap.BarcodeMap`. Ensuring this is the 
    right :py:class:`~seqlib.barcodemap.BarcodeMap` is the responsibility 
    of the caller.
    """

    treeview_class_name = "Barcoded ID SeqLib"

    def __init__(self):
        BarcodeSeqLib.__init__(self)
        self.barcode_map = None
        self.identifier_min_count = None
        self.add_label("identifiers")
        self.logger = logging.getLogger("{}.{}".format(__name__, self.__class__))

    def configure(self, cfg, barcode_map=None):
        """
        Set up the object using the config object *cfg*, usually derived from 
        a ``.json`` file.
        """
        BarcodeSeqLib.configure(self, cfg)
        self.logger = logging.getLogger(
            "{}.{} - {}".format(__name__, self.__class__.__name__, self.name)
        )
        try:
            if "min count" in cfg["identifiers"]:
                self.identifier_min_count = int(cfg["identifiers"]["min count"])
            else:
                self.identifier_min_count = 0

            if barcode_map is not None:
                if barcode_map.filename == cfg["barcodes"]["map file"]:
                    self.barcode_map = barcode_map
                else:
                    raise ValueError(
                        "Attempted to assign non-matching barcode map [{}]".format(
                            self.name
                        )
                    )
            else:
                self.barcode_map = BarcodeMap(
                    cfg["barcodes"]["map file"], is_variant=False
                )
        except KeyError as key:
            raise KeyError(
                "Missing required config value {key} [{name}]".format(
                    key=key, name=self.name
                )
            )

    def serialize(self):
        """
        Format this object (and its children) as a config object suitable for dumping to a config file.
        """
        cfg = BarcodeSeqLib.serialize(self)

        cfg["identifiers"] = dict()
        if self.identifier_min_count > 0:
            cfg["identifiers"]["min count"] = self.identifier_min_count

        if self.barcode_map is not None:  # required for creating new objects in GUI
            cfg["barcodes"]["map file"] = self.barcode_map.filename

        return cfg

    def calculate(self):
        """
        Counts the barcodes using :py:meth:`BarcodeSeqLib.count` and combines them into 
        identifier counts using the :py:class:`BarcodeMap`.
        """
        if not self.check_store("/main/identifiers/counts"):
            BarcodeSeqLib.calculate(self)  # count the barcodes
            df_dict = dict()
            barcode_identifiers = dict()

            self.logger.info("Converting barcodes to identifiers")
            # store mapped barcodes
            self.save_filtered_counts(
                "barcodes",
                "index in self.barcode_map.keys() & count >= self.barcode_min_count",
            )

            # count identifiers associated with the barcodes
            for bc, count in self.store["/main/barcodes/counts"].iterrows():
                count = count["count"]
                identifier = self.barcode_map[bc]
                try:
                    df_dict[identifier] += count
                except KeyError:
                    df_dict[identifier] = count
                barcode_identifiers[bc] = identifier

            # save counts, filtering based on the min count
            self.save_counts(
                "identifiers",
                {
                    k: v
                    for k, v in df_dict.iteritems()
                    if v >= self.identifier_min_count
                },
                raw=False,
            )
            del df_dict

            # write the active subset of the BarcodeMap to the store
            barcodes = barcode_identifiers.keys()
            barcode_identifiers = pd.DataFrame(
                {"value": [barcode_identifiers[bc] for bc in barcodes]}, index=barcodes
            )
            del barcodes
            barcode_identifiers.sort_values("value", inplace=True)
            self.store.put(
                "/raw/barcodemap",
                barcode_identifiers,
                data_columns=barcode_identifiers.columns,
                format="table",
            )
            del barcode_identifiers

            # self.report_filter_stats()
            self.save_filter_stats()

    def make_plots(self):
        """
        Make plots for :py:class:`~seqlib.seqlib.BcidSeqLib` objects.

        Creates plot of the number of barcodes mapping to each identifier.
        """
        if self.plots_requested:
            SeqLib.make_plots(self)
            # open the PDF file
            pdf = PdfPages(os.path.join(self.plot_dir, "barcodes_per_identifier.pdf"))
            barcodemap_plot(self, pdf)
            pdf.close() class BasicSeqLib(VariantSeqLib):
    """
    Class for count data from sequencing libraries with a single read for
    each variant. Creating a :py:class:`BasicSeqLib` requires a valid
    *config* object, usually from a ``.json`` configuration file.
    """

    treeview_class_name = "Basic SeqLib"

    def __init__(self):
        VariantSeqLib.__init__(self)
        self.reads = None
        self.revcomp_reads = None
        self.trim_start = None
        self.trim_length = None
        self.logger = logging.getLogger("{}.{}".format(__name__, self.__class__))

    def configure(self, cfg):
        """
        Set up the object using the config object *cfg*, usually derived from
        a ``.json`` file.
        """
        VariantSeqLib.configure(self, cfg)
        self.logger = logging.getLogger(
            "{}.{} - {}".format(__name__, self.__class__.__name__, self.name)
        )

        # if counts are specified, copy them later
        # else handle the FASTQ config options and check the files
        if self.counts_file is None:
            self.configure_fastq(cfg)
            try:
                if split_fastq_path(self.reads) is None:
                    raise IOError(
                        "FASTQ file error: unrecognized extension "
                        "[{}]".format(self.name)
                    )
            except IOError as fqerr:
                raise IOError("FASTQ file error [{}]: {}".format(self.name, fqerr))

    def serialize(self):
        """
        Format this object (and its children) as a config object suitable for
        dumping to a config file.
        """
        cfg = VariantSeqLib.serialize(self)

        cfg["fastq"] = self.serialize_fastq()

        return cfg

    def configure_fastq(self, cfg):
        """
        Set up the object's FASTQ_ file handling and filtering options.
        """
        try:
            self.reads = cfg["fastq"]["reads"]

            if "reverse" in cfg["fastq"]:
                self.revcomp_reads = cfg["fastq"]["reverse"]
            else:
                self.revcomp_reads = False

            if "start" in cfg["fastq"]:
                self.trim_start = cfg["fastq"]["start"]
            else:
                self.trim_start = 1

            if "length" in cfg["fastq"]:
                self.trim_length = cfg["fastq"]["length"]
            else:
                self.trim_length = sys.maxsize

            self.filters = cfg["fastq"]["filters"]
        except KeyError as key:
            raise KeyError(
                "Missing required config value {key} [{name}]"
                "".format(key=key, name=self.name)
            )

    def serialize_fastq(self):
        """
        Serialize this object's FASTQ_ file handling and filtering options.
        """
        fastq = {"filters": self.serialize_filters()}
        fastq["reads"] = self.reads

        if self.revcomp_reads:
            fastq["reverse"] = True
        else:
            fastq["reverse"] = False

        if self.trim_start > 1:
            fastq["start"] = self.trim_start

        if self.trim_length < sys.maxsize:
            fastq["length"] = self.trim_length

        return fastq

    def counts_from_reads(self):
        """
        Reads the forward or reverse FASTQ_ file (reverse reads are
        reverse-complemented), performs quality-based filtering, and counts
        the variants.
        """
        df_dict = dict()

        self.logger.info("Counting variants")
        max_mut_variants = 0
        for fq in read_fastq(self.reads):
            fq.trim_length(self.trim_length, start=self.trim_start)
            if self.revcomp_reads:
                fq.revcomp()

            if self.read_quality_filter(fq):
                mutations = self.count_variant(fq.sequence)
                if mutations is None:  # too many mutations
                    max_mut_variants += 1
                    if self.report_filtered:
                        self.report_filtered_variant(fq.sequence, 1)
                else:
                    try:
                        df_dict[mutations] += 1
                    except KeyError:
                        df_dict[mutations] = 1

        self.save_counts("variants", df_dict, raw=True)
        del df_dict

        if self.aligner is not None:
            self.logger.info("Aligned {} variants".format(self.aligner.calls))
            self.aligner_cache = None
        self.logger.info(
            "Removed {} total variants with excess mutations"
            "".format(max_mut_variants)
        )
        self.save_filter_stats()

    def calculate(self):
        """
        Counts variants from counts file or FASTQ.
        """
        if not self.check_store("/main/variants/counts"):
            if not self.check_store("/raw/variants/counts"):
                if self.counts_file is not None:
                    self.counts_from_file(self.counts_file)
                else:
                    self.counts_from_reads()
            self.save_filtered_counts("variants", "count >= self.variant_min_count")
        self.count_synonymous() class IdOnlySeqLib(SeqLib):
    """
    Class for counting data with non-variant identifiers and no associated
    FASTQ_ data.
    """

    treeview_class_name = "ID-only SeqLib"

    def __init__(self):
        SeqLib.__init__(self)
        self.identifier_min_count = None
        self.add_label("identifiers")
        self.logger = logging.getLogger("{}.{}".format(__name__, self.__class__))

    def configure(self, cfg):
        """
        Set up the object using the config object *cfg*, usually derived from
        a ``.json`` file.
        """
        SeqLib.configure(self, cfg)
        self.logger = logging.getLogger(
            "{}.{} - {}".format(__name__, self.__class__.__name__, self.name)
        )
        try:
            if "min count" in cfg["identifiers"]:
                self.identifier_min_count = int(cfg["identifiers"]["min count"])
            else:
                self.identifier_min_count = 0
        except KeyError as key:
            raise KeyError(
                "Missing required config value {key} [{name}]"
                "".format(key=key, name=self.name)
            )

    def serialize(self):
        """
        Format this object (and its children) as a config object suitable for
        dumping to a config file.
        """
        cfg = SeqLib.serialize(self)

        cfg["identifiers"] = dict()
        if self.identifier_min_count > 0:
            cfg["identifiers"]["min count"] = self.identifier_min_count

        return cfg

    def calculate(self):
        """
        Get the identifier counts from the counts file.
        """
        if not self.check_store("/main/identifiers/counts"):
            if self.counts_file is not None:
                self.counts_from_file(self.counts_file)
            else:
                raise ValueError("Missing counts file [{}]".format(self.name))
            self.save_filtered_counts(
                "identifiers", "count >= self.identifier_min_count"
            ) class OverlapSeqLib(VariantSeqLib):
    """
    Class for count data from sequencing libraries with overlapping paired-end 
    reads for each variant. Creating a 
    :py:class:`~seqlib.overlap.OverlapSeqLib` requires a valid *config* object 
    with an ``'overlap'`` entry.

    The ``"fastq"`` config entry must contain two read files, with the keys 
    ``"forward"`` and ``"reverse"``. Information about how to combine these 
    reads is in the ``"overlap"`` config entry.

    The ``"overlap"`` config entry contains the following keys:

    * ``"forward start"`` --- position in the forward read where the \
        overlapping region begins
    * ``"reverse start"`` --- position in the reverse read where the \
        overlapping region begins (before being reverse-complemented)
    * ``"length"`` --- number of bases in the overlapping region
    * ``"max mismatches"`` --- maximum number of mismatches tolerated in the \
        overlapping region before discarding the read
    * ``"overlap only"`` --- whether to trim the merged read to contain only \
        the overlapping region (optional, defaults to ``False``)

    Here is a schematic of the case in the above JSON example::

        forward ---> 1   
                     CGACGCAAGGA
                       |||||||||
                       ACTCCTTGCGTCG
                                   1 <--- reverse

    Note that the merged sequence is identical to the wild type sequence given 
    in the JSON file.
    """

    treeview_class_name = "Overlap SeqLib"

    def __init__(self):
        VariantSeqLib.__init__(self)
        self.forward = None
        self.reverse = None
        self.fwd_start = None
        self.rev_start = None
        self.overlap_length = None
        self.trim = None
        self.max_overlap_mismatches = None
        self.merge_mismatches = None
        self.default_filters.update({"merge failure": True})
        self.default_filters.update({"remove unresolvable": False})
        self.logger = logging.getLogger("{}.{}".format(__name__, self.__class__))

    def configure(self, cfg):
        """
        Set up the object using the config object *cfg*, usually derived from 
        a ``.json`` file.
        """
        VariantSeqLib.configure(self, cfg)
        self.logger = logging.getLogger(
            "{}.{} - {}".format(__name__, self.__class__.__name__, self.name)
        )

        # if counts are specified, copy them later
        # else handle the FASTQ config options and check the files
        if self.counts_file is None:
            self.configure_fastq(cfg)
            try:
                self.fwd_start = int(cfg["overlap"]["forward start"])
                self.rev_start = int(cfg["overlap"]["reverse start"])
                self.overlap_length = int(cfg["overlap"]["length"])
                self.trim = cfg["overlap"]["trim"]
                self.max_overlap_mismatches = int(cfg["overlap"]["max mismatches"])

                forward_error = False
                reverse_error = False
                if split_fastq_path(self.forward) is None:
                    forward_error = True
                if split_fastq_path(self.reverse) is None:
                    reverse_error = True
                if forward_error and reverse_error:
                    raise IOError(
                        "FASTQ file error: unrecognized extension (forward and reverse) [{}]".format(
                            self.name
                        )
                    )
                elif forward_error:
                    raise IOError(
                        "FASTQ file error: unrecognized extension (forward) [{}]".format(
                            self.name
                        )
                    )
                elif reverse_error:
                    raise IOError(
                        "FASTQ file error: unrecognized extension (reverse) [{}]".format(
                            self.name
                        )
                    )
            except IOError as fqerr:
                raise IOError("FASTQ file error [{}]: {}".format(self.name, fqerr))
            except KeyError as key:
                raise KeyError(
                    "Missing required config value {key} [{name}]".format(
                        key=key, name=self.name
                    )
                )
            except ValueError as value:
                raise ValueError(
                    "Invalid parameter value {value} [{name}]".format(
                        value=value, name=self.name
                    )
                )

    def serialize(self):
        """
        Format this object (and its children) as a config object suitable for dumping to a config file.
        """
        cfg = VariantSeqLib.serialize(self)

        cfg["fastq"] = self.serialize_fastq()
        cfg["overlap"] = {
            "forward start": self.fwd_start,
            "reverse start": self.rev_start,
            "length": self.overlap_length,
            "trim": self.trim,
            "max mismatches": self.max_overlap_mismatches,
        }

        return cfg

    def configure_fastq(self, cfg):
        """
        Set up the object's FASTQ_ file handling and filtering options.
        """
        try:
            self.forward = cfg["fastq"]["forward reads"]
            self.reverse = cfg["fastq"]["reverse reads"]

            if "merge failure" in cfg["fastq"]["filters"]:
                raise ValueError(
                    "'merge failure' is not user-configurable [{}]".format(self.name)
                )
            self.filters = cfg["fastq"]["filters"]
        except KeyError as key:
            raise KeyError(
                "Missing required config value {key} [{name}]".format(
                    key=key, name=self.name
                )
            )

    def serialize_fastq(self):
        """
        Serialize this object's FASTQ_ file handling and filtering options.
        """
        fastq = {
            "forward reads": self.forward,
            "reverse reads": self.reverse,
            "filters": self.serialize_filters(),
        }

        return fastq

    def merge_reads(self, fwd, rev):
        """
        Combines the *fwd* and *rev* :py:class:`~fqread.FQRead` objects into a 
        single :py:class:`~fqread.FQRead` with the same header information as 
        *fwd*. Mismatches are resolved by taking the highest quality base. If 
        discrepant bases have the same quality value, this position is 
        unresolvable and an ``'X'`` is inserted. Quality values in the 
        resulting :py:class:`~fqread.FQRead` are the maximum quality for the 
        given base at that position. Returns ``None`` if the maximum number of 
        mismatches in the overlap region is exceded.
        """
        rev.revcomp()

        # print(fwd.sequence, "-" * (self.rev_start - 1), sep="")
        # print("-" * (self.fwd_start - 1), rev.sequence, sep="")
        rev_extra_start = len(rev) - self.rev_start + 1
        fwd_end = self.fwd_start + self.overlap_length - 1
        merge = FQRead(
            header=fwd.header + "|" + rev.header,
            sequence="A",
            header2=fwd.header2 + "|" + rev.header2,
            quality="#",
            qbase=fwd.qbase,
        )
        merge.sequence = list(fwd.sequence[:fwd_end] + rev.sequence[rev_extra_start:])
        merge.quality = fwd.quality[:fwd_end] + rev.quality[rev_extra_start:]

        mismatches = 0
        first = True
        for i in xrange(self.overlap_length):
            a = self.fwd_start - 1 + i
            b = len(rev) - self.rev_start - self.overlap_length + i + 1
            try:
                if fwd.sequence[a] == rev.sequence[b]:
                    # take the highest quality value
                    if rev.quality[b] > fwd.quality[a]:
                        merge.quality[a] = rev.quality[b]
                else:
                    if fwd.quality[a] == rev.quality[b]:
                        merge.sequence[a] = "X"  # unresolvable
                        self.merge_mismatches.iloc[i]["unresolved"] += 1
                    elif rev.quality[b] > fwd.quality[a]:
                        merge.sequence[a] = rev.sequence[b]
                        merge.quality[a] = rev.quality[b]
                        self.merge_mismatches.iloc[i]["resolved"] += 1
                    else:
                        # overlap region already same as fwd
                        self.merge_mismatches.iloc[i]["resolved"] += 1
                    mismatches += 1
                    if first:
                        self.merge_mismatches.iloc[i]["first"] += 1
                        first = False
            except IndexError:
                raise IndexError(
                    "Failed to calculate overlap (a={a}, len(a)={lena}, b={b}, len(b)={lenb}) [{name}]".format(
                        a=a,
                        b=b,
                        lena=len(fwd.sequence),
                        lenb=len(rev.sequence),
                        name=self.name,
                    )
                )

        if mismatches > self.max_overlap_mismatches:
            return None  # merge failed

        merge.sequence = "".join(merge.sequence)
        if self.trim:
            merge.trim_length(self.overlap_length, self.fwd_start)
        return merge

    def counts_from_reads(self):
        df_dict = dict()

        self.merge_mismatches = pd.DataFrame(
            data=0,
            index=[
                x + self.fwd_start + self.wt.dna_offset
                for x in xrange(0, self.overlap_length)
            ],
            columns=["resolved", "unresolved", "first"],
        )

        self.logger.info("Counting variants")
        max_mut_variants = 0
        for fwd, rev in read_fastq_multi([self.forward, self.reverse]):
            # filter on chastity before merge
            chaste = True
            if self.filters["chastity"]:
                if not fwd.is_chaste():
                    chaste = False
                    if self.report_filtered:
                        self.report_filtered_read(fwd, "chastity")
                if not rev.is_chaste():
                    chaste = False
                    if self.report_filtered:
                        self.report_filtered_read(rev, "chastity")
                if not chaste:
                    self.filter_stats["chastity"] += 1
                    self.filter_stats["total"] += 1
                    continue

            merge = self.merge_reads(fwd, rev)
            if merge is None:  # merge failed
                self.filter_stats["merge failure"] += 1
                self.filter_stats["total"] += 1
                if self.report_filtered:
                    self.report_filtered_read(fwd, filter_flags)
                    self.report_filtered_read(rev, filter_flags)
            else:
                if self.read_quality_filter(merge):
                    mutations = self.count_variant(merge.sequence)
                    if mutations is None:  # merge read has too many mutations
                        max_mut_variants += 1
                        if self.report_filtered:
                            self.report_filtered_variant(merge.sequence, 1)
                    else:
                        try:
                            df_dict[mutations] += 1
                        except KeyError:
                            df_dict[mutations] = 1

        self.store.put(
            "/raw/overlap_mismatches",
            self.merge_mismatches,
            format="table",
            data_columns=self.merge_mismatches.columns,
        )
        self.merge_mismatches = None
        self.save_counts("variants", df_dict, raw=True)
        del df_dict

        if self.aligner is not None:
            self.logger.info("Aligned {} variants".format(self.aligner.calls))
            self.aligner_cache = None
        self.logger.info(
            "Removed {} total variants with excess mutations"
            "".format(max_mut_variants)
        )
        self.save_filter_stats()

    def calculate(self):
        """
        Reads the forward and reverse reads, merges them, performs 
        quality-based filtering, and counts the variants.
        """
        if not self.check_store("/main/variants/counts"):
            if not self.check_store("/raw/variants/counts"):
                if self.counts_file is not None:
                    self.counts_from_file(self.counts_file)
                else:  # count everything
                    self.counts_from_reads()
            self.save_filtered_counts("variants", "count >= self.variant_min_count")

        self.count_synonymous()

    def make_plots(self):
        """
        Make plots for :py:class:`~seqlib.seqlib.OverlapSeqLib` objects.

        Creates plots of the location of merged read mismatches.
        """
        if self.plots_requested:
            SeqLib.make_plots(self)
            pdf = PdfPages(os.path.join(self.plot_dir, "overlap_mismatches.pdf"))
            overlap_merge_plot(self, pdf)
            pdf.close() class Selection(StoreManager):
    """
    Class for a single selection replicate, consisting of multiple 
    timepoints. This class coordinates :py:class:`~seqlib.seqlib.SeqLib` 
    objects.
    """

    store_suffix = "sel"
    treeview_class_name = "Selection"

    def __init__(self):
        StoreManager.__init__(self)
        self.libraries = dict()
        self.barcode_maps = dict()
        self._wt = None
        self.logger = logging.getLogger("{}.{}".format(__name__, self.__class__))

    def _children(self):
        """
        Return the :py:class:`~seqlib.seqlib.SeqLib` objects as a list, 
        sorted by timepoint and then by name.
        """
        libs = list()
        for tp in self.timepoints:
            libs.extend(sorted(self.libraries[tp], key=lambda x: x.name))
        return libs

    def remove_child_id(self, tree_id):
        """
        Remove the reference to a :py:class:`~seqlib.seqlib.SeqLib` with 
        Treeview id *tree_id*. Deletes empty time points.
        """
        empty = None
        for tp in self.libraries:
            tp_ids = [lib.treeview_id for lib in self.libraries[tp]]
            if tree_id in tp_ids:
                del self.libraries[tp][tp_ids.index(tree_id)]
                if len(self.libraries[tp]) == 0:
                    empty = tp
                break  # found the id, stop looking
        if empty is not None:
            del self.libraries[empty]

    @property
    def timepoints(self):
        return sorted(self.libraries.keys())

    @property
    def wt(self):
        if self.has_wt_sequence():
            if self._wt is None:
                self._wt = self.children[0].wt.duplicate(self.name)
            return self._wt
        else:
            if self._wt is not None:
                raise ValueError(
                    "Selection should not contain wild type sequence [{}]".format(
                        self.name
                    )
                )
            else:
                return None

    def configure(self, cfg, configure_children=True):
        """
        Set up the :py:class:`~selection.Selection` using the *cfg* object, 
        usually from a ``.json`` configuration file.

        If *configure_children* is false, do not configure the children in 
        *cfg*.
        """
        StoreManager.configure(self, cfg)
        self.logger = logging.getLogger(
            "{}.{} - {}".format(__name__, self.__class__.__name__, self.name)
        )
        if configure_children:
            if "libraries" not in cfg:
                raise KeyError(
                    "Missing required config value {} [{}]".format(
                        "libraries", self.name
                    )
                )

            for lib_cfg in cfg["libraries"]:
                libtype = seqlib_type(lib_cfg)
                if libtype is None:
                    raise ValueError("Unrecognized SeqLib config")
                elif libtype in ("BcvSeqLib", "BcidSeqLib"):
                    lib = SEQLIB_CLASSES[libtype]()
                    # don't re-parse the barcode maps if possible
                    mapfile = lib_cfg["barcodes"]["map file"]
                    if mapfile in self.barcode_maps.keys():
                        lib.configure(lib_cfg, barcode_map=self.barcode_maps[mapfile])
                    else:
                        lib.configure(lib_cfg)
                        self.barcode_maps[mapfile] = lib.barcode_map
                    self.add_child(lib)
                else:
                    # requires that the SeqLib derived classes be imported into the
                    # module namespace using "from x import y" style
                    lib = SEQLIB_CLASSES[libtype]()
                    lib.configure(lib_cfg)
                    self.add_child(lib)

    def validate(self):
        """
        Raises an informative ``ValueError`` if the time points in the analysis are not suitable.

        Calls validate method on all child SeqLibs.
        """
        # check the time points
        if 0 not in self.timepoints:
            raise ValueError("Missing timepoint 0 [{}]".format(self.name))
        if self.timepoints[0] != 0:
            raise ValueError("Invalid negative timepoint [{}]".format(self.name))
        if len(self.timepoints) < 2:
            raise ValueError("Multiple timepoints required [{}]".format(self.name))
        elif len(self.timepoints) < 3 and self.scoring_method in ("WLS", "OLS"):
            raise ValueError(
                "Insufficient number of timepoints for regression scoring [{}]".format(
                    self.name
                )
            )

        # check the wild type sequences
        if self.has_wt_sequence():
            for child in self.children[1:]:
                if self.children[0].wt != child.wt:
                    self.logger.warning("Inconsistent wild type sequences")
                    break

        # check that we're not doing wild type normalization on something with no wild type
        # if not self.has_wt_sequence() and self.logr_method == "wt":
        #    raise ValueError("No wild type sequence for wild type normalization [{}]".format(self.name))

        # validate children
        for child in self.children:
            child.validate()

    def serialize(self):
        """
        Format this object (and its children) as a config object suitable for dumping to a config file.
        """
        cfg = StoreManager.serialize(self)
        cfg["libraries"] = [child.serialize() for child in self.children]
        return cfg

    def add_child(self, child):
        if child.name in self.child_names():
            raise ValueError(
                "Non-unique sequencing library name '{}' [{}]".format(
                    child.name, self.name
                )
            )

        child.parent = self

        # add it to the libraries dictionary
        try:
            self.libraries[child.timepoint].append(child)
        except KeyError:
            self.libraries[child.timepoint] = [child]

    def is_barcodevariant(self):
        """
        Return ``True`` if all :py:class:`~seqlib.seqlib.SeqLib` in the 
        :py:class:`~selection.Selection` are 
        :py:class:`~barcodevariant.BcvSeqLib` objects with 
        the same barcode map, else ``False``.
        """
        return (
            all(isinstance(lib, BcvSeqLib) for lib in self.children)
            and len(self.barcode_maps.keys()) == 1
        )

    def is_barcodeid(self):
        """
        Return ``True`` if all :py:class:`~seqlib.SeqLib` in the 
        :py:class:`~selection.Selection` are 
        :py:class:`~barcodeid.BcidSeqLib` objects with 
        the same barcode map, else ``False``.
        """
        return (
            all(isinstance(lib, BcidSeqLib) for lib in self.children)
            and len(self.barcode_maps.keys()) == 1
        )

    def is_coding(self):
        """
        Return ``True`` if the all :py:class:`~seqlib.seqlib.SeqLib` in the 
        :py:class:`~selection.Selection` count protein-coding variants, else 
        ``False``.
        """
        return all(x.is_coding() for x in self.children)

    def has_wt_sequence(self):
        """
        Return ``True`` if the all :py:class:`~seqlib.seqlib.SeqLib` in the 
        :py:class:`~selection.Selection` have a wild type sequence, else 
        ``False``.
        """
        return all(x.has_wt_sequence() for x in self.children)

    def merge_counts_unfiltered(self, label):
        """
        Counts :py:class:`~seqlib.seqlib.SeqLib` objects and tabulates counts 
        for each timepoint. :py:class:`~seqlib.seqlib.SeqLib` objects from 
        the same timepoint are combined by summing the counts.

        Stores the unfiltered counts under ``/main/label/counts_unfiltered``.
        """
        if self.check_store("/main/{}/counts_unfiltered".format(label)):
            return

        # calculate counts for each SeqLib
        self.logger.info("Counting for each time point ({})".format(label))
        for lib in self.children:
            lib.calculate()

        # combine all libraries for a given timepoint
        self.logger.info("Aggregating SeqLib data")

        destination = "/main/{}/counts_unfiltered".format(label)
        if destination in self.store.keys():
            # need to remove the current destination table because we are using append
            # append is required because it takes the "min_itemsize" argument, and put doesn't
            self.logger.info("Replacing existing '{}'".format(destination))
            self.store.remove(destination)

        # seqlib count table name for this element type
        lib_table = "/main/{}/counts".format(label)

        # create an index of all elements in the analysis
        complete_index = pd.Index([])
        for tp in self.timepoints:
            for lib in self.libraries[tp]:
                complete_index = complete_index.union(
                    pd.Index(lib.store.select_column(lib_table, "index"))
                )
        self.logger.info(
            "Created shared index for count data ({} {})".format(
                len(complete_index), label
            )
        )

        # min_itemsize value
        max_index_length = complete_index.map(len).max()

        # perform operation in chunks
        tp_frame = None
        for i in xrange(0, len(complete_index), self.chunksize):
            # don't duplicate the index if the chunksize is large
            if self.chunksize < len(complete_index):
                index_chunk = complete_index[i : i + self.chunksize]
            else:
                index_chunk = complete_index
            self.logger.info(
                "Merging counts for chunk {} ({} rows)".format(
                    i / self.chunksize + 1, len(index_chunk)
                )
            )

            for tp in self.timepoints:
                c = self.libraries[tp][0].store.select(lib_table, "index = index_chunk")
                for lib in self.libraries[tp][1:]:
                    c = c.add(
                        lib.store.select(lib_table, "index = index_chunk"), fill_value=0
                    )
                c.columns = ["c_{}".format(tp)]
                if tp == 0:
                    tp_frame = c
                else:
                    tp_frame = tp_frame.join(c, how="outer")

            # save the unfiltered counts
            if "/main/{}/counts_unfiltered".format(label) not in self.store:
                self.store.append(
                    "/main/{}/counts_unfiltered".format(label),
                    tp_frame.astype(float),
                    min_itemsize={"index": max_index_length},
                    data_columns=list(tp_frame.columns),
                )
            else:
                self.store.append(
                    "/main/{}/counts_unfiltered".format(label), tp_frame.astype(float)
                )

    def filter_counts(self, label):
        """
        Converts unfiltered counts stored in ``/main/label/counts_unfiltered`` 
        into filtered counts calculated from complete cases (elements with a 
        non-zero count in each time point).

        For the most basic element type (variant or barcode, depending on the 
        experimental design), the result of this operation simply drops any 
        rows that have missing counts. For other element types, such as 
        synonymous variants, the counts are re-aggregated using only the 
        complete cases in the underlying element type.
        """
        if (self.is_barcodeid() or self.is_barcodevariant()) and label != "barcodes":
            # calculate proper combined counts
            # df = self.store.select("/main/barcodes/counts") # this should exist because of the order of label calculations
            # redo the barcode->variant/id mapping with the filtered counts
            # NOT YET IMPLEMENTED
            df = self.store.select(
                "/main/{}/counts_unfiltered".format(label)
            )  # just do this for now
        else:
            df = self.store.select("/main/{}/counts_unfiltered".format(label))
        df.dropna(axis="index", how="any", inplace=True)
        self.store.put(
            "/main/{}/counts".format(label),
            df.astype(float),
            format="table",
            data_columns=df.columns,
        )

    def combine_barcode_maps(self):
        if self.check_store("/main/barcodemap"):
            return

        bcm = None
        for lib in self.children:
            if bcm is None:
                bcm = lib.store["/raw/barcodemap"]
            else:
                bcm = bcm.join(
                    lib.store["/raw/barcodemap"], rsuffix=".drop", how="outer"
                )
                new = bcm.loc[pd.isnull(bcm)["value"]]
                bcm.loc[new.index, "value"] = new["value.drop"]
                bcm.drop("value.drop", axis="columns", inplace=True)
        bcm.sort_values("value", inplace=True)
        self.store.put(
            "/main/barcodemap", bcm, format="table", data_columns=bcm.columns
        )

    def calculate(self):
        """
        Wrapper method to calculate counts and enrichment scores 
        for all data in the :py:class:`~selection.Selection`.
        """
        if len(self.labels) == 0:
            raise ValueError(
                "No data present across all sequencing libraries [{}]".format(self.name)
            )
        for label in self.labels:
            self.merge_counts_unfiltered(label)
            self.filter_counts(label)
        if self.is_barcodevariant() or self.is_barcodeid():
            self.combine_barcode_maps()
        if self.scoring_method == "counts":
            pass
        elif self.scoring_method == "ratios":
            for label in self.labels:
                self.calc_ratios(label)
        elif self.scoring_method == "simple":
            for label in self.labels:
                self.calc_simple_ratios(label)
        elif self.scoring_method in ("WLS", "OLS"):
            if len(self.timepoints) <= 2:
                raise ValueError(
                    "Regression-based scoring requires three or more time points."
                )
            for label in self.labels:
                self.calc_log_ratios(label)
                if self.scoring_method == "WLS":
                    self.calc_weights(label)
                self.calc_regression(label)
        else:
            raise ValueError(
                'Invalid scoring method "{}" [{}]'.format(
                    self.scoring_method, self.name
                )
            )

        if self.scoring_method in ("ratios", "WLS", "OLS") and self.component_outliers:
            if self.is_barcodevariant() or self.is_barcodeid():
                self.calc_outliers("barcodes")
            if self.is_coding():
                self.calc_outliers("variants")

    def calc_simple_ratios(self, label):
        """
        Calculate simplified (original Enrich) ratios scores. This method does not produce standard errors.
        """
        if self.check_store("/main/{}/scores".format(label)):
            return

        self.logger.info("Calculating simple ratios ({})".format(label))
        c_last = "c_{}".format(self.timepoints[-1])
        df = self.store.select(
            "/main/{}/counts".format(label), "columns in ['c_0', c_last]"
        )

        # perform operations on the numpy values of the data frame for easier broadcasting
        ratios = (df[c_last].values.astype("float") / df[c_last].sum(axis="index")) / (
            df["c_0"].values.astype("float") / df["c_0"].sum(axis="index")
        )
        # make it a data frame again
        ratios = pd.DataFrame(data=ratios, index=df.index, columns=["ratio"])
        ratios["score"] = np.log2(ratios["ratio"])
        ratios["SE"] = np.nan
        ratios = ratios[["score", "SE", "ratio"]]  # re-order columns

        self.store.put(
            "/main/{}/scores".format(label),
            ratios,
            format="table",
            data_columns=ratios.columns,
        )

    def calc_ratios(self, label):
        """
        Calculate frequency ratios and standard errors between the last timepoint and the input.

        Ratios can be calculated using one of three methods:

        * wt
        * complete
        * full

        """
        if self.check_store("/main/{}/scores".format(label)):
            return

        self.logger.info("Calculating ratios ({})".format(label))
        c_last = "c_{}".format(self.timepoints[-1])
        df = self.store.select(
            "/main/{}/counts".format(label), "columns in ['c_0', c_last]"
        )

        if self.logr_method == "wt":
            if "variants" in self.labels:
                wt_label = "variants"
            elif "identifiers" in self.labels:
                wt_label = "identifiers"
            else:
                raise ValueError(
                    "Failed to use wild type log ratio method, suitable data table not present [{}]".format(
                        self.name
                    )
                )
            shared_counts = self.store.select(
                "/main/{}/counts".format(wt_label),
                "columns in ['c_0', c_last] & index='{}'".format(WILD_TYPE_VARIANT),
            )
            if len(shared_counts) == 0:  # wild type not found
                raise ValueError(
                    "Failed to use wild type log ratio method, wild type sequence not present [{}]".format(
                        self.name
                    )
                )
            shared_counts = shared_counts.values + 0.5
        elif self.logr_method == "complete":
            shared_counts = (
                self.store.select(
                    "/main/{}/counts".format(label), "columns in ['c_0', c_last]"
                )
                .sum(axis="index")
                .values
                + 0.5
            )
        elif self.logr_method == "full":
            shared_counts = (
                self.store.select(
                    "/main/{}/counts_unfiltered".format(label),
                    "columns in ['c_0', c_last]",
                )
                .sum(axis="index", skipna=True)
                .values
                + 0.5
            )
        else:
            raise ValueError(
                'Invalid log ratio method "{}" [{}]'.format(self.logr_method, self.name)
            )

        ratios = np.log(df[["c_0", c_last]].values + 0.5) - np.log(shared_counts)
        ratios = ratios[:, 1] - ratios[:, 0]  # selected - input

        ratios = pd.DataFrame(data=ratios, index=df.index, columns=["logratio"])

        shared_variance = np.sum(1.0 / shared_counts)
        ratios["variance"] = (
            np.sum(1.0 / (df[["c_0", c_last]].values + 0.5), axis=1) + shared_variance
        )

        ratios["score"] = ratios["logratio"]
        ratios["SE"] = np.sqrt(ratios["variance"])
        ratios = ratios[["score", "SE", "logratio", "variance"]]  # re-order columns

        self.store.put(
            "/main/{}/scores".format(label),
            ratios,
            format="table",
            data_columns=ratios.columns,
        )

    def calc_log_ratios(self, label):
        """
        Calculate the log ratios that will be fit using the linear model.
        """
        if self.check_store("/main/{}/log_ratios".format(label)):
            return

        self.logger.info("Calculating log ratios ({})".format(label))
        ratios = self.store.select("/main/{}/counts".format(label))
        index = ratios.index
        c_n = ["c_{}".format(x) for x in self.timepoints]
        ratios = np.log(ratios + 0.5)

        # perform operations on the numpy values of the data frame for easier broadcasting
        ratios = ratios[c_n].values
        if self.logr_method == "wt":
            if "variants" in self.labels:
                wt_label = "variants"
            elif "identifiers" in self.labels:
                wt_label = "identifiers"
            else:
                raise ValueError(
                    "Failed to use wild type log ratio method, suitable data table not present [{}]".format(
                        self.name
                    )
                )
            wt_counts = self.store.select(
                "/main/{}/counts".format(wt_label),
                "columns=c_n & index='{}'".format(WILD_TYPE_VARIANT),
            )
            if len(wt_counts) == 0:  # wild type not found
                raise ValueError(
                    "Failed to use wild type log ratio method, wild type sequence not present [{}]".format(
                        self.name
                    )
                )
            ratios = ratios - np.log(wt_counts.values + 0.5)
        elif self.logr_method == "complete":
            ratios = ratios - np.log(
                self.store.select("/main/{}/counts".format(label), "columns=c_n")
                .sum(axis="index")
                .values
                + 0.5
            )
        elif self.logr_method == "full":
            ratios = ratios - np.log(
                self.store.select(
                    "/main/{}/counts_unfiltered".format(label), "columns=c_n"
                )
                .sum(axis="index", skipna=True)
                .values
                + 0.5
            )
        else:
            raise ValueError(
                'Invalid log ratio method "{}" [{}]'.format(self.logr_method, self.name)
            )

        # make it a data frame again
        ratios = pd.DataFrame(
            data=ratios,
            index=index,
            columns=["L_{}".format(x) for x in self.timepoints],
        )
        self.store.put(
            "/main/{}/log_ratios".format(label),
            ratios,
            format="table",
            data_columns=ratios.columns,
        )

    def calc_weights(self, label):
        """
        Calculate the regression weights (1 / variance).
        """
        if self.check_store("/main/{}/weights".format(label)):
            return

        self.logger.info("Calculating regression weights ({})".format(label))
        variances = self.store.select("/main/{}/counts".format(label))
        c_n = ["c_{}".format(x) for x in self.timepoints]
        index = variances.index
        # perform operations on the numpy values of the data frame for easier broadcasting
        # variances = 1.0 / (variances[c_n].values + 0.5) + 1.0 / (variances[['c_0']].values + 0.5)
        variances = 1.0 / (variances[c_n].values + 0.5)
        if self.logr_method == "wt":
            if "variants" in self.labels:
                wt_label = "variants"
            elif "identifiers" in self.labels:
                wt_label = "identifiers"
            else:
                raise ValueError(
                    "Failed to use wild type log ratio method, suitable data table not present [{}]".format(
                        self.name
                    )
                )
            wt_counts = self.store.select(
                "/main/{}/counts".format(wt_label),
                "columns=c_n & index='{}'".format(WILD_TYPE_VARIANT),
            )
            if len(wt_counts) == 0:  # wild type not found
                raise ValueError(
                    "Failed to use wild type log ratio method, wild type sequence not present [{}]".format(
                        self.name
                    )
                )
            variances = variances + 1.0 / (wt_counts.values + 0.5)
        elif self.logr_method == "complete":
            variances = variances + 1.0 / (
                self.store.select("/main/{}/counts".format(label), "columns=c_n")
                .sum(axis="index")
                .values
                + 0.5
            )
        elif self.logr_method == "full":
            variances = variances + 1.0 / (
                self.store.select(
                    "/main/{}/counts_unfiltered".format(label), "columns=c_n"
                )
                .sum(axis="index", skipna=True)
                .values
                + 0.5
            )
        else:
            raise ValueError(
                'Invalid log ratio method "{}" [{}]'.format(self.logr_method, self.name)
            )
        variances = 1.0 / variances  # weights are reciprocal of variances
        # make it a data frame again
        variances = pd.DataFrame(
            data=variances,
            index=index,
            columns=["W_{}".format(x) for x in self.timepoints],
        )
        self.store.put(
            "/main/{}/weights".format(label),
            variances,
            format="table",
            data_columns=variances.columns,
        )

    def calc_regression(self, label):
        """
        Calculate least squares regression for *label*. If *weighted* is ``True``, calculates weighted least squares; else ordinary least squares.

        Regression results are stored in ``'/main/label/scores'``

        """
        if self.check_store("/main/{}/scores".format(label)):
            return
        elif "/main/{}/scores".format(label) in self.store.keys():
            # need to remove the current keys because we are using append
            self.store.remove("/main/{}/scores".format(label))

        self.logger.info(
            "Calculating {} regression coefficients ({})".format(
                self.scoring_method, label
            )
        )
        # append is required because it takes the "min_itemsize" argument, and put doesn't
        longest = (
            self.store.select("/main/{}/log_ratios".format(label), "columns='index'")
            .index.map(len)
            .max()
        )
        chunk = 1
        if self.scoring_method == "WLS":
            for data in self.store.select_as_multiple(
                ["/main/{}/log_ratios".format(label), "/main/{}/weights".format(label)],
                chunksize=self.chunksize,
            ):
                self.logger.info(
                    "Calculating weighted least squares for chunk {} ({} rows)".format(
                        chunk, len(data.index)
                    )
                )
                result = data.apply(
                    regression_apply, args=[self.timepoints, True], axis="columns"
                )
                self.store.append(
                    "/main/{}/scores".format(label),
                    result,
                    min_itemsize={"index": longest},
                )
                chunk += 1
        elif self.scoring_method == "OLS":
            for data in self.store.select(
                "/main/{}/log_ratios".format(label), chunksize=self.chunksize
            ):
                self.logger.info(
                    "Calculating ordinary least squares for chunk {} ({} rows)".format(
                        chunk, len(data.index)
                    )
                )
                result = data.apply(
                    regression_apply, args=[self.timepoints, False], axis="columns"
                )
                self.store.append(
                    "/main/{}/scores".format(label),
                    result,
                    min_itemsize={"index": longest},
                )
                chunk += 1
        else:
            raise ValueError(
                'Invalid regression scoring method "{}" [{}]'.format(
                    self.scoring_method, self.name
                )
            )

        # need to read from the file, calculate percentiles, and rewrite it
        self.logger.info(
            "Calculating slope standard error percentiles ({})".format(label)
        )
        data = self.store["/main/{}/scores".format(label)]
        data["score"] = data["slope"]
        data["SE"] = data["SE_slope"]
        data["SE_pctile"] = [
            stats.percentileofscore(data["SE"], x, "weak") for x in data["SE"]
        ]
        data = data[
            [
                "score",
                "SE",
                "SE_pctile",
                "slope",
                "intercept",
                "SE_slope",
                "t",
                "pvalue_raw",
            ]
        ]  # reorder columns
        self.store.put(
            "/main/{}/scores".format(label),
            data,
            format="table",
            data_columns=data.columns,
        )

    def wt_plot(self, pdf):
        """
        Create a plot of the linear fit of the wild type variant.

        *pdf* is an open PdfPages instance.

        Only created for selections that use WLS or OLS scoring and have a wild type specified. 
        Uses :py:func:`~plots.fit_axes` for the plotting.
        """
        self.logger.info("Creating wild type fit plot")

        # get the data and calculate log ratios
        if "variants" in self.labels:
            wt_label = "variants"
        elif "identifiers" in self.labels:
            wt_label = "identifiers"
        else:
            raise ValueError(
                "Invalid selection type for wild type fit plot [{}]".format(self.name)
            )
        data = self.store.select(
            "/main/{}/counts".format(wt_label),
            where='index = "{}"'.format(WILD_TYPE_VARIANT),
        ).ix[0]
        sums = self.store["/main/{}/counts".format(wt_label)].sum(
            axis="index"
        )  # sum of complete cases (N')
        yvalues = np.log(data + 0.5) - np.log(sums + 0.5)
        xvalues = [tp / float(max(self.timepoints)) for tp in self.timepoints]

        # fit the line
        X = sm.add_constant(xvalues)  # fit intercept
        if self.scoring_method == "WLS":
            W = 1 / (1 / (data + 0.5) + 1 / (sums + 0.5))
            fit = sm.WLS(yvalues, X, weights=W).fit()
        elif self.scoring_method == "OLS":
            fit = sm.OLS(yvalues, X).fit()
        else:
            raise ValueError(
                'Invalid regression scoring method "{}" [{}]'.format(
                    self.scoring_method, self.name
                )
            )
        intercept, slope = fit.params
        slope_se = fit.bse["x1"]

        # make the plot
        fig, ax = plt.subplots()
        fig.set_tight_layout(True)
        fit_axes(ax, xvalues, yvalues, slope, intercept, xlabels=self.timepoints)
        fit_axes_text(ax, cornertext="Slope {:3.2f}\nSE {:.1f}".format(slope, slope_se))
        ax.set_title("Wild Type Shape\n{}".format(self.name))
        ax.set_ylabel("Log Ratio (Complete Cases)")

        pdf.savefig(fig)
        plt.close(fig)

    def se_pctile_plot(self, label, pdf):
        """
        Create plots of the linear fit of 21 selected variants, evenly distributed based on their standard error percentiles (0, 5, 10, ..., 100).

        *label* is the data label (barcode, variant, etc.)

        *pdf* is an open PdfPages instance.

        Uses :py:func:`~plots.fit_axes` for the plotting.
        """
        self.logger.info("Creating representative fit plots ({})".format(label))

        se_data = self.store.select(
            "/main/{}/scores".format(label),
            where="columns in ['slope', 'intercept', 'SE_pctile'] & index!='{}' & index!='{}'".format(
                WILD_TYPE_VARIANT, SYNONYMOUS_VARIANT
            ),
        )
        se_data.sort_values("SE_pctile", inplace=True)
        se_data.dropna(axis="index", how="any", inplace=True)

        indices = np.linspace(0, len(se_data.index) - 1, 21).astype("int")
        se_data = se_data.ix[indices]

        # retrieves the whole DF because one case was hanging when trying to use select
        # totally unexplained, should fix later
        ratio_data = self.store.select("/main/{}/log_ratios".format(label)).loc[
            se_data.index
        ]

        fig, axarr = plt.subplots(7, 3, sharex="all", sharey="all")
        fig.subplots_adjust(
            hspace=0, wspace=0
        )  # eliminate white space between the subplots
        fig.set_size_inches((10, 17))
        fig.suptitle(
            "Representative {} Fits\n{} ({})".format(
                self.scoring_method, self.name, label.title()
            )
        )
        fig.subplots_adjust(top=0.958)  # eliminate white space after the title

        # turn off all tick labels
        plt.setp([ax.get_xticklabels() for ax in axarr.reshape(-1)], visible=False)
        plt.setp([ax.get_yticklabels() for ax in axarr.reshape(-1)], visible=False)

        # tick labels back on for some plots
        plt.setp(
            [ax.get_xticklabels() for ax in (axarr[-1, 0], axarr[-1, 2])], visible=True
        )
        plt.setp(
            [
                ax.get_yticklabels()
                for ax in (axarr[0, 0], axarr[2, 0], axarr[4, 0], axarr[6, 0])
            ],
            visible=True,
        )
        plt.setp(
            [ax.get_yticklabels() for ax in (axarr[1, -1], axarr[3, -1], axarr[5, -1])],
            visible=True,
        )
        plt.setp(
            [ax.yaxis for ax in (axarr[1, -1], axarr[3, -1], axarr[5, -1])],
            ticks_position="right",
        )

        # create the fit plots and add text to the individual plots
        for i, ax in enumerate(axarr.reshape(-1)):
            index = se_data.index[i]
            fit_axes(
                ax,
                xvalues=[x / float(max(self.timepoints)) for x in self.timepoints],
                yvalues=ratio_data.loc[index],
                slope=se_data.loc[index, "slope"],
                intercept=se_data.loc[index, "intercept"],
                xlabels=self.timepoints,
            )
            fit_axes_text(
                ax,
                cornertext="Slope {:3.2f}\nSE Pctile {:.1f}".format(
                    se_data.loc[index, "slope"], se_data.loc[index, "SE_pctile"]
                ),
                centertext=index,
            )

        # turn off the subplot axis labels
        [ax.set_xlabel("") for ax in axarr.reshape(-1)]
        [ax.set_ylabel("") for ax in axarr.reshape(-1)]

        # add x and y labels in the middle
        axarr[-1, 1].set_xlabel("Time Point")
        axarr[3, 0].set_ylabel("Log Ratio")

        pdf.savefig(fig, bbox_inches="tight")
        plt.close(fig)

    def timepoint_counts_plot(self, label, pdf):
        """
        Create barplot of the number of items counted for each time point.

        *label* is the data label (barcode, variant, etc.)

        *pdf* is an open PdfPages instance.
        """
        self.logger.info("Creating time point count plots ({})".format(label))

        counts = self.store["/main/{}/counts".format(label)].sum(axis="index")

        fig, ax = plt.subplots()
        fig.set_tight_layout(True)
        configure_axes(ax)

        xpos = np.arange(len(self.timepoints))
        width = 0.8
        ax.bar(xpos, counts, width, color=plot_colors["bright1"])
        ax.set_title("Total {}\n{}".format(label.title(), self.name))
        ax.set_ylabel("Count")
        ax.set_xlabel("Timepoint")
        ax.set_xticks(xpos + width / 2.0)
        ax.set_xticklabels(self.timepoints)

        pdf.savefig(fig)
        plt.close(fig)

    def volcano_plot(self, label, pdf, colors="YlGnBu_r", log_bins=True):
        """
        Create a volcano plot (p-value vs. functional score).

        *label* is the data label (barcode, variant, etc.)

        *pdf* is an open PdfPages instance.

        The p-values used are the regression p-values (p-value of non-zero slope). Due to the large number of points, we use a hexbin plot showing the density instead of a scatter plot.
        """
        self.logger.info("Creating volcano plot ({})".format(label))

        # get the data
        data = self.store.select(
            "/main/{}/scores".format(label), "columns=['score', 'pvalue_raw']"
        )
        volcano_plot(
            data,
            pdf,
            title="{} ({})".format(self.name, label.title()),
            colors=colors,
            log_bins=log_bins,
        )

    def make_plots(self):
        """
        Create plots for this entity.

        This function handles opening and closing the various PDF files for multi-page plots, as well as plotting similar plots for different data labels.
        """
        if self.plots_requested:
            self.logger.info("Creating plots")

            # counts per time point
            pdf = PdfPages(os.path.join(self.plot_dir, "timepoint_counts.pdf"))
            for label in self.labels:
                self.timepoint_counts_plot(label, pdf)
            pdf.close()

            # wild type shape
            if self.logr_method == "wt" and self.scoring_method in ("WLS", "OLS"):
                pdf = PdfPages(os.path.join(self.plot_dir, "wt_shape.pdf"))
                self.wt_plot(pdf)
                pdf.close()

            # regression weights
            if self.scoring_method == "WLS":
                pdf = PdfPages(os.path.join(self.plot_dir, "regression_weights.pdf"))
                for label in self.labels:
                    weights_plot(self, label, pdf)
                pdf.close()

            # linear fits by standard error percentile
            if self.scoring_method in ("WLS", "OLS"):
                pdf = PdfPages(os.path.join(self.plot_dir, "se_pctile.pdf"))
                for label in self.labels:
                    self.se_pctile_plot(label, pdf)
                pdf.close()

            # volcano plots
            # if self.scoring_method in ("WLS", "OLS", "ratios"):
            if self.scoring_method in ("WLS", "OLS") and "variants" in self.labels:
                pdf = PdfPages(os.path.join(self.plot_dir, "volcano.pdf"))
                for label in self.labels:
                    self.volcano_plot(label, pdf, log_bins=True)
                pdf.close()

            # library diversity for each time point (amino acid)
            if "synonymous" in self.labels:
                pdf = PdfPages(os.path.join(self.plot_dir, "diversity_aa.pdf"))
                for tp in self.timepoints:
                    self.sfmap_wrapper(
                        cname="c_{}".format(tp), pdf=pdf, coding=True, log10=True
                    )
                pdf.close()

            # library diversity for each time point (nucleotide)
            if "variants" in self.labels:
                pdf = PdfPages(os.path.join(self.plot_dir, "diversity_nt.pdf"))
                for tp in self.timepoints:
                    self.sfmap_wrapper(
                        cname="c_{}".format(tp), pdf=pdf, coding=False, log10=True
                    )
                pdf.close()

            # sequence-function maps
            if self.scoring_method != "counts":
                if "synonymous" in self.labels:
                    pdf = PdfPages(
                        os.path.join(self.plot_dir, "sequence_function_map_aa.pdf")
                    )
                    self.sfmap_wrapper(cname="score", pdf=pdf, coding=True)
                    pdf.close()
                if "variants" in self.labels:
                    pdf = PdfPages(
                        os.path.join(self.plot_dir, "sequence_function_map_nt.pdf")
                    )
                    self.sfmap_wrapper(cname="score", pdf=pdf, coding=False)
                    pdf.close()

        # SeqLib plots
        for lib in self.children:
            lib.make_plots()

    def write_tsv(self):
        """
        Write each table from the store to its own tab-separated file.

        Files are written to a ``tsv`` directory in the default output location. 
        File names are the HDF5 key with ``'_'`` substituted for ``'/'``.
        """
        if self.tsv_requested:
            self.logger.info("Generating tab-separated output files")
            for k in self.store.keys():
                self.write_table_tsv(k)
        for lib in self.children:
            lib.write_tsv()

    def synonymous_variants(self):
        """
        Populate and return a dictionary mapping synonymous variants to the 
        list of associated variants in ``/main/variants/counts``.
        """
        mapping = dict()
        try:
            variants = self.store.select_column("/main/variants/counts", "index")
        except KeyError:
            raise KeyError("No variant counts found [{}]".format(self.name))
        for v in variants:
            pv = protein_variant(v)
            try:
                mapping[pv].append(v)
            except KeyError:
                mapping[pv] = [v]
        return mapping

    def sfmap_wrapper(self, cname, pdf, coding, log10=False):
        """
        Create a sequence function map for either scores or library diversity.

        Uses :py:func:`~sfmap.sfmap_plot` for the plotting.
        """
        plot_options = self.get_root().plot_options

        if cname.startswith("c_"):
            counts = True
        elif cname == "score":
            counts = False
        else:
            raise ValueError("Invalid sequence-function map data column.")

        if coding:
            label = "amino acid"
        else:
            label = "nucleotide"

        if counts:
            self.logger.info("Creating diversity map ({})".format(label))
        else:
            self.logger.info("Creating sequence-function map ({})".format(label))

        # build the data frame name and get the data
        df_name = "/main/"
        if coding:
            df_name += "synonymous/"
        else:
            df_name += "variants/"
        if counts:
            df_name += "counts_unfiltered"
        else:
            df_name += "scores"
        if plot_options is not None:
            data, wtseq = singleton_dataframe(
                self.store[df_name][cname],
                self.wt,
                coding=coding,
                aa_list=plot_options["aa_list"],
            )
        else:
            data, wtseq = singleton_dataframe(
                self.store[df_name][cname], self.wt, coding=coding
            )
        if counts:
            data_se = None
        else:
            if plot_options is not None:
                data_se, _ = singleton_dataframe(
                    self.store[df_name]["SE"],
                    self.wt,
                    coding=coding,
                    aa_list=plot_options["aa_list"],
                )
            else:
                data_se, _ = singleton_dataframe(
                    self.store[df_name]["SE"], self.wt, coding=coding
                )

        # format the title
        if coding:
            title = "Amino Acid"
        else:
            title = "Nucleotide"
        if counts:
            title += " Diversity Map\n{} (Time {})".format(
                self.name, cname[2:]
            )  # trim the "c_"
        else:
            if self.scoring_method in ("WLS", "OLS"):
                title += " Sequence-Function Map\n{} ({} Slope)".format(
                    self.name, self.scoring_method
                )
            elif self.scoring_method == "ratios":
                title += " Sequence-Function Map\n{} ({})".format(
                    self.name, "Enrich2 Ratio"
                )
            elif self.scoring_method == "simple":
                title += " Sequence-Function Map\n{} ({})".format(
                    self.name, "Simplified Ratio"
                )
            else:
                raise ValueError("Invalid scoring method", self.name)

        if counts and log10:
            style = "logcounts"
        elif counts:
            style = "counts"
        else:
            style = "scores"

        if plot_options is not None:
            sfmap_plot(
                df=data,
                pdf=pdf,
                style=style,
                df_se=data_se,
                dimensions="tall",
                wt=wtseq,
                title=title,
                aa_list=plot_options["aa_list"],
                aa_label_groups=plot_options["aa_label_groups"],
            )
        else:
            sfmap_plot(
                df=data,
                pdf=pdf,
                style=style,
                df_se=data_se,
                dimensions="tall",
                wt=wtseq,
                title=title,
            )

    def barcodemap_mapping(self):
        mapping = dict()
        for bc, v in self.store["/main/barcodemap"].iterrows():
            v = v["value"]
            try:
                mapping[v].update([bc])
            except KeyError:
                mapping[v] = {bc}
        return mapping

    def calc_outliers(self, label, minimum_components=4, log_chunksize=20000):
        """
        Test whether an element's individual components have significantly different 
        scores from the element. Results are stored in ``'/main/<label>/outliers'``.

        Args:
            label (str): label for the component (``'variants'`` or ``'barcodes'``)

            minimum_components (int): minimum number of componenents required for any statistics to be calculated

            log_chunksize (int): will output a log message every *n* rows

        """
        if self.check_store("/main/{}/outliers".format(label)):
            return

        if label == "variants":
            label2 = "synonymous"
        elif label == "barcodes":
            if self.is_barcodevariant():
                label2 = "variants"
            elif self.is_barcodeid():
                label2 = "identifiers"
            else:
                # this should never happen
                raise AttributeError(
                    "Failed to identify parent label when calculating barcode outliers [{}]".format(
                        self.name
                    )
                )
        else:
            raise KeyError(
                "Invalid label '{}' for calc_outliers [{}]".format(label, self.name)
            )

        self.logger.info("Identifying outliers ({}-{})".format(label, label2))

        self.logger.info("Mapping {} to {}".format(label, label2))
        if label == "variants":
            mapping = self.synonymous_variants()
        elif label == "barcodes":
            mapping = self.barcodemap_mapping()
        else:
            raise KeyError(
                "Invalid label '{}' for calc_outliers [{}]".format(label, self.name)
            )

        # get the scores
        df1 = self.store.select(
            "/main/{}/scores".format(label), "columns=['score', 'SE']"
        )
        df2 = self.store.select(
            "/main/{}/scores".format(label2), "columns=['score', 'SE']"
        )

        # pre-calculate variances
        df1["var"] = df1["SE"] ** 2
        df1.drop("SE", axis=1, inplace=True)
        df2["var"] = df2["SE"] ** 2
        df2.drop("SE", axis=1, inplace=True)

        # set up empty results DF
        result_df = pd.DataFrame(
            np.nan, index=df1.index, columns=["z", "pvalue_raw", "parent"]
        )

        # because this step can be slow, output chunk-style logging messages
        # pre-calculate the lengths for the log messages
        log_chunksize_list = [log_chunksize] * (len(df2) / log_chunksize) + [
            len(df2) % log_chunksize
        ]
        log_chunk = 1

        for i, x in enumerate(df2.index):
            if i % log_chunksize == 0:
                self.logger.info(
                    "Calculating outlier p-values for chunk {} ({} rows) ({}-{})".format(
                        log_chunk, log_chunksize_list[log_chunk - 1], label, label2
                    )
                )
                log_chunk += 1
            try:
                components = df1.loc[mapping[x]].dropna(axis="index", how="all")
            except KeyError:
                # none of the components were in the index
                continue
            if len(components.index) >= minimum_components:
                for c in components.index:
                    zvalue = np.absolute(
                        df2.loc[x, "score"] - df1.loc[c, "score"]
                    ) / np.sqrt(df2.loc[x, "var"] + df1.loc[c, "var"])
                    result_df.loc[c, "z"] = zvalue
                    result_df.loc[c, "pvalue_raw"] = 2 * stats.norm.sf(zvalue)
                    result_df.loc[c, "parent"] = x
        if WILD_TYPE_VARIANT in result_df.index:
            result_df.loc[WILD_TYPE_VARIANT, "z"] = np.nan
            result_df.loc[WILD_TYPE_VARIANT, "pvalue_raw"] = np.nan
        result_df["z"] = result_df["z"].astype(float)
        result_df["pvalue_raw"] = result_df["pvalue_raw"].astype(float)

        self.store.put(
            "/main/{}/outliers".format(label),
            result_df,
            format="table",
            data_columns=result_df.columns,
        ) class Experiment(StoreManager):
    """
    Class for a coordinating multiple :py:class:`~.selection.Selection`
    objects. Creating an
    :py:class:`~experiment.Experiment` requires a valid *config* object,
    usually from a ``.json`` configuration file.
    """

    store_suffix = "exp"
    treeview_class_name = "Experiment"

    def __init__(self):
        StoreManager.__init__(self)
        self.conditions = list()
        self._wt = None
        self.logger = logging.getLogger("{}.{}".format(__name__, self.__class__))

    @property
    def wt(self):
        if self.has_wt_sequence():
            if self._wt is None:
                self._wt = self.selection_list()[0].wt.duplicate(self.name)
            return self._wt
        else:
            if self._wt is not None:
                raise ValueError(
                    "Experiment should not contain wild type "
                    "sequence [{}]".format(self.name)
                )
            else:
                return None

    def configure(self, cfg, configure_children=True):
        """
        Set up the :py:class:`~experiment.Experiment` using the *cfg* object,
        usually from a ``.json`` configuration file.
        """
        StoreManager.configure(self, cfg)
        self.logger = logging.getLogger(
            "{}.{} - {}".format(__name__, self.__class__.__name__, self.name)
        )
        if configure_children:
            if "conditions" not in cfg:
                raise KeyError(
                    "Missing required config value {} [{}]"
                    "".format("conditions", self.name)
                )
            for cnd_cfg in cfg["conditions"]:
                cnd = Condition()
                cnd.configure(cnd_cfg)
                self.add_child(cnd)

        selection_names = [x.name for x in self.selection_list()]
        if len(set(selection_names)) != len(selection_names):
            raise ValueError("Non-unique selection names [{}]" "".format(self.name))

    def serialize(self):
        """
        Format this object (and its children) as a config object suitable for
        dumping to a config file.
        """
        cfg = StoreManager.serialize(self)
        cfg["conditions"] = [child.serialize() for child in self.children]
        return cfg

    def _children(self):
        """
        Method bound to the ``children`` property. Returns a list of all
        :py:class:`~condition.Condition` objects belonging to this object,
        sorted by name.
        """
        return sorted(self.conditions, key=lambda x: x.name)

    def add_child(self, child):
        """
        Add a selection.
        """
        if child.name in self.child_names():
            raise ValueError(
                "Non-unique condition name '{}' [{}]" "".format(child.name, self.name)
            )
        child.parent = self
        self.conditions.append(child)

    def remove_child_id(self, tree_id):
        """
        Remove the reference to a :py:class:`~condition.Condition` with
        Treeview id *tree_id*.
        """
        self.conditions = [x for x in self.conditions if x.treeview_id != tree_id]

    def selection_list(self):
        """
        Return the :py:class:`~selection.Selection` objects as a list.
        """
        selections = list()
        for cnd in self.children:
            selections.extend(cnd.children)
        return selections

    def validate(self):
        """
        Calls validate on all child Conditions. Also checks the wild type
        sequence status.
        """
        # check the wild type sequences
        if self.has_wt_sequence():
            for child in self.selection_list()[1:]:
                if self.selection_list()[0].wt != child.wt:
                    self.logger.warning("Inconsistent wild type sequences")
                    break

        for child in self.children:
            child.validate()

    def is_coding(self):
        """
        Return ``True`` if the all :py:class:`~selection.Selection` in the 
        :py:class:`~experiment.Experiment` count protein-coding variants, else 
        ``False``.
        """
        return all(x.is_coding() for x in self.selection_list())

    def has_wt_sequence(self):
        """
        Return ``True`` if the all :py:class:`~selection.Selection` in the 
        :py:class:`~experiment.Experiment` have a wild type sequence, else 
        ``False``.
        """
        return all(x.has_wt_sequence() for x in self.selection_list())

    def calculate(self):
        """
        Calculate scores for all :py:class:`~selection.Selection` objects.
        """
        if len(self.labels) == 0:
            raise ValueError(
                "No data present across all conditions [{}]" "".format(self.name)
            )
        for s in self.selection_list():
            s.calculate()
        self.combine_barcode_maps()
        for label in self.labels:
            self.calc_counts(label)
            if self.scoring_method != "counts":
                self.calc_shared_full(label)
                self.calc_shared(label)
                self.calc_scores(label)
                if label != "barcodes":
                    self.calc_pvalues_wt(label)

    def combine_barcode_maps(self):
        """
        Combine all barcode maps for :py:class:`~selection.Selection` objects
        into a single data frame and store it in ``'/main/barcodemap'``.

        If multiple variants or IDs map to the same barcode, only the first one
        will be present in the barcode map table.

        The ``'/main/barcodemap'`` table is not created if no 
        :py:class:`~selection.Selection` has barcode map information.
        """
        if self.check_store("/main/barcodemap"):
            return

        bcm = None
        for sel in self.selection_list():
            if "/main/barcodemap" in sel.store.keys():
                if bcm is None:
                    bcm = sel.store["/main/barcodemap"]
                else:
                    bcm = bcm.join(
                        sel.store["/main/barcodemap"], rsuffix=".drop", how="outer"
                    )
                    new = bcm.loc[pd.isnull(bcm)["value"]]
                    bcm.loc[new.index, "value"] = new["value.drop"]
                    bcm.drop("value.drop", axis="columns", inplace=True)
        if bcm is not None:
            bcm.sort_values("value", inplace=True)
            self.store.put(
                "/main/barcodemap", bcm, format="table", data_columns=bcm.columns
            )

    def calc_counts(self, label):
        """
        Create a data frame of all counts in this Experiment. This data frame
        is not used for any calculations, but is provided to facilitate
        exploration of the data set.
        """
        if self.check_store("/main/{}/counts".format(label)):
            return

        # create columns multi-index
        # has to be lex-sorted for multi-slicing to work
        self.logger.info("Creating column multi-index for counts ({})".format(label))
        conditions_index = list()
        selections_index = list()
        values_index = list()

        for cnd in self.children:
            for sel in cnd.children:
                conditions_index.extend([cnd.name] * len(sel.timepoints))
                selections_index.extend([sel.name] * len(sel.timepoints))
                values_index.extend(["c_{}".format(x) for x in sorted(sel.timepoints)])
        columns = pd.MultiIndex.from_tuples(
            zip(conditions_index, selections_index, values_index),
            names=["condition", "selection", "timepoint"],
        )

        # create union index
        self.logger.info("Creating row index for counts ({})".format(label))
        combined = None
        first = True
        for s in self.selection_list():
            if first:
                combined = s.store.select(
                    "/main/{}/counts_unfiltered" "".format(label), "columns='index'"
                ).index
                first = False
            else:
                combined = combined.join(
                    s.store.select(
                        "/main/{}/counts_unfiltered".format(label), "columns='index'"
                    ).index,
                    how="outer",
                )

        # create and fill the data frames
        self.logger.info(
            "Populating Experiment data frame with counts ({})".format(label)
        )
        data = pd.DataFrame(index=combined, columns=columns)
        for cnd in self.children:
            for sel in cnd.children:
                sel_data = sel.store.select(
                    "/main/{}/counts_unfiltered" "".format(label)
                )
                for tp in sel.timepoints:
                    data.loc[:][cnd.name, sel.name, "c_{}".format(tp)] = sel_data[
                        "c_{}".format(tp)
                    ]
        self.store.put("/main/{}/counts".format(label), data, format="table")

    def calc_shared_full(self, label):
        """
        Use joins to create a data frame containing all scores across all
        Selections in the Experiment.
        """
        if self.check_store("/main/{}/scores_shared_full".format(label)):
            return

        # create columns multi-index
        # has to be lex-sorted for multi-slicing to work
        self.logger.info("Creating column multi-index for scores ({})")
        conditions_index = list()
        selections_index = list()
        values_index = list()

        if self.scoring_method == "simple":
            values_list = ["score"]
        else:
            values_list = ["score", "SE"]

        for cnd in self.children:
            for sel in cnd.children:
                conditions_index.extend([cnd.name] * len(values_list))
                selections_index.extend([sel.name] * len(values_list))
                values_index.extend(sorted(values_list))
        columns = pd.MultiIndex.from_tuples(
            zip(conditions_index, selections_index, values_index),
            names=["condition", "selection", "value"],
        )

        # create union index
        self.logger.info("Creating row index for scores ({})".format(label))
        combined = None
        first = True
        for s in self.selection_list():
            if first:
                combined = s.store.select(
                    "/main/{}/scores".format(label), "columns='index'"
                ).index
                first = False
            else:
                combined = combined.join(
                    s.store.select(
                        "/main/{}/scores".format(label), "columns='index'"
                    ).index,
                    how="outer",
                )

        # create and fill the data frames
        self.logger.info(
            "Populating Experiment data frame with scores ({})".format(label)
        )
        data = pd.DataFrame(index=combined, columns=columns)
        for cnd in self.children:
            for sel in cnd.children:
                sel_data = sel.store.select("/main/{}/scores".format(label))
                for v in values_list:
                    data.loc[:, (cnd.name, sel.name, v)] = sel_data[v]
        self.store.put(
            "/main/{}/scores_shared_full".format(label), data, format="table"
        )

    def calc_shared(self, label):
        """
        Get the subset of scores that are shared across all Selections in each
        Condition.
        """
        if self.check_store("/main/{}/scores_shared".format(label)):
            return

        idx = pd.IndexSlice

        self.logger.info(
            "Identifying subset shared across all Selections ({})".format(label)
        )
        data = self.store.select("/main/{}/scores_shared_full".format(label))

        # identify variants found in all selections in at least one condition
        complete = np.full(len(data.index), False, dtype=bool)
        for cnd in data.columns.levels[0]:
            complete = np.logical_or(
                complete, data.loc[:, idx[cnd, :, :]].notnull().all(axis="columns")
            )

        data = data.loc[complete]

        self.store.put("/main/{}/scores_shared".format(label), data, format="table")

    def calc_scores(self, label):
        """
        Combine the scores and standard errors within each condition.
        """
        if self.check_store("/main/{}/scores".format(label)):
            return

        self.logger.info("Calculating per-condition scores ({})".format(label))

        # set up new data frame
        shared_index = self.store.select(
            "/main/{}/scores_shared" "".format(label), columns="index"
        ).index
        columns = pd.MultiIndex.from_product(
            [sorted(self.child_names()), sorted(["score", "SE", "epsilon"])],
            names=["condition", "value"],
        )
        data = pd.DataFrame(np.nan, index=shared_index, columns=columns)
        del shared_index
        del columns

        # set up local variables
        idx = pd.IndexSlice

        score_df = self.store.select("/main/{}/scores_shared".format(label))
        if self.scoring_method == "simple":
            # special case for simple ratios that have no SE
            # calculates the average score
            for cnd in score_df.columns.levels[0]:
                data.loc[:, idx[cnd, "score"]] = score_df.loc[
                    :, idx[cnd, :, "score"]
                ].mean(axis=1)
        else:
            for cnd in score_df.columns.levels[0]:
                y = np.array(score_df.loc[:, idx[cnd, :, "score"]].values).T
                sigma2i = np.array(score_df.loc[:, idx[cnd, :, "SE"]].values ** 2).T

                # single replicate of the condition
                if y.shape[0] == 1:
                    data.loc[:, idx[cnd, "score"]] = y.ravel()
                    data.loc[:, idx[cnd, "SE"]] = np.sqrt(sigma2i).ravel()
                    data.loc[:, idx[cnd, "epsilon"]] = 0.0
                # multiple replicates
                else:
                    betaML, var_betaML, eps = rml_estimator(y, sigma2i)
                    data.loc[:, idx[cnd, "score"]] = betaML
                    data.loc[:, idx[cnd, "SE"]] = np.sqrt(var_betaML)
                    data.loc[:, idx[cnd, "epsilon"]] = eps

                # special case for normalized wild type variant
                if self.logr_method == "wt" and WILD_TYPE_VARIANT in data.index:
                    data.loc[WILD_TYPE_VARIANT, idx[:, "SE"]] = 0.0
                    data.loc[WILD_TYPE_VARIANT, idx[:, "score"]] = 0.0
                    data.loc[WILD_TYPE_VARIANT, idx[:, "epsilon"]] = 0.0

        # store the data
        self.store.put("/main/{}/scores".format(label), data, format="table")

    def calc_pvalues_wt(self, label):
        """
        Calculate uncorrected pvalue for each variant compared to wild type.
        """
        if self.check_store("/main/{}/scores_pvalues_wt".format(label)):
            return

        idx = pd.IndexSlice

        wt = self.store.select(
            "/main/{}/scores".format(label), "index=WILD_TYPE_VARIANT"
        )
        if len(wt) == 0:  # no wild type score
            self.logger.info(
                "Failed to find wild type score, skipping wild type p-value calculations"
            )
            return
        data = self.store.select(
            "/main/{}/scores".format(label), "index!=WILD_TYPE_VARIANT"
        )

        columns = pd.MultiIndex.from_product(
            [sorted(self.child_names()), sorted(["z", "pvalue_raw"])],
            names=["condition", "value"],
        )
        result_df = pd.DataFrame(index=data.index, columns=columns)

        condition_labels = data.columns.levels[0]
        for cnd in condition_labels:
            result_df.loc[:, idx[cnd, "z"]] = np.absolute(
                wt.loc[WILD_TYPE_VARIANT, idx[cnd, "score"]]
                - data.loc[:, idx[cnd, "score"]]
            ) / np.sqrt(
                wt.loc[WILD_TYPE_VARIANT, idx[cnd, "SE"]] ** 2
                + data.loc[:, idx[cnd, "SE"]] ** 2
            )
            result_df.loc[:, idx[cnd, "pvalue_raw"]] = 2 * stats.norm.sf(
                result_df.loc[:, idx[cnd, "z"]]
            )

        self.store.put(
            "/main/{}/scores_pvalues_wt".format(label), result_df, format="table"
        )

    def calc_pvalues_pairwise(self, label):
        """
        Calculate pvalues for each variant in each pair of Conditions.
        """
        if self.check_store("/main/{}/scores_pvalues".format(label)):
            return

        data = self.store["/main/{}/scores".format(label)]

        cnd1_index = list()
        cnd2_index = list()
        values_index = list()

        values_list = ["z", "pvalue_raw"]

        condition_labels = data.columns.levels[0]
        for i, cnd1 in enumerate(condition_labels):
            for cnd2 in condition_labels[i + 1 :]:
                cnd1_index.extend([cnd1] * len(values_list))
                cnd2_index.extend([cnd2] * len(values_list))
                values_index.extend(sorted(values_list))
        columns = pd.MultiIndex.from_tuples(
            zip(cnd1_index, cnd2_index, values_index),
            names=["condition1", "condition2", "value"],
        )

        idx = pd.IndexSlice
        result_df = pd.DataFrame(np.nan, index=data.index, columns=columns)
        for i, cnd1 in enumerate(condition_labels):
            for cnd2 in condition_labels[i + 1 :]:
                result_df.loc[:, idx[cnd1, cnd2, "z"]] = np.absolute(
                    data.loc[:, idx[cnd1, "score"]] - data.loc[:, idx[cnd2, "score"]]
                ) / np.sqrt(
                    data.loc[:, idx[cnd1, "SE"]] ** 2
                    + data.loc[:, idx[cnd2, "SE"]] ** 2
                )
                result_df.loc[:, idx[cnd1, cnd2, "pvalue_raw"]] = 2 * stats.norm.sf(
                    result_df.loc[:, idx[cnd1, cnd2, "z"]]
                )

        self.store.put(
            "/main/{}/scores_pvalues".format(label), result_df, format="table"
        )

    def make_plots(self):
        if self.plots_requested:
            self.logger.info("Creating plots")

            # sequence-function maps
            if self.scoring_method != "counts":
                if "synonymous" in self.labels:
                    pdf = PdfPages(
                        os.path.join(self.plot_dir, "sequence_function_map_aa.pdf")
                    )
                    for condition in self.children:
                        self.sfmap_wrapper(
                            condition=condition.name, pdf=pdf, coding=True
                        )
                    pdf.close()
                if "variants" in self.labels:
                    pdf = PdfPages(
                        os.path.join(self.plot_dir, "sequence_function_map_nt.pdf")
                    )
                    for condition in self.children:
                        self.sfmap_wrapper(
                            condition=condition.name, pdf=pdf, coding=False
                        )
                    pdf.close()

        for s in self.selection_list():
            s.make_plots()

    def write_tsv(self):
        """
        Write each table from the store to its own tab-separated file.

        Files are written to a ``tsv`` directory in the default output
        location. File names are the HDF5 key with ``'_'`` substituted for
        ``'/'``.
        """
        if self.tsv_requested:
            self.logger.info("Generating tab-separated output files")
            for k in self.store.keys():
                self.write_table_tsv(k)
        for s in self.selection_list():
            s.write_tsv()

    def sfmap_wrapper(self, condition, pdf, coding):
        """
        Create a sequence function map for scores in *condition*.

        Uses :py:func:`~sfmap.sfmap_plot` for the plotting.
        """
        plot_options = self.get_root().plot_options

        if coding:
            label = "amino acid"
        else:
            label = "nucleotide"

        self.logger.info(
            "Creating sequence-function map ({}, {})".format(condition, label)
        )

        idx = pd.IndexSlice

        if coding:
            df_name = "/main/synonymous/scores"
        else:
            df_name = "/main/variants/scores"

        if plot_options is not None:
            data, wtseq = singleton_dataframe(
                self.store[df_name][idx[condition, "score"]],
                self.wt,
                coding=coding,
                aa_list=plot_options["aa_list"],
            )
            data_se, _ = singleton_dataframe(
                self.store[df_name][idx[condition, "SE"]],
                self.wt,
                coding=coding,
                aa_list=plot_options["aa_list"],
            )
        else:
            data, wtseq = singleton_dataframe(
                self.store[df_name][idx[condition, "score"]], self.wt, coding=coding
            )
            data_se, _ = singleton_dataframe(
                self.store[df_name][idx[condition, "SE"]], self.wt, coding=coding
            )

        # format the title
        if coding:
            title = "Amino Acid"
        else:
            title = "Nucleotide"

        if self.scoring_method in ("WLS", "OLS"):
            title += " Sequence-Function Map\n{} ({} Slope)".format(
                condition, self.scoring_method
            )
        elif self.scoring_method == "ratios":
            title += " Sequence-Function Map\n{} ({})".format(
                condition, "Enrich2 Ratio"
            )
        elif self.scoring_method == "simple":
            title += " Sequence-Function Map\n{} ({})".format(
                condition, "Simplified Ratio"
            )
        else:
            raise ValueError("Invalid scoring method", self.name)

        if plot_options is not None:
            sfmap_plot(
                df=data,
                pdf=pdf,
                style="scores",
                df_se=data_se,
                dimensions="tall",
                wt=wtseq,
                title=title,
                aa_list=plot_options["aa_list"],
                aa_label_groups=plot_options["aa_label_groups"],
            )
        else:
            sfmap_plot(
                df=data,
                pdf=pdf,
                style="scores",
                df_se=data_se,
                dimensions="tall",
                wt=wtseq,
                title=title,
            )

    def correlation_plot(self, pdf, label):
        """
        Create a triangular heatmap showing the Pearson correlation coefficient
        for each pairwise comparison of replicate scores.
        """
        pass
import Tkinter as tk
import ttk
import tkSimpleDialog
from .dialog_elements import FileEntry, StringEntry, DEFAULT_COLUMNS
from .create_seqlib_dialog import SEQLIB_LABEL_TEXT
from ..barcode import BarcodeSeqLib
from ..barcodevariant import BcvSeqLib
from ..barcodeid import BcidSeqLib
from ..basic import BasicSeqLib
from ..idonly import IdOnlySeqLib
from ..overlap import OverlapSeqLib
from ..selection import Selection
from ..experiment import Experiment
from __future__ import print_function


#: map class names to class definitions to avoid use of globals()
ELEMENT_CLASSES = {
    "BarcodeSeqLib": BarcodeSeqLib,
    "BcvSeqLib": BcvSeqLib,
    "BcidSeqLib": BcidSeqLib,
    "BasicSeqLib": BasicSeqLib,
    "IdOnlySeqLib": IdOnlySeqLib,
    "OverlapSeqLib": OverlapSeqLib,
    "Selection": Selection,
    "Experiment": Experiment,
}


class CreateRootDialog(tkSimpleDialog.Dialog):
    """
    Dialog box for creating a new root element.
    """

    def __init__(self, parent_window, title="Create Root Object"):
        self.element_tkstring = tk.StringVar()
        self.cfg_dict = dict()
        self.output_directory_tk = FileEntry(
            "Output Directory",
            self.cfg_dict,
            "output directory",
            optional=False,
            directory=True,
        )

        self.name_tk = StringEntry("Name", self.cfg_dict, "name", optional=False)