#!/usr/bin/env python2
# accession_analysis 0.0.1
# Generated by dx-app-wizard.
#
# Basic execution pattern: Your app will run on a single machine from
# beginning to end.
#
# See https://wiki.dnanexus.com/Developer-Portal for documentation and
# tutorials on how to modify this file.
#
# DNAnexus Python Bindings (dxpy) documentation:
#   http://autodoc.dnanexus.com/bindings/python/current/

import os
import sys
import subprocess
import logging
import traceback
import re
import urlparse
import time
import pprint
import csv
import json
import copy
from base64 import b64encode
import tempfile

import dxpy
import common

logging.basicConfig()
# logging.getLogger("requests").setLevel(logging.WARNING)
logger = logging.getLogger(__name__)
logger.addHandler(dxpy.DXLogHandler())
logger.propagate = False
logger.setLevel(logging.INFO)
logger.info('Logging from the applet is activated')

DCC_CREDENTIALS_PROJECT = 'project-F30FzF0048K9JZKxPvB3Y563'
DCC_CREDENTIALS_FOLDER = '/credentials'

COMMON_METADATA = {
    'lab': 'encode-processing-pipeline',
    'award': 'U41HG006992'
}

DEPRECATED = ['deleted', 'replaced', 'revoked']

# This is a mapping from pipeline version to step version names and
# ENCODEd analysis_step_version aliases
# Version dates are only used to infer the pipeline version if it cannot be
# inferred directly from the workflow (for instance if the workflow has been
# deleted)
# As seconds from the epoch
VERSION_TIMES = {
    '1': 0,
    '1.2': 1471460400  # "Wed Aug 17 12:00:00 2016"
}
STEP_VERSION_ALIASES = {
    'default': {
        'bwa-alignment-step':                      'ced7264f-9206-4838-99ea-3badb3fb124c',
        'bwa-indexing-step':                       '58894736-6c6b-430a-ac19-4f86cede236a',
        'bwa-raw-alignment-step':                  'b9e5856e-b6bb-4e5c-8564-3e6c79ec9b14',
        'histone-peak-calling-step':               'ba41eaf0-db75-44b6-8045-11c0d84d2a5b',
        'histone-overlap-peaks-step':              '5907a64c-1a1c-4662-a429-5aa6c1f77bea',
        'histone-peaks-to-bigbed-step':            'a7d5cbf0-c9f8-446a-9db8-a50b443e6628',
        'histone-replicated-peaks-to-bigbed-step': 'b67cd935-d665-4ba7-bf34-17273db5451e',
        'tf-macs2-signal-calling-step':            'd2619376-aeb4-4fe4-9374-785133195864',
        'tf-spp-peak-calling-step':                '8250c103-b2a7-4ac7-842c-65e2ef651fd6',
        'tf-idr-step':                             'ad55c9a7-76f4-483d-9feb-ff47c067c523',
        'tf-peaks-to-bigbed-step':                 '511ba3d7-4b15-49a8-901e-c7516f3312c5',
        'tf-idr-peaks-to-bigbed-step':             '30b7cd3c-55fd-42c2-8175-19bc98d58752'
    },
    '1': {
        'bwa-alignment-step':                      'ced7264f-9206-4838-99ea-3badb3fb124c',
        'bwa-indexing-step':                       '58894736-6c6b-430a-ac19-4f86cede236a',
        'bwa-raw-alignment-step':                  'b9e5856e-b6bb-4e5c-8564-3e6c79ec9b14',
        'histone-peak-calling-step':               'ba41eaf0-db75-44b6-8045-11c0d84d2a5b',
        'histone-overlap-peaks-step':              '5907a64c-1a1c-4662-a429-5aa6c1f77bea',
        'histone-peaks-to-bigbed-step':            'a7d5cbf0-c9f8-446a-9db8-a50b443e6628',
        'histone-replicated-peaks-to-bigbed-step': 'b67cd935-d665-4ba7-bf34-17273db5451e',
        'tf-macs2-signal-calling-step':            'd2619376-aeb4-4fe4-9374-785133195864',
        'tf-spp-peak-calling-step':                '8250c103-b2a7-4ac7-842c-65e2ef651fd6',
        'tf-idr-step':                             'ad55c9a7-76f4-483d-9feb-ff47c067c523',
        'tf-peaks-to-bigbed-step':                 '511ba3d7-4b15-49a8-901e-c7516f3312c5',
        'tf-idr-peaks-to-bigbed-step':             '30b7cd3c-55fd-42c2-8175-19bc98d58752'
    },
    '1.2': {
        'bwa-alignment-step':                      'dc6dac54-0a5f-4356-815c-975c4b052249',
        'bwa-indexing-step':                       '58894736-6c6b-430a-ac19-4f86cede236a',
        'bwa-raw-alignment-step':                  '353a6f3b-999c-4f01-84ce-5de9e83dab22',

        'histone-peak-calling-step':               '6ef1d1f5-b075-4178-bad8-bc1b8319f680',
        'histone-overlap-peaks-step':              '06a53965-b1b1-4d2c-999f-189c10716096',
        'histone-peaks-to-bigbed-step':            'a7d5cbf0-c9f8-446a-9db8-a50b443e6628',
        'histone-replicated-peaks-to-bigbed-step': 'b67cd935-d665-4ba7-bf34-17273db5451e',

        'histone-unreplicated-peak-calling-step':                           '5f9d0122-9974-4750-a8ff-1dad98485c28',
        'histone-unreplicated-partition-concordance-step':                  'fb4a04ac-a13c-4a6f-86b0-f4781e7e2232',
        'histone-unreplicated-peaks-to-bigbed-step':                        '2b8daae2-cdf6-41b1-bcd5-6ce99b0970d2',
        'histone-unreplicated-partition-concordance-peaks-to-bigbed-step':  'db3b685c-8462-49b8-a29c-5f6b180a9294',

        'tf-macs2-signal-calling-step':            '37284851-f58c-4c4e-898d-0b595699cd20',
        'tf-spp-peak-calling-step':                '63326924-00a7-414f-ba7d-923c0f25da92',
        'tf-idr-step':                             'ad55c9a7-76f4-483d-9feb-ff47c067c523',
        'tf-peaks-to-bigbed-step':                 '511ba3d7-4b15-49a8-901e-c7516f3312c5',
        'tf-idr-peaks-to-bigbed-step':             '30b7cd3c-55fd-42c2-8175-19bc98d58752',

        'tf-unreplicated-macs2-signal-calling-step': 'e2902ced-8567-4b82-aa6a-99ff5bbcc916',
        'tf-unreplicated-spp-peak-calling-step':     '48d55aea-f922-4102-b93e-a57c4af26f6c',
        'tf-unreplicated-idr-step':                  'abd4401d-96aa-4c37-894d-8e0fee9273ec',
        'tf-unreplicated-peaks-to-bigbed-step':      'abcc5287-39ae-456a-b587-3af3fc202e00',
        'tf-unreplicated-idr-peaks-to-bigbed-step':  'fa044b11-91ff-4320-a212-5f98f9f7e165'

    }
}


class AccessioningError(Exception):
    def __init__(self, value):
        self.value = value

    def __str__(self):
        return repr(self.value)


def flat(l):
    result = []
    for el in l:
        if hasattr(el, "__iter__") and not isinstance(el, basestring):
            result.extend(flat(el))
        else:
            result.append(el)
    return result


def dup_parse(dxlink):
    desc = dxpy.describe(dxlink)
    with dxpy.DXFile(desc['id'], mode='r') as dup_file:
        if not dup_file:
            return None

        lines = iter(dup_file.read().splitlines())

        for line in lines:
            if line.startswith('## METRICS CLASS'):
                headers = lines.next().rstrip('\n').lower()
                metrics = lines.next().rstrip('\n')
                break

        headers = headers.split('\t')
        metrics = metrics.split('\t')
        headers.pop(0)
        metrics.pop(0)

        dup_qc = dict(zip(headers, metrics))
    return dup_qc


def xcor_parse(dxlink):
    desc = dxpy.describe(dxlink)
    with dxpy.DXFile(desc['id'], mode='r') as xcor_file:
        if not xcor_file:
            return None

        lines = xcor_file.read().splitlines()
        line = lines[0].rstrip('\n')
        # CC_SCORE FILE format:
        #   Filename <tab>
        #   numReads <tab>
        #   estFragLen <tab>
        #   corr_estFragLen <tab>
        #   PhantomPeak <tab>
        #   corr_phantomPeak <tab>
        #   argmin_corr <tab>
        #   min_corr <tab>
        #   phantomPeakCoef <tab>
        #   relPhantomPeakCoef <tab>
        #   QualityTag

        headers = ['Filename',
                   'numReads',
                   'estFragLen',
                   'corr_estFragLen',
                   'PhantomPeak',
                   'corr_phantomPeak',
                   'argmin_corr',
                   'min_corr',
                   'phantomPeakCoef',
                   'relPhantomPeakCoef',
                   'QualityTag']
        metrics = line.split('\t')
        headers.pop(0)
        metrics.pop(0)

        xcor_qc = dict(zip(headers, metrics))
    return xcor_qc


def pbc_parse(dxlink):
    desc = dxpy.describe(dxlink)
    with dxpy.DXFile(desc['id'], mode='r') as pbc_file:
        if not pbc_file:
            return None

        lines = pbc_file.read().splitlines()
        line = lines[0].rstrip('\n')
        # PBC File output:
        #   TotalReadPairs <tab>
        #   DistinctReadPairs <tab>
        #   OneReadPair <tab>
        #   TwoReadPairs <tab>
        #   NRF=Distinct/Total <tab>
        #   PBC1=OnePair/Distinct <tab>
        #   PBC2=OnePair/TwoPair

        headers = ['TotalReadPairs',
                   'DistinctReadPairs',
                   'OneReadPair',
                   'TwoReadPairs',
                   'NRF',
                   'PBC1',
                   'PBC2']
        metrics = line.split('\t')

        pbc_qc = dict(zip(headers, metrics))
    return pbc_qc


def flagstat_parse(dxlink):
    desc = dxpy.describe(dxlink)
    with dxpy.DXFile(desc['id'], mode='r') as flagstat_file:
        if not flagstat_file:
            return None
        flagstat_lines = flagstat_file.read().splitlines()

    qc_dict = {
        # values are regular expressions,
        # will be replaced with scores [hiq, lowq]
        'in_total': 'in total',
        'duplicates': 'duplicates',
        'mapped': 'mapped',
        'paired_in_sequencing': 'paired in sequencing',
        'read1': 'read1',
        'read2': 'read2',
        'properly_paired': 'properly paired',
        'with_self_mate_mapped': 'with itself and mate mapped',
        'singletons': 'singletons',
        # i.e. at the end of the line
        'mate_mapped_different_chr': 'with mate mapped to a different chr$',
        # RE so must escape
        'mate_mapped_different_chr_hiQ':
            'with mate mapped to a different chr \(mapQ>=5\)'
    }

    for (qc_key, qc_pattern) in qc_dict.items():
        qc_metrics = next(re.split(qc_pattern, line)
                          for line in flagstat_lines
                          if re.search(qc_pattern, line))
        (hiq, lowq) = qc_metrics[0].split(' + ')
        qc_dict[qc_key] = [int(hiq.rstrip()), int(lowq.rstrip())]

    return qc_dict


def get_attachment(dxlink):
    desc = dxpy.describe(dxlink)
    filename = desc['name']
    mime_type = desc['media']
    if mime_type == 'text/plain' and not filename.endswith(".txt"):
        filename += ".txt"
    with dxpy.DXFile(desc['id'], mode='r') as stream:
        obj = {
            'download': filename,
            'type': mime_type,
            'href': 'data:%s;base64,%s' % (mime_type, b64encode(stream.read()))
        }
    return obj


# these are stopgaps until proper QC metrics are available
# and mapping_report uses them.

def qc(stages):
    raw_mapping_stage = next(
        stages[stage_name]['stage_metadata']
        for stage_name in stages.keys()
        if stage_name.startswith("Map ENCSR"))
    return flagstat_parse(raw_mapping_stage['output']['mapping_statistics'])


def dup_qc(stages):
    qc_stage = next(
        stages[stage_name]['stage_metadata']
        for stage_name in stages.keys()
        if stage_name.startswith("Filter and QC"))
    return dup_parse(qc_stage['output']['dup_file_qc'])


def pbc_qc(stages):
    qc_stage = next(
        stages[stage_name]['stage_metadata']
        for stage_name in stages.keys()
        if stage_name.startswith("Filter and QC"))
    return pbc_parse(qc_stage['output']['pbc_file_qc'])


def filtered_qc(stages):
    qc_stage = next(
        stages[stage_name]['stage_metadata']
        for stage_name in stages.keys()
        if stage_name.startswith("Filter and QC"))
    return flagstat_parse(qc_stage['output']['filtered_mapstats'])


def xcor_qc(stages):
    xcor_stage = next(
        stages[stage_name]['stage_metadata']
        for stage_name in stages.keys()
        if stage_name.startswith("Calculate cross-correlation"))
    return xcor_parse(xcor_stage['output']['CC_scores_file'])


def chipseq_filter_quality_metric(step_run, stages, files):
    # this is currently a mix of deduplication and cross-correlation stats
    # maybe break out all the xcor stuff to its own object
    logger.debug(
        "in chip_seq_filter_quality_metric with "
        "step_run %s stages.keys() %s output files %s"
        % (step_run, stages.keys(), files))

    file_accessions = list(set(flat([
        resolve_name_to_accessions(stages, output_name)
        for output_name in files])))

    qc_stage = next(
        stages[stage_name]['stage_metadata']
        for stage_name in stages.keys()
        if stage_name.startswith("Filter and QC"))
    xcor_stage = next(
        stages[stage_name]['stage_metadata']
        for stage_name in stages.keys()
        if stage_name.startswith("Calculate cross-correlation"))

    xcor_plot = get_attachment(xcor_stage['output']['CC_plot_file'])
    xcor_scores = get_attachment(xcor_stage['output']['CC_scores_file'])

    pbc_qc = pbc_parse(qc_stage['output']['pbc_file_qc'])
    xcor_qc = xcor_parse(xcor_stage['output']['CC_scores_file'])

    obj = {
        # 'assay_term_id': 'OBI:0000716',
        'assay_term_name': 'ChIP-seq',
        'step_run': step_run,
        'quality_metric_of': file_accessions,
        'cross_correlation_plot': xcor_plot,
        'attachment': xcor_scores,
        'NSC': float(xcor_qc['phantomPeakCoef']),
        'RSC': float(xcor_qc['relPhantomPeakCoef']),
        'fragment length': int(xcor_qc['estFragLen']),
        'PBC1': float(pbc_qc['PBC1']),
        'PBC2': float(pbc_qc['PBC2']),
        'NRF': float(pbc_qc['NRF'])
    }
    # for very small, very complex fastq's, there will be exactly no overlap
    # in such cases, PBC2 will be infinite, but there is no infinity in JSON,
    # so we change that to the string "Infinity"
    if obj['PBC2'] == float('inf'):
        obj['PBC2'] = 'Infinity'
    obj.update(COMMON_METADATA)
    return [obj]


def get_flagstat_obj(step_run, stage, file_accessions):

    if 'filtered_mapstats' in stage['output']:
        flagstat_qc = flagstat_parse(stage['output']['filtered_mapstats'])
        processing_stage = 'filtered'
        attachment = get_attachment(stage['output']['filtered_mapstats'])
    elif 'mapping_statistics' in stage['output']:
        flagstat_qc = flagstat_parse(stage['output']['mapping_statistics'])
        processing_stage = 'unfiltered'
        attachment = get_attachment(stage['output']['mapping_statistics'])
    else:
        logger.error(
            'get_flagstat_obj: No filtered_mapstats or mapping_statistics')
        logger.debug(
            'get_flagstat_obj: Stage %s' % (stage.get('name')))
        return None

    obj = {
        # 'assay_term_id': 'OBI:0000716',
        'assay_term_name': 'ChIP-seq',
        'step_run': step_run,
        'quality_metric_of': file_accessions,
        'processing_stage': processing_stage,
        'attachment': attachment,
        'total':                int(flagstat_qc['in_total'][0]),
        'total_qc_failed':      int(flagstat_qc['in_total'][1]),
        'duplicates':           int(flagstat_qc['duplicates'][0]),
        'duplicates_qc_failed': int(flagstat_qc['duplicates'][1]),
        'mapped':               int(flagstat_qc['mapped'][0]),
        'mapped_qc_failed':     int(flagstat_qc['mapped'][1]),
        'mapped_pct':           '{:.2%}'.format(
            float(flagstat_qc['mapped'][0]) /
            float(flagstat_qc['in_total'][0]))
    }
    if int(flagstat_qc['paired_in_sequencing'][0]) or \
       int(flagstat_qc['paired_in_sequencing'][1]):
        obj.update({
            'paired':
                int(flagstat_qc['paired_in_sequencing'][0]),

            'paired_qc_failed':
                int(flagstat_qc['paired_in_sequencing'][1]),

            'read1':
                int(flagstat_qc['read1'][0]),

            'read1_qc_failed':
                int(flagstat_qc['read1'][1]),

            'read2':
                int(flagstat_qc['read2'][0]),

            'read2_qc_failed':
                int(flagstat_qc['read2'][1]),

            'paired_properly':
                int(flagstat_qc['properly_paired'][0]),

            'paired_properly_qc_failed':
                int(flagstat_qc['properly_paired'][1]),

            'paired_properly_pct':
                '{:.2%}'.format(
                    float(flagstat_qc['properly_paired'][0]) /
                    float(flagstat_qc['in_total'][0])),

            'with_itself':
                int(flagstat_qc['with_self_mate_mapped'][0]),

            'with_itself_qc_failed':
                int(flagstat_qc['with_self_mate_mapped'][1]),

            'singletons':
                int(flagstat_qc['singletons'][0]),

            'singletons_qc_failed':
                int(flagstat_qc['singletons'][1]),

            'singletons_pct':   '{:.2%}'.format(
                float(flagstat_qc['singletons'][0]) /
                float(flagstat_qc['in_total'][0])),

            'diff_chroms':
                int(flagstat_qc['mate_mapped_different_chr_hiQ'][0]),

            'diff_chroms_qc_failed':
                int(flagstat_qc['mate_mapped_different_chr_hiQ'][1])
        })
    obj.update(COMMON_METADATA)

    return obj


def samtools_flagstats_quality_metric(step_run, stages, files):
    logger.debug(
        "in samtools_flagstats_quality_metric with "
        "step_run %s stages.keys() %s output files %s"
        % (step_run, stages.keys(), files))

    file_accessions = list(set(flat([
        resolve_name_to_accessions(stages, output_name)
        for output_name in files])))

    quality_metric_objects = []

    if any([stage_name.startswith('Map ENCSR')
            for stage_name in stages.keys()]):
        quality_metric_objects.append(get_flagstat_obj(
            step_run,
            next(stages[stage_name]['stage_metadata']
                 for stage_name in stages.keys()
                 if stage_name.startswith("Map ENCSR")),
            file_accessions))

    if any([stage_name.startswith('Filter and QC')
            for stage_name in stages.keys()]) \
       and 'scrubbed_unfiltered_bam' not in files:
        quality_metric_objects.append(get_flagstat_obj(
                step_run,
                next(stages[stage_name]['stage_metadata']
                     for stage_name in stages.keys()
                     if stage_name.startswith("Filter and QC")),
                file_accessions))

    return quality_metric_objects


def idr_quality_metric(step_run, stages, files):

    logger.debug(
        "in idr_filter_quality_metric with "
        "step_run %s stages.keys() %s output files %s"
        % (step_run, stages.keys(), files))

    file_accessions = list(set(flat([
        resolve_name_to_accessions(stages, output_name)
        for output_name in files])))

    final_idr_stage_output = \
        stages['Final IDR peak calls']['stage_metadata']['output']

    def IDR_plot(stage_name):
        return get_attachment(
            stages[stage_name]['stage_metadata']['output']['IDR2_plot'])

    def IDR_params(stage_name):
        return get_attachment(
            stages[stage_name]['stage_metadata']['output']['EM_parameters_log'])

    def IDR_threshold(stage_name):
        return float(
            stages[stage_name]['stage_metadata']['originalInput']['idr_threshold'])

    obj = {
        # 'assay_term_id':     'OBI:0000716',
        'assay_term_name':   'ChIP-seq',
        'step_run':          step_run,
        'quality_metric_of': file_accessions,
    }

    # this is just a cheap way to detect replicated vs unreplicated experiment
    if final_idr_stage_output.get('rescue_ratio'):
        obj.update({
            'N1': int(final_idr_stage_output['N1']),
            'N2': int(final_idr_stage_output['N2']),
            'Np': int(final_idr_stage_output['Np']),
            'Nt': int(final_idr_stage_output['Nt']),

            'self_consistency_ratio':
                float(final_idr_stage_output['self_consistency_ratio']),
            'rescue_ratio':
                float(final_idr_stage_output['rescue_ratio']),
            'reproducibility_test':
                str(final_idr_stage_output['reproducibility_test']),

            'IDR_plot_true':    IDR_plot('IDR True Replicates'),
            'IDR_plot_rep1_pr': IDR_plot('IDR Rep 1 Self-pseudoreplicates'),
            'IDR_plot_rep2_pr': IDR_plot('IDR Rep 2 Self-pseudoreplicates'),
            'IDR_plot_pool_pr': IDR_plot('IDR Pooled Pseudoreplicates'),

            'IDR_parameters_true':    IDR_params('IDR True Replicates'),
            'IDR_parameters_rep1_pr': IDR_params('IDR Rep 1 Self-pseudoreplicates'),
            'IDR_parameters_rep2_pr': IDR_params('IDR Rep 2 Self-pseudoreplicates'),
            'IDR_parameters_pool_pr': IDR_params('IDR Pooled Pseudoreplicates')
        })
        # Only accession FRiP scores if calculated
        if final_idr_stage_output.get('F1'):
            obj.update({
                'F1': float(final_idr_stage_output['F1']),
                'F2': float(final_idr_stage_output['F2']),
                'Fp': float(final_idr_stage_output['Fp']),
                'Ft': float(final_idr_stage_output['Ft'])
            })
        # these were not surfaced as outputs in earlier versions of the
        # ENCODE IDR applet, so need to check first if they're there
        if 'No' in final_idr_stage_output:
            obj.update({'N_optimal': final_idr_stage_output['No']})
        if 'Nc' in final_idr_stage_output:
            obj.update({'N_conservative': final_idr_stage_output['Nc']})

        # IDR cutoff should be the same for all IDR stages
        idr_cutoffs = \
            [IDR_threshold(stage_name) for stage_name in
                ['IDR True Replicates', 'IDR Rep 1 Self-pseudoreplicates',
                 'IDR Rep 2 Self-pseudoreplicates', 'IDR Pooled Pseudoreplicates']]
        if all(x == idr_cutoffs[0] for x in idr_cutoffs):
            obj.update({'IDR_cutoff': idr_cutoffs[0]})
        else:  # this is a serious enough error to block creation of the object
            logger.error('Unequal IDR cutoffs: %s' % (idr_cutoffs))
            return None

    else:
        obj.update({
            'N1': int(final_idr_stage_output['N1']),
            'IDR_plot_rep1_pr': IDR_plot('IDR Rep 1 Self-pseudoreplicates'),
            'IDR_parameters_rep1_pr': IDR_params('IDR Rep 1 Self-pseudoreplicates'),
            'IDR_cutoff': IDR_threshold('IDR Rep 1 Self-pseudoreplicates')
        })
        # Only accession FRiP scores if calculated
        if final_idr_stage_output.get('F1'):
            obj.update({
                'F1': float(final_idr_stage_output['F1'])
            })

    obj.update(COMMON_METADATA)

    return [obj]


def histone_chipseq_quality_metric(step_run, stages, files):
    logger.debug(
        "in histone_chipseq_quality_metric with "
        "step_run %s stages.keys() %s output files %s"
        % (step_run, stages.keys(), files))

    file_accessions = list(set(flat([
        resolve_name_to_accessions(stages, output_name)
        for output_name in files])))

    histone_overlap_stage_output = \
        stages['Final narrowpeaks']['stage_metadata']['output']

    obj = {
        # 'assay_term_id':     'OBI:0000716',
        'assay_term_name':   'ChIP-seq',
        'step_run':          step_run,
        'quality_metric_of': file_accessions,
    }

    if not histone_overlap_stage_output.get('rep1_signal'):
        raise AccessioningError(
            "Expected stage 'Final narrowpeaks' to have rep1_signal output.  Found these outputs:\n%s"
            % (histone_overlap_stage_output.keys()))
    replicated_analysis = all([
        histone_overlap_stage_output.get(s) for s in [
            'rep1_signal',
            'rep2_signal']
        ])
    logger.debug("replicated_analysis %s" % (replicated_analysis))
    qc_data_mapping_DX_ENCODE = {
        'frip_score': ('Ft' if replicated_analysis else 'F1', float),
        'frip_nreads': ('nreads', int),
        'frip_nreads_in_peaks': ('nreads_in_peaks', int),
        'npeaks_out': ('npeak_overlap', int)
    }

    for dx_output_name, encode_property in qc_data_mapping_DX_ENCODE.iteritems():
        # Only accession FRiP scores if calculated
        encode_property_name = encode_property[0]
        encode_property_type = encode_property[1]
        if histone_overlap_stage_output.get(dx_output_name):
            obj.update({
                encode_property_name:
                    encode_property_type(histone_overlap_stage_output[dx_output_name])
            })

    obj.update(COMMON_METADATA)

    return [obj]


def dxf_md5(dx_fh):
    logger.debug(
        "in dxf_md5 with handler %s with name %s"
        % (dx_fh, dx_fh.name))
    if 'md5sum' in dx_fh.get_properties():
        md5sum = dx_fh.get_properties()['md5sum']
    else:
        with tempfile.NamedTemporaryFile(delete=True) as tmpfile:
            dxpy.download_dxfile(dx_fh.get_id(), tmpfile.name)
            md5sum = common.md5(tmpfile.name)
        try:
            set_property(dx_fh, {'md5sum': md5sum})
        except Exception as e:
            logger.warning(
                '%s: skipping adding md5sum property to %s.' % (e, dx_fh.name))
    logger.debug('exiting dxf_md5 with %s' % (md5sum))
    return md5sum


def dxf_content_md5(dx_fh):
    logger.debug(
        "in dxf_content_md5 with handler %s with name %s"
        % (dx_fh, dx_fh.name))
    with tempfile.NamedTemporaryFile(delete=True) as tmpfile:
        dxpy.download_dxfile(dx_fh.get_id(), tmpfile.name)
        from magic import from_file
        compressed_mimetypes = [
            "application/x-compress",
            "application/x-bzip2",
            "application/x-gzip"
            ]
        mime_type = from_file(tmpfile.name, mime=True)
        if mime_type in compressed_mimetypes:
            with tempfile.NamedTemporaryFile(delete=True) as uncompressed_tmpfile:
                out, err = common.run_pipe([
                    'cat %s' % (tmpfile.name),
                    'gzip -d'], uncompressed_tmpfile.name)
                md5sum = common.md5(uncompressed_tmpfile.name)
        else:
            md5sum = common.md5(tmpfile.name)
    logger.debug('exiting dxf_content_md5 with %s' % (md5sum))
    return md5sum


def get_rep_fastqs(experiment, keypair, server, repn):
    fastq_valid_status = ['released', 'in progress', 'uploaded']
    logger.debug('in get_rep_fastqs with experiment[accession] %s rep %d'
                 % (experiment.get('accession'), repn))

    original_files = \
        [common.encoded_get(urlparse.urljoin(server, '%s' % (uri)), keypair)
         for uri in experiment.get('original_files')]

    fastqs = \
        [f for f in original_files
         if f.get('file_format') in ['fastq', 'fasta'] and
         f.get('status') in fastq_valid_status]

    # resolve the biorep_n for each fastq
    rep_fastqs = \
        [f for f in fastqs
         if common.encoded_get(
            urlparse.urljoin(
                server, '%s'
                % (f.get('replicate'))), keypair).get(
                'biological_replicate_number') == repn]
    logger.debug('get_rep_fastqs returning %s'
                 % ([f.get('accession') for f in rep_fastqs]))
    return rep_fastqs


def get_stage_name(pattern, stages):
    if not isinstance(stages, list):
        stages = [stages]
    logger.debug('in get_stage_name with stages %s and pattern %s'
                 % ([stage.get('name') for stage in stages], pattern))
    return next(
        re.match(pattern, stage['name']).group(0)
        for stage in stages
        if re.match(pattern, stage['name'])) or None


def get_stage_metadata(analysis, stage_pattern):
    logger.debug('in get_stage_metadata with analysis %s and stage_pattern %s'
                 % (analysis.get('name'), stage_pattern))
    # unfortunately, very early runs of the IDR pipeline mispelled
    # a stage.  We still need to go back to those runs to harvest QC or for
    # other reasons, so here we just special-case it then over-ride the
    # error.
    try:
        stage_metadata = \
            next(s['execution'] for s in analysis.get('stages')
                 if re.match(stage_pattern, s['execution']['name']))
    except StopIteration:
        if stage_pattern == "IDR Pooled Pseudoreplicates":
            tmp_metadata = \
                get_stage_metadata(analysis, "IDR Pooled Pseudoeplicates")
            tmp_metadata['name'] = "IDR Pooled Pseudoreplicates"
            return tmp_metadata
        else:
            raise
    except:
        raise
    else:
        return stage_metadata


def get_experiment_accession(analysis_or_id):

    if isinstance(analysis_or_id, dict):
        analysis = analysis_or_id
    else:  # assumed to be an analysis id
        analysis = dxpy.describe(analysis_or_id)

    m_executableName = \
        re.search('(ENCSR[0-9]{3}[A-Z]{3})', analysis['executableName'])

    m_name = \
        re.search('(ENCSR[0-9]{3}[A-Z]{3})', analysis['name'])

    if not (m_executableName or m_name):
        logger.error("No experiment accession in name %s or executableName %s."
                     % (analysis['name'], analysis['executableName']))
        return
    elif (m_executableName and m_name):
        executableName_accession = m_executableName.group(1)
        name_accession = m_name.group(1)
        if executableName_accession == name_accession:
            return executableName_accession
        else:
            logger.error(
                'Different experiment accessions: name %s, executableName %s.'
                % (analysis['name'], analysis['executableName']))
            return None
    else:
        m = (m_executableName or m_name)
        experiment_accession = m.group(1)
        logger.debug("get_experiment_accession returning %s"
                     % (experiment_accession))
        return experiment_accession


def get_encoded_repn(mapping_analysis):
    # this is a fragile way to infer the rep number.  It depends on the name of
    # the mapping analysis.  But since there is nowhere in the anlysis input
    # to put it, the only alternative is to infer it from the list of fastq
    # accessions in the inputs.
    # Probably the best thing is to add that as an input into the analysis,
    # then it would be recorded in the analysis metadata
    logger.debug("in get_encoded_repn with analysis[name] %s"
                 % (mapping_analysis['name']))

    m_name = re.search(
        'Map ENCSR[0-9]{3}[A-Z]{3} rep(\d+)', mapping_analysis['name'])

    if not m_name:
        logger.error("Could not infer ENCODED repn from analysis name: %s"
                     % (mapping_analysis['name']))
        return
    else:
        logger.debug("in get_encoded_repn and found repn to be %s"
                     % (m_name.group(1)))
        encoded_repn = int(m_name.group(1))
        return encoded_repn


def is_unreplicated_analysis(analysis):
    return (
        analysis['properties'].get('unreplicated_experiment') in ['True', 'true']
        or analysis['properties'].get('simplicate_experiment') in ['True', 'true'])


def is_unary_control(analysis):
    return analysis['properties'].get('unary_control') in ['True', 'true']


def scrubbed_stage(stage):
    logger.debug('in scrubbed_stage with stage %s', pprint.pformat(stage, depth=3))
    return stage['input'].get('scrub')


def get_raw_mapping_stages(mapping_analysis, keypair, server, fqcheck, repn):
    logger.debug(
        'in get_raw_mapping_stages with mapping analysis %s and rep %s'
        % (mapping_analysis['id'], repn))

    experiment_accession = get_experiment_accession(mapping_analysis)

    experiment = common.encoded_get(
        urlparse.urljoin(server, '/experiments/%s'
                                 % (experiment_accession)), keypair)

    # This encoded_repn is the biological_replicate_number at ENCODEd,
    # which needs to be puzzled out from the mapping analysis name
    # or, better, by inferring the rep number from the fastqs actually
    # imported into the analysis
    encoded_repn = get_encoded_repn(mapping_analysis)

    experiment_fastqs = \
        get_rep_fastqs(experiment, keypair, server, encoded_repn)

    experiment_fastq_accessions = \
        [f.get('accession') for f in experiment_fastqs]

    logger.info('%s: Found accessioned experiment fastqs with accessions %s'
                % (experiment_accession, experiment_fastq_accessions))

    analysis_stages = \
        [stage['execution'] for stage in mapping_analysis.get('stages')]

    input_stage = \
        next(stage for stage in analysis_stages
             if stage['name'].startswith("Gather inputs"))
    logger.debug("input_stage['input'] JSON:")
    logger.debug(pprint.pformat(input_stage['input']))

    input_fastq_accessions = []
    input_fastq_accessions.extend(input_stage['input']['reads1'])
    logger.debug(
        'reads1 only input_fastq_accessions %s'
        % (input_fastq_accessions))

    if input_stage['input']['reads2']:
        # Coerce into a list here because in earlier versions of the pipeline
        # code, reads2 was just a string.
        reads2 = input_stage['input']['reads2']
        logger.debug('found reads2 %s' % (reads2))
        if type(reads2) is list:
            input_fastq_accessions.extend(reads2)
        else:
            input_fastq_accessions.extend([reads2])

    logger.debug(
        'reads1 and reads2 input_fastq_accessions %s'
        % (input_fastq_accessions))

    fastqs = []
    for acc in input_fastq_accessions:
        fobj = common.encoded_get(
            urlparse.urljoin(server, 'files/%s' % (acc)), keypair)
        fastqs.append(fobj)

    logger.info('Found input fastq objects with accessions %s'
                % ([f.get('accession') for f in fastqs]))

    # Error if it appears we're trying to accession an out-dated analysis
    # (i.e. one not derived from proper fastqs ... maybe some added or revoked)
    if fqcheck:
        if cmp(sorted(flat(experiment_fastq_accessions)),
               sorted(flat(input_fastq_accessions))):
            fastqs_match = False
            assert fastqs_match, (
                '%s rep%s: Accessioned experiment fastqs differ from analysis.'
                % (experiment_accession, repn) +
                'Suppress with fqcheck=False')
            return None
    else:
        logger.warning(
            '--fqcheck is False, '
            'so not checking to see if experiment and mapped fastqs match')

    raw_mapping_stage = next(
        stage for stage in analysis_stages
        if stage['name'].startswith("Map ENCSR"))
    filter_qc_stage = next(
        (stage for stage in analysis_stages
         if stage['name'].startswith("Filter and QC")), None)
    if filter_qc_stage:
        scrubbed = any([scrubbed_stage(stage) for stage in analysis_stages])
    else:
        scrubbed = False
    logger.debug('scrubbed is %s' % (scrubbed))
    crop_length = raw_mapping_stage['output'].get('crop_length')

    if not crop_length or crop_length == 'native':
        logger.warning(
            'crop_length %s. Inferring mapped_read_length from fastqs'
            % (crop_length))
        native_lengths = set([fq.get('read_length') for fq in fastqs])
        try:
            assert (len(native_lengths) == 1 and
                    all([isinstance(rl, int) for rl in native_lengths])), \
                   ('fastqs with different or non-integer read_lengths: %s'
                    % ([(fq.get('accession'), fq.get('read_length'))
                        for fq in fastqs]))
        except AssertionError:
            if fqcheck:
                raise
            else:
                logger.warning(
                    'fastqs with different or non-integer read_lengths: %s But fqcheck is False so ignoring'
                    % ([(fq.get('accession'), fq.get('read_length'))
                        for fq in fastqs]))
        except:
            raise
        mapped_read_length = int(next(l for l in native_lengths))
    else:
        mapped_read_length = int(crop_length)

    # here we get the actual DNAnexus file that was used as the reference
    # need to remain backwards-compatible with analyses that used output_JSON
    input_stage_output = \
        input_stage['output'].get('output_JSON') or input_stage['output']
    reference_file = dxpy.describe(input_stage_output['reference_tar'])

    # and construct the alias to find the corresponding file at ENCODEd
    reference_alias = "dnanexus:" + reference_file.get('id')

    logger.debug('looking for reference file with alias %s'
                 % (reference_alias))

    reference = common.encoded_get(
        urlparse.urljoin(server, 'files/%s' % (reference_alias)), keypair)
    assert reference, "Reference file %s not found on Portal" % (reference_alias) 
    logger.debug('found reference file %s' % (reference.get('accession')))

    bam_metadata = common.merge_dicts({
        'file_format': 'bam',
        'output_type': 'unfiltered alignments' if not scrubbed else 'redacted unfiltered alignments',
        'assembly': reference.get('assembly'),
        'mapped_read_length': mapped_read_length
    }, COMMON_METADATA)

    if scrubbed:
        mapping_stages = {
            get_stage_name("Map ENCSR.*", analysis_stages): {
                'input_files': [],
                'output_files': [],
                'qc': [],
                'stage_metadata': {}  # initialized below
            },
            get_stage_name("Filter and QC.*", analysis_stages): {
                'input_files': [

                    {'name': 'rep%s_fastqs' % (repn),
                     'derived_from': None,
                     'metadata': None,
                     'encode_object': fastqs},

                    {'name': 'reference',
                     'derived_from': None,
                     'metadata': None,
                     'encode_object': reference}

                ],

                'output_files': [
                    {'name': 'scrubbed_unfiltered_bam',
                     'derived_from': ['rep%s_fastqs' % (repn), 'reference'],
                     'metadata': bam_metadata}
                ],
                'qc': [qc],

                'stage_metadata': {}  # initialized below
            }
        }
    else:
        mapping_stages = {
            get_stage_name("Map ENCSR.*", analysis_stages): {
                'input_files': [

                    {'name': 'rep%s_fastqs' % (repn),
                     'derived_from': None,
                     'metadata': None,
                     'encode_object': fastqs},

                    {'name': 'reference',
                     'derived_from': None,
                     'metadata': None,
                     'encode_object': reference}

                ],

                'output_files': [

                    {'name': 'mapped_reads',
                     'derived_from': ['rep%s_fastqs' % (repn), 'reference'],
                     'metadata': bam_metadata}

                ],

                'qc': [qc],

                'stage_metadata': {}  # initialized below
            }
        }

    for stage_name in mapping_stages:
        if not stage_name.startswith('_'):
            mapping_stages[stage_name].update(
                {'stage_metadata': get_stage_metadata(
                    mapping_analysis, stage_name)})

    return mapping_stages


def get_mapping_stages(mapping_analysis, keypair, server, fqcheck, repn):
    logger.debug(
        'in get_mapping_stages with mapping analysis %s and rep %s'
        % (mapping_analysis['id'], repn))

    if not mapping_analysis:
        logger.warning(
            'get_mapping_stages got empty mapping_analysis, returning None')
        return None
    experiment_accession = get_experiment_accession(mapping_analysis)
    url = urlparse.urljoin(
            server, '/experiments/%s' % (experiment_accession))
    r = common.encoded_get(url, keypair, return_response=True)
    r.raise_for_status()
    experiment = r.json()

    # This encoded_repn is the biological_replicate_number at ENCODEd, which
    # needs to be puzzled out from the mapping analysis name or, better, by
    # inferring the rep number from the fastqs actually imported into the
    # analysis
    encoded_repn = get_encoded_repn(mapping_analysis)

    experiment_fastqs = \
        get_rep_fastqs(experiment, keypair, server, encoded_repn)

    experiment_fastq_accessions = \
        [f.get('accession') for f in experiment_fastqs]

    logger.info('%s: Found accessioned experiment fastqs with accessions %s'
                % (experiment_accession, experiment_fastq_accessions))

    analysis_stages = \
        [stage['execution'] for stage in mapping_analysis.get('stages')]

    input_stage = \
        next(stage for stage in analysis_stages
             if stage['name'].startswith("Gather inputs"))
    logger.debug("input_stage['input'] JSON:")
    logger.debug(pprint.pformat(input_stage['input']))

    input_fastq_accessions = []
    input_fastq_accessions.extend(input_stage['input']['reads1'])
    logger.debug(
        'reads1 only input_fastq_accessions %s'
        % (input_fastq_accessions))

    if input_stage['input']['reads2']:
        # Coerce into a list here because in earlier versions of the pipeline
        # code, reads2 was just a string.
        reads2 = input_stage['input']['reads2']
        logger.debug('found reads2 %s' % (reads2))
        if type(reads2) is list:
            input_fastq_accessions.extend(reads2)
        else:
            input_fastq_accessions.extend([reads2])

    logger.debug(
        'reads1 and reads2 input_fastq_accessions %s'
        % (input_fastq_accessions))

    fastqs = []
    for acc in input_fastq_accessions:
        fobj = common.encoded_get(
            urlparse.urljoin(server, 'files/%s' % (acc)), keypair)
        fastqs.append(fobj)

    logger.info('Found input fastq objects with accessions %s'
                % ([f.get('accession') for f in fastqs]))

    # Error if it appears we're trying to accession an out-dated analysis
    # (i.e. one not derived from proper fastqs ... maybe some added or revoked)
    if fqcheck:
        if cmp(sorted(flat(experiment_fastq_accessions)),
               sorted(flat(input_fastq_accessions))):
            fastqs_match = False
            assert fastqs_match, (
                '%s rep%s: Accessioned experiment fastqs differ from analysis.'
                % (experiment_accession, repn) +
                'Suppress with fqcheck=False')
            return None
    else:
        logger.warning(
            '--fqcheck is False, '
            'so not checking to see if experiment and mapped fastqs match')

    raw_mapping_stage = next(
        stage for stage in analysis_stages
        if stage['name'].startswith("Map ENCSR"))
    filter_qc_stage = next(
        stage for stage in analysis_stages
        if stage['name'].startswith("Filter and QC"))
    scrubbed = any([scrubbed_stage(stage) for stage in analysis_stages])
    crop_length = raw_mapping_stage['output'].get('crop_length')

    if not crop_length or crop_length == 'native':
        logger.warning(
            'crop_length %s. Inferring mapped_read_length from fastqs'
            % (crop_length))
        native_lengths = set([fq.get('read_length') for fq in fastqs])
        try:
            assert (len(native_lengths) == 1 and
                    all([isinstance(rl, int) for rl in native_lengths])), \
                   ('fastqs with different or non-integer read_lengths: %s'
                    % ([(fq.get('accession'), fq.get('read_length'))
                        for fq in fastqs]))
        except AssertionError:
            if fqcheck:
                raise
            else:
                logger.warning(
                    'fastqs with different or non-integer read_lengths: %s But fqcheck is False so ignoring'
                    % ([(fq.get('accession'), fq.get('read_length'))
                        for fq in fastqs]))
        except:
            raise
        mapped_read_length = int(next(l for l in native_lengths))
    else:
        mapped_read_length = int(crop_length)

    # here we get the actual DNAnexus file that was used as the reference
    # need to remain backwards-compatible with analyses that used output_JSON
    input_stage_output = \
        input_stage['output'].get('output_JSON') or input_stage['output']
    reference_file = dxpy.describe(input_stage_output['reference_tar'])

    # and construct the alias to find the corresponding file at ENCODEd
    reference_alias = "dnanexus:" + reference_file.get('id')

    logger.debug('looking for reference file with alias %s'
                 % (reference_alias))

    reference = common.encoded_get(
        urlparse.urljoin(server, 'files/%s' % (reference_alias)), keypair)
    assert reference, "Reference file %s not found on Portal" % (reference_alias) 
    logger.debug('found reference file %s' % (reference.get('accession')))

    bam_metadata = common.merge_dicts({
        'file_format': 'bam',
        'output_type': 'alignments' if not scrubbed else 'redacted alignments',
        'mapped_read_length': mapped_read_length,
        'assembly': reference.get('assembly')
    }, COMMON_METADATA)

    mapping_stages = {

        get_stage_name("Map ENCSR.*", analysis_stages): {
            'input_files': [],
            'output_files': [],
            'qc': [],
            'stage_metadata': {}  # initialized below
        },

        get_stage_name("Filter and QC.*", analysis_stages): {
            'input_files': [

                {'name': 'rep%s_fastqs' % (repn),
                 'derived_from': None,
                 'metadata': None,
                 'encode_object': fastqs},

                {'name': 'reference',
                 'derived_from': None,
                 'metadata': None,
                 'encode_object': reference}

            ],

            'output_files': [
                {'name': 'scrubbed_filtered_bam' if scrubbed else 'filtered_bam',
                 'derived_from': ['rep%s_fastqs' % (repn), 'reference'],
                 'metadata': bam_metadata}
            ],
            'qc': [qc, dup_qc, pbc_qc, filtered_qc, xcor_qc],

            'stage_metadata': {}  # initialized below
        },

        get_stage_name("Calculate cross-correlation.*", analysis_stages): {
            'input_files': [],
            'output_files': [],
            'qc': [],
            'stage_metadata': {}
        }
    }

    for stage_name in mapping_stages:
        if not stage_name.startswith('_'):
            mapping_stages[stage_name].update(
                {'stage_metadata': get_stage_metadata(
                    mapping_analysis, stage_name)})

    return mapping_stages


def get_control_mapping_stages(peaks_analysis, keypair, server, fqcheck):
    # Find the control inputs
    logger.debug(
        'in get_control_mapping_stages with peaks_analysis %s'
        % (peaks_analysis['id']))

    peaks_stages = \
        [stage['execution'] for stage in peaks_analysis.get('stages')]

    peaks_stage = next(
        stage for stage in peaks_stages
        if stage['name'] == "ENCODE Peaks")

    # here reps is always 1,2 because the peaks job rep numbers are 1,2 ...
    # these are not the ENCODEd biological_replicate_numbers, which are
    # only known to the analysis via its name, or by going back to ENCODEd and
    # figuring out where the fastqs came from
    reps = sorted([
        re.match("ctl(\d+)_ta", input_key).group(1)
        for input_key in peaks_stage['input'].keys()
        if re.match("ctl(\d+)_ta", input_key)])

    mapping_analyses = []
    for repn in reps:

        ta = dxpy.describe(peaks_stage['input']['ctl%s_ta' % (repn)])
        if not ta:
            raise AccessioningError("Could not find contol tagAlign for control rep %s" % (repn))

        mapping_job_id = ta['createdBy']['job']
        if not mapping_job_id:
            raise AccessioningError("Failed to find job id that created tagAlign for control rep %s.  If this analysis used a custom control consider running with --skip_control" % (repn))

        mapping_job = dxpy.describe(mapping_job_id)
        mapping_analysis_id = mapping_job.get('analysis')
        if not mapping_analysis_id:
            raise AccessioningError("Failed to find mapping analysis ids that created tagAlign for control rep %s.  If this analysis used a custom control consider running with --skip_control" % (repn))

        mapping_analyses.append(dxpy.describe(mapping_analysis_id))

    mapping_stages = []

    for (i, repn) in enumerate(reps):
        mapping_stage = get_mapping_stages(
            mapping_analyses[i], keypair, server, fqcheck, repn)

        if not mapping_stage:
            logger.warning(
                '%s: failed to find mapping stages for rep%d'
                % (peaks_analysis['id'], repn))
        mapping_stages.append(mapping_stage)

    return mapping_stages


def get_peak_mapping_stages(peaks_analysis, keypair, server,
                            fqcheck):

    # Find the tagaligns actually used as inputs into the analysis
    # Find the mapping analyses that produced those tagaligns
    # Find the filtered bams from those analyses
    # Build the stage dict and return it

    logger.debug(
        'in get_peak_mapping_stages: peaks_analysis is %s named %s'
        % (peaks_analysis.get('id'), peaks_analysis.get('name')))

    if is_unreplicated_analysis(peaks_analysis):
        reps = [1]
    else:
        reps = [1, 2]

    peaks_stages = \
        [stage['execution'] for stage in peaks_analysis.get('stages')]

    peaks_stage = next(
        stage for stage in peaks_stages
        if stage['name'] == "ENCODE Peaks")

    tas = [dxpy.describe(peaks_stage['input']['rep%s_ta' % (n)])
           for n in reps]

    mapping_jobs = \
        [dxpy.describe(ta['createdBy']['job'])
         for ta in tas]

    mapping_analyses = \
        [dxpy.describe(mapping_job['analysis'])
         for mapping_job in mapping_jobs if mapping_job]

    mapping_stages = []
    for (i, repn) in enumerate(reps):
        mapping_stage = \
            get_mapping_stages(
                mapping_analyses[i], keypair, server, fqcheck, repn)
        if not mapping_stage:
            logger.error('%s: failed to find mapping stages for rep%d'
                         % (peaks_analysis['id'], repn))
            return None
        else:
            mapping_stages.append(mapping_stage)

    return mapping_stages


def pooled_controls(peaks_analysis, rep):
    # this is not surfaced explicitly so must be inferred
    # General:  get the id's of the files actually used for the specified rep
    # and pooled controls.  If the id is the same as the pooled control id then
    # return true.
    # Specifically:
    # starting with the peaks_analysis, get its stages
    # get "ENCODE Peaks" stage
    # get the job id for "ENCODE Peaks"
    # get the control and experiment file ID's for the specified rep
    # find the child jobs of the "ENCODE Peaks" job
    # find the child job where macs2 was run with the experiment file
    # corresponding to the experiment file for this rep
    # get from that child job the file ID of the control
    # if the contol file ID for this rep from ENCODE Peaks is the same as in
    # macs2 then return False else return True
    # Could double-check the log output of ENCODE Peaks to search for the
    # strings "Using pooled controls for replicate 1.", "Using pooled controls
    # for replicate 2." and "Using pooled controls."
    # But there is no corresponding "Not pooling controls"
    # message, so it's underdertermined.

    logger.debug('in pooled_controls with peaks_analysis %s; rep %s'
                 % (peaks_analysis['id'], rep))

    peaks_stages = peaks_analysis.get('stages')

    ENCODE_Peaks_stage = next(
        stage['execution'] for stage in peaks_stages
        if stage['execution']['name'] == "ENCODE Peaks")

    ENCODE_Peaks_exp_file = \
        ENCODE_Peaks_stage['input']['rep%s_ta' % (rep)]

    ENCODE_Peaks_ctl_file = \
        ENCODE_Peaks_stage['input']['ctl%s_ta' % (rep)]

    child_jobs = dxpy.find_jobs(
        parent_job=ENCODE_Peaks_stage['id'],
        name="MACS2",
        project=ENCODE_Peaks_stage['project'],
        describe=True)

    rep_job = next(
        job for job in child_jobs
        if job['describe']['input']['experiment'] == ENCODE_Peaks_exp_file)

    # for job in child_jobs:
    #   #pprint.pprint(job)
    #   if job['describe']['input']['experiment'] == ENCODE_Peaks_exp_file:
    #       rep_job = job

    rep_job_ctl_file = rep_job['describe']['input']['control']

    logger.info("Rep%s input control file %s; actually used %s"
                % (rep, ENCODE_Peaks_ctl_file, rep_job_ctl_file))

    if ENCODE_Peaks_ctl_file == rep_job_ctl_file:
        logger.info('Inferred controls not pooled for rep%s' % (rep))
        return False
    else:
        logger.info('Inferred pooled controls for rep%s' % (rep))
        return True


def get_assembly(stage_output_tuple):
    stages, output_key = stage_output_tuple
    logger.debug('stages %s' % (stages))
    logger.debug(
        'in get_assembly with output_key %s and stages:\n%s'
        % (output_key, pprint.pformat([stage for stage in stages or []])))
    if not stages:
        return None
    for stage in stages.itervalues():
        output_files = stage.get('output_files')
        if output_files:
            logger.debug('output_files: %s' % (pprint.pformat([output_file for output_file in output_files])))
            output_file_metadata = next(
                output_file.get('metadata')
                for output_file in output_files
                if output_file.get('name') == output_key)
            return output_file_metadata.get('assembly')


def get_histone_peak_stages(peaks_analysis, mapping_stages, control_stages,
                            experiment, keypair, server):

    logger.debug(
        'in get_histone_peak_stages with peaks_analysis %s;'
        % (peaks_analysis['id']) +
        'experiment %s and len(mapping_stages) %d len(control_stages) %d'
        % (experiment['accession'], len(mapping_stages), len(control_stages)))
    unreplicated_analysis = is_unreplicated_analysis(peaks_analysis)

    bams = \
        [(mapping_stage,
          'scrubbed_filtered_bam' if any([scrubbed_stage(stage) for stage in [mapping_stage.get(stage_name).get('stage_metadata') for stage_name in mapping_stage.keys()]])
          else
          'filtered_bam') for mapping_stage in mapping_stages]

    ctl_bams = \
        [(control_stage,
          'scrubbed_filtered_bam' if any([scrubbed_stage(stage) for stage in [control_stage.get(stage_name).get('stage_metadata') for stage_name in control_stage.keys()]])
          else
          'filtered_bam') for control_stage in control_stages]

    assemblies = \
        [get_assembly(bam)
         for bam in bams + ctl_bams
         if get_assembly(bam)]
    observed_assemblies = set(assemblies)
    assert len(observed_assemblies) == 1, "Different bam assemblies for rep1,2 and control rep1,2 bams: %s" % (assemblies)
    assembly = observed_assemblies.pop()

    if not ctl_bams:
        rep1_ctl = []
        rep2_ctl = []
    else:
        if pooled_controls(peaks_analysis, rep=1):
            rep1_ctl = ctl_bams
        else:
            rep1_ctl = [ctl_bams[0]]

        if not unreplicated_analysis:
            if pooled_controls(peaks_analysis, rep=2):
                rep2_ctl = ctl_bams
            else:
                rep2_ctl = [ctl_bams[1]]

    analysis_stages = \
        [stage['execution'] for stage in peaks_analysis.get('stages')]

    common_file_metadata = copy.copy(COMMON_METADATA)
    common_file_metadata.update({'assembly': assembly})

    narrowpeak_metadata = common.merge_dicts(
        {'file_format': 'bed',
         'file_format_type': 'narrowPeak',
         'file_format_specifications': ['encode:narrowPeak.as'],
         'output_type': 'peaks'},
        common_file_metadata)

    narrowpeak_bb_metadata = common.merge_dicts({
        'file_format': 'bigBed',
        'file_format_type': 'narrowPeak',
        'file_format_specifications': ['encode:narrowPeak.as'],
        'output_type': 'peaks'},
        common_file_metadata)

    replicated_narrowpeak_metadata = common.merge_dicts(
        {'file_format': 'bed',
         'file_format_type': 'narrowPeak',
         'file_format_specifications': ['encode:narrowPeak.as'],
         'output_type': 'replicated peaks'},
        common_file_metadata)

    replicated_narrowpeak_bb_metadata = common.merge_dicts({
        'file_format': 'bigBed',
        'file_format_type': 'narrowPeak',
        'file_format_specifications': ['encode:narrowPeak.as'],
        'output_type': 'replicated peaks'},
        common_file_metadata)

    stable_narrowpeak_metadata = common.merge_dicts(
        {'file_format': 'bed',
         'file_format_type': 'narrowPeak',
         'file_format_specifications': ['encode:narrowPeak.as'],
         'output_type': 'stable peaks'},
        common_file_metadata)

    stable_narrowpeak_bb_metadata = common.merge_dicts({
        'file_format': 'bigBed',
        'file_format_type': 'narrowPeak',
        'file_format_specifications': ['encode:narrowPeak.as'],
        'output_type': 'stable peaks'},
        common_file_metadata)

    fc_signal_metadata = common.merge_dicts({
        'file_format': 'bigWig',
        'output_type': 'fold change over control'},
        common_file_metadata)

    pvalue_signal_metadata = common.merge_dicts({
        'file_format': 'bigWig',
        'output_type': 'signal p-value'},
        common_file_metadata)

    # This is lame because the repns are hard-wired.
    # Needs to be made more general
    rep1_bam = bams[0]
    if not unreplicated_analysis:
        rep2_bam = bams[1]
    # This is lame because it assumes pooled controls is all the controls
    # Needs to be made to get the controls that were actually pooled
    pooled_ctl_bams = ctl_bams
    peak_stages = {
        # derived_from is by name here, will be patched into the file metadata
        # after all files are accessioned
        # derived_from can also be a tuple of (stages,name) to connect to#
        # files outside of this set of stages
        get_stage_name("ENCODE Peaks", analysis_stages): {
            'output_files': [
                {'name': 'rep1_narrowpeaks',
                 'derived_from': [rep1_bam] + rep1_ctl,
                 'metadata': narrowpeak_metadata},

                {'name': 'rep1_narrowpeaks_bb',
                 'derived_from': ['rep1_narrowpeaks'],
                 'metadata': narrowpeak_bb_metadata},

                {'name': 'rep1_pvalue_signal',
                 'derived_from': [rep1_bam] + rep1_ctl,
                 'metadata': pvalue_signal_metadata},

                {'name': 'rep1_fc_signal',
                 'derived_from': [rep1_bam] + rep1_ctl,
                 'metadata': fc_signal_metadata}

            ] if unreplicated_analysis else [

                {'name': 'rep1_narrowpeaks',
                 'derived_from': [rep1_bam] + rep1_ctl,
                 'metadata': narrowpeak_metadata},

                {'name': 'rep2_narrowpeaks',
                 'derived_from': [rep2_bam] + rep2_ctl,
                 'metadata': narrowpeak_metadata},

                {'name': 'pooled_narrowpeaks',
                 'derived_from': [rep1_bam, rep2_bam] + pooled_ctl_bams,
                 'metadata': narrowpeak_metadata},

                {'name': 'rep1_narrowpeaks_bb',
                 'derived_from': ['rep1_narrowpeaks'],
                 'metadata': narrowpeak_bb_metadata},

                {'name': 'rep2_narrowpeaks_bb',
                 'derived_from': ['rep2_narrowpeaks'],
                 'metadata': narrowpeak_bb_metadata},

                {'name': 'pooled_narrowpeaks_bb',
                 'derived_from': ['pooled_narrowpeaks'],
                 'metadata': narrowpeak_bb_metadata},

                {'name': 'rep1_pvalue_signal',
                 'derived_from': [rep1_bam] + rep1_ctl,
                 'metadata': pvalue_signal_metadata},

                {'name': 'rep2_pvalue_signal',
                 'derived_from': [rep2_bam] + rep2_ctl,
                 'metadata': pvalue_signal_metadata},

                {'name': 'pooled_pvalue_signal',
                 'derived_from': [rep1_bam, rep2_bam] + pooled_ctl_bams,
                 'metadata': pvalue_signal_metadata},

                {'name': 'rep1_fc_signal',
                 'derived_from': [rep1_bam] + rep1_ctl,
                 'metadata': fc_signal_metadata},

                {'name': 'rep2_fc_signal',
                 'derived_from': [rep2_bam] + rep2_ctl,
                 'metadata': fc_signal_metadata},

                {'name': 'pooled_fc_signal',
                 'derived_from': [rep1_bam, rep2_bam] + pooled_ctl_bams,
                 'metadata': fc_signal_metadata}
            ],

            'qc': [],

            'stage_metadata': {}  # initialized below
        },

        # older pipeline versions called this "Overlap" so need a pattern for
        # backward compatibility
        get_stage_name("(Overlap|Final) narrowpeaks", analysis_stages): {
            'output_files': [

                {'name': 'overlapping_peaks',
                 'derived_from': ['rep1_narrowpeaks', 'rep2_narrowpeaks',
                                  'pooled_narrowpeaks'],
                 'metadata': replicated_narrowpeak_metadata},

                {'name': 'overlapping_peaks_bb',
                 'derived_from': ['overlapping_peaks'],
                 'metadata': replicated_narrowpeak_bb_metadata}
            ] if not unreplicated_analysis else [
                {'name': 'overlapping_peaks',
                 'derived_from': ['rep1_narrowpeaks'],
                 'metadata': stable_narrowpeak_metadata},

                {'name': 'overlapping_peaks_bb',
                 'derived_from': ['overlapping_peaks'],
                 'metadata': stable_narrowpeak_bb_metadata}
            ],

            'qc': [
                'npeaks_in', 'npeaks_out', 'npeaks_rejected'
            ],

            'stage_metadata': {}  # initialized below
        }
    }

    for stage_name in peak_stages:
        if not stage_name.startswith('_'):
            peak_stages[stage_name].update(
                {'stage_metadata': get_stage_metadata(
                    peaks_analysis, stage_name)})

    return [peak_stages]


def idr_sets_same(analysis_stages):
    if not isinstance(analysis_stages, list):
        stages = [analysis_stages]
    else:
        stages = analysis_stages

    logger.debug('in idr_sets_same with stages %s'
                 % ([stage.get('name') for stage in stages]))
    pattern = "Final IDR peak calls"
    final_idr_stage = next(
        stage
        for stage in stages
        if re.match(pattern, stage['name'])) or None
    if not final_idr_stage:
        raise AccessioningError("In idr_sets_same: no peak stage named 'Final IDR peak calls'")
    optimal_set = dxpy.DXFile(final_idr_stage['output']['optimal_set'])
    conservative_set = dxpy.DXFile(final_idr_stage['output']['conservative_set'])
    if dxf_content_md5(optimal_set) == dxf_content_md5(conservative_set):
        logger.debug("IDR sets are the same")
        return True
    else:
        logger.debug("IDR sets differ")
        return False


def get_tf_peak_stages(peaks_analysis, mapping_stages, control_stages,
                       experiment, keypair, server, signal_only):

    logger.debug(
        'in get_tf_peak_stages with peaks_analysis %s;'
        % (peaks_analysis['id']) +
        'experiment %s and len(mapping_stages) %d len(control_stages) %d'
        % (experiment['accession'], len(mapping_stages), len(control_stages)))
    unreplicated_analysis = is_unreplicated_analysis(peaks_analysis)

    bams = \
        [(mapping_stage,
          'scrubbed_filtered_bam' if any([scrubbed_stage(stage) for stage in [mapping_stage.get(stage_name).get('stage_metadata') for stage_name in mapping_stage.keys()]])
          else
          'filtered_bam') for mapping_stage in mapping_stages]

    ctl_bams = \
        [(control_stage,
          'scrubbed_filtered_bam' if any([scrubbed_stage(stage) for stage in [control_stage.get(stage_name).get('stage_metadata') for stage_name in control_stage.keys()]])
          else
          'filtered_bam') for control_stage in control_stages]

    assemblies = \
        [get_assembly(bam)
         for bam in bams + ctl_bams
         if get_assembly(bam)]
    observed_assemblies = set(assemblies)
    assert len(observed_assemblies) == 1, "Different (or no) bam assemblies found for rep1,2 and control rep1,2 bams: %s" % (assemblies)
    assembly = observed_assemblies.pop()

    if not ctl_bams:
        rep1_ctl = []
        rep2_ctl = []
    else:
        if pooled_controls(peaks_analysis, rep=1):
            rep1_ctl = ctl_bams
        else:
            rep1_ctl = [ctl_bams[0]]

        if not unreplicated_analysis:
            if pooled_controls(peaks_analysis, rep=2):
                rep2_ctl = ctl_bams
            else:
                rep2_ctl = [ctl_bams[1]]

    analysis_stages = \
        [stage['execution'] for stage in peaks_analysis.get('stages')]

    common_file_metadata = copy.copy(COMMON_METADATA)
    common_file_metadata.update({'assembly': assembly})

    narrowpeak_metadata = common.merge_dicts(
        {'file_format': 'bed',
         'file_format_type': 'narrowPeak',
         'file_format_specifications': ['encode:narrowPeak.as'],
         'output_type': 'peaks and background as input for IDR'},
        common_file_metadata)

    idr_optimal_narrowpeak_metadata = common.merge_dicts(
        {'file_format': 'bed',
         'file_format_type': 'narrowPeak',
         'file_format_specifications': ['encode:narrowPeak.as'],
         'output_type': 'optimal idr thresholded peaks'},
        common_file_metadata)

    idr_conservative_narrowpeak_metadata = common.merge_dicts(
        {'file_format': 'bed',
         'file_format_type': 'narrowPeak',
         'file_format_specifications': ['encode:narrowPeak.as'],
         'output_type': 'conservative idr thresholded peaks'},
        common_file_metadata)

    idr_stable_narrowpeak_metadata = common.merge_dicts(
        {'file_format': 'bed',
         'file_format_type': 'narrowPeak',
         'file_format_specifications': ['encode:narrowPeak.as'],
         'output_type': 'pseudoreplicated idr thresholded peaks'},
        common_file_metadata)

    narrowpeak_bb_metadata = common.merge_dicts(
        {'file_format': 'bigBed',
         'file_format_type': 'narrowPeak',
         'file_format_specifications': ['encode:narrowPeak.as'],
         'output_type': 'peaks and background as input for IDR'},
        common_file_metadata)

    idr_optimal_narrowpeak_bb_metadata = common.merge_dicts(
        {'file_format': 'bigBed',
         'file_format_type': 'narrowPeak',
         'file_format_specifications': ['encode:narrowPeak.as'],
         'output_type': 'optimal idr thresholded peaks'},
        common_file_metadata)

    idr_conservative_narrowpeak_bb_metadata = common.merge_dicts(
        {'file_format': 'bigBed',
         'file_format_type': 'narrowPeak',
         'file_format_specifications': ['encode:narrowPeak.as'],
         'output_type': 'conservative idr thresholded peaks'},
        common_file_metadata)

    idr_stable_narrowpeak_bb_metadata = common.merge_dicts(
        {'file_format': 'bigBed',
         'file_format_type': 'narrowPeak',
         'file_format_specifications': ['encode:narrowPeak.as'],
         'output_type': 'pseudoreplicated idr thresholded peaks'},
        common_file_metadata)

    fc_signal_metadata = common.merge_dicts(
        {'file_format': 'bigWig',
         'output_type': 'fold change over control'},
        common_file_metadata)

    pvalue_signal_metadata = common.merge_dicts(
        {'file_format': 'bigWig',
         'output_type': 'signal p-value'},
        common_file_metadata)

    # This is lame because the repns are hard-wired.
    # Needs to be made more general
    rep1_bam = bams[0]
    if not unreplicated_analysis:
        rep2_bam = bams[1]
    # This is lame because it assumes pooled controls is all the controls
    # Needs to be made to get the controls that were actually pooled
    pooled_ctl_bams = ctl_bams
    peak_stages = {
        get_stage_name("ENCODE Peaks", analysis_stages): {

            'output_files': [

                {'name': 'rep1_pvalue_signal',
                 'derived_from': [rep1_bam] + rep1_ctl,
                 'metadata': pvalue_signal_metadata},

                {'name': 'rep1_fc_signal',
                 'derived_from': [rep1_bam] + rep1_ctl,
                 'metadata': fc_signal_metadata}

            ] if unreplicated_analysis else [

                {'name': 'rep1_pvalue_signal',
                 'derived_from': [rep1_bam] + rep1_ctl,
                 'metadata': pvalue_signal_metadata},

                {'name': 'rep1_fc_signal',
                 'derived_from': [rep1_bam] + rep1_ctl,
                 'metadata': fc_signal_metadata},

                {'name': 'rep2_pvalue_signal',
                 'derived_from': [rep2_bam] + rep2_ctl,
                 'metadata': pvalue_signal_metadata},

                {'name': 'pooled_pvalue_signal',
                 'derived_from': [rep1_bam, rep2_bam] + pooled_ctl_bams,
                 'metadata': pvalue_signal_metadata},

                {'name': 'rep2_fc_signal',
                 'derived_from': [rep2_bam] + rep2_ctl,
                 'metadata': fc_signal_metadata},

                {'name': 'pooled_fc_signal',
                 'derived_from': [rep1_bam, rep2_bam] + pooled_ctl_bams,
                 'metadata': fc_signal_metadata}
            ],

            'qc': [],

            'stage_metadata': {}  # initialized below
        }
    }

    if not signal_only:
        peak_stages.update({
            # derived_from is by name here, will be patched into the file
            # metadata after all files are accessioned
            # derived_from can also be a tuple of (stages,name) to connect to
            # files outside of this set of stages

            get_stage_name("SPP Peaks", analysis_stages): {
                'output_files': [

                    {'name': 'rep1_peaks',
                     'derived_from': [rep1_bam] + rep1_ctl,
                     'metadata': narrowpeak_metadata},

                    {'name': 'rep1_peaks_bb',
                     'derived_from': ['rep1_peaks'],
                     'metadata': narrowpeak_bb_metadata}

                ] if unreplicated_analysis else [

                    {'name': 'rep1_peaks',
                     'derived_from': [rep1_bam] + rep1_ctl,
                     'metadata': narrowpeak_metadata},

                    {'name': 'rep2_peaks',
                     'derived_from': [rep2_bam] + rep2_ctl,
                     'metadata': narrowpeak_metadata},

                    {'name': 'pooled_peaks',
                     'derived_from': [rep1_bam, rep2_bam] + pooled_ctl_bams,
                     'metadata': narrowpeak_metadata},

                    {'name': 'rep1_peaks_bb',
                     'derived_from': ['rep1_peaks'],
                     'metadata': narrowpeak_bb_metadata},

                    {'name': 'rep2_peaks_bb',
                     'derived_from': ['rep2_peaks'],
                     'metadata': narrowpeak_bb_metadata},

                    {'name': 'pooled_peaks_bb',
                     'derived_from': ['pooled_peaks'],
                     'metadata': narrowpeak_bb_metadata}
                ],

                'qc': [],

                'stage_metadata': {}  # initialized below
            }
        })

        peak_stages.update({
            get_stage_name("IDR Rep 1 Self-pseudoreplicates", analysis_stages): {
                'output_files': [],
                'qc': [],
                'stage_metadata': {}  # initialized below
            },

            get_stage_name("IDR True Replicates", analysis_stages): {
                'output_files': [],
                'qc': [],
                'stage_metadata': {}  # initialized below
            },

            get_stage_name("IDR Rep 2 Self-pseudoreplicates", analysis_stages): {
                'output_files': [],
                'qc': [],
                'stage_metadata': {}  # initialized below
            },

            get_stage_name("IDR Pooled Pseudor?eplicates", analysis_stages): {
                'output_files': [],
                'qc': [],
                'stage_metadata': {}  # initialized below
            }
        } if not unreplicated_analysis else {
            get_stage_name("IDR Rep 1 Self-pseudoreplicates", analysis_stages): {
                'output_files': [],
                'qc': [],
                'stage_metadata': {}  # initialized below
            }
        })

        if unreplicated_analysis:
            peak_stages.update({
                get_stage_name("Final IDR peak calls", analysis_stages): {
                    'output_files': [
                        {'name': 'stable_set',
                         'derived_from': ['rep1_peaks'],
                         'metadata': idr_stable_narrowpeak_metadata},

                        {'name': 'stable_set_bb',
                         'derived_from': ['stable_set'],
                         'metadata': idr_stable_narrowpeak_bb_metadata},
                    ],
                    'qc': ['N1', 'Ns'],
                    'stage_metadata': {}  # initialized below
                }
            })
        else:  # replicated analysis
            if idr_sets_same(analysis_stages):
                peak_stages.update({
                    get_stage_name("Final IDR peak calls", analysis_stages): {

                        'output_files': [
                            {'name': 'optimal_set',
                             'derived_from': ['rep1_peaks', 'rep2_peaks', 'pooled_peaks'],
                             'metadata': idr_optimal_narrowpeak_metadata},

                            {'name': 'optimal_set_bb',
                             'derived_from': ['optimal_set'],
                             'metadata': idr_optimal_narrowpeak_bb_metadata}

                        ],
                        'qc': [
                            'reproducibility_test', 'rescue_ratio', 'Np', 'N1',
                            'N2', 'Nt', 'self_consistency_ratio'
                        ],
                        'stage_metadata': {}  # initialized below
                    }
                })
            else:  # idr sets not same and a replicated analysis
                peak_stages.update({
                    get_stage_name("Final IDR peak calls", analysis_stages): {
                        'output_files': [
                            {'name': 'conservative_set',
                             'derived_from': ['rep1_peaks', 'rep2_peaks', 'pooled_peaks'],
                             'metadata': idr_conservative_narrowpeak_metadata},

                            {'name': 'conservative_set_bb',
                             'derived_from': ['conservative_set'],
                             'metadata': idr_conservative_narrowpeak_bb_metadata},

                            {'name': 'optimal_set',
                             'derived_from': ['rep1_peaks', 'rep2_peaks', 'pooled_peaks'],
                             'metadata': idr_optimal_narrowpeak_metadata},

                            {'name': 'optimal_set_bb',
                             'derived_from': ['optimal_set'],
                             'metadata': idr_optimal_narrowpeak_bb_metadata}
                        ],
                        'qc': [
                            'reproducibility_test', 'rescue_ratio', 'Np', 'N1',
                            'N2', 'Nt', 'self_consistency_ratio'
                        ],
                        'stage_metadata': {}  # initialized below
                    }
                })

        final_idr_stage_name = \
            get_stage_name("Final IDR peak calls", analysis_stages)
        final_idr_stage_metadata = \
            get_stage_metadata(peaks_analysis, "Final IDR peak calls")
        final_idr_stage_input = final_idr_stage_metadata.get('input')
        blacklist = final_idr_stage_input.get('blacklist')
        if blacklist:
            blacklist_file = dxpy.describe(blacklist)
            # and construct the alias to find the corresponding file at ENCODEd
            blacklist_alias = "dnanexus:" + blacklist_file.get('id')
            encoded_blacklist = common.encoded_get(
                urlparse.urljoin(
                    server, 'files/%s' % (blacklist_alias)), keypair)
            assert encoded_blacklist, "Blacklist file %s not found on Portal" % (blacklist_alias)
            logger.info(
                'DNAnexus blacklist %s matched to ENCODE file %s'
                % (blacklist_file.get('name'), encoded_blacklist.get('accession')))
            peak_stages[final_idr_stage_name].update({
                'input_files': [
                    {'name': 'blacklist',
                     'derived_from': None,
                     'metadata': None,
                     'encode_object': encoded_blacklist}
                ]
                })
            for output_file in peak_stages[final_idr_stage_name]['output_files']:
                if output_file['name'] in ['conservative_set', 'optimal_set', 'stable_set']:
                    output_file['derived_from'].append('blacklist')
        else:
            logger.info('No blacklist used in this analysis.')

    for stage_name in peak_stages:
        if not stage_name.startswith('_'):
            peak_stages[stage_name].update(
                {'stage_metadata': get_stage_metadata(
                    peaks_analysis, stage_name)})

    return [peak_stages]


def resolve_name_to_accessions(stages, stage_file_name):
    # given a dict of named stages, and the name of one of the stages' outputs,
    # return that output's ENCODE accession number

    logger.debug(
        "in resolve_name_to_accessions with stage_file_name %s"
        % (pprint.pformat(stage_file_name)))
    logger.debug(
        "in resolve_name_to_accessions with stages\n%s"
        % (pprint.pformat(stages.keys())))

    accessions = []
    if not stages:
        return [None]
    for stage_name in [s for s in stages if s]:
        logger.debug("stage name: %s" % (stage_name))
        logger.debug(
            "%d input files:" % (len(stages[stage_name].get('input_files', []))))
        logger.debug(
            "%d output files:" % (len(stages[stage_name].get('output_files', []))))
        if stages[stage_name].get('input_files'):
            all_files = \
                stages[stage_name].get('output_files') + \
                stages[stage_name].get('input_files')
        else:
            all_files = stages[stage_name].get('output_files')
        for stage_file in all_files:
            logger.debug('stage_file[name]: %s' % (stage_file['name']))
            if stage_file['name'] == stage_file_name:
                encode_object = stage_file.get('encode_object')
                if encode_object:
                    if isinstance(encode_object, list):
                        for obj in encode_object:
                            accessions.append(obj.get('accession'))
                    else:
                        accessions.append(encode_object.get('accession'))
    if accessions:
        logger.debug('resolve_name_to_accessons returning:')
        logger.debug('%s' % (pprint.pformat(accessions)))
        return accessions
    else:
        logger.warning('Failed to resolve to accessions, stage_file_name:')
        logger.warning('%s' % (pprint.pformat(stage_file_name)))
        return None


def patch_file(payload, keypair, server, dryrun):
    logger.debug('in patch_file with %s' % (pprint.pformat(payload)))
    accession = payload.pop('accession')
    # old_file_object = common.encoded.get(server + )        
    url = urlparse.urljoin(server, 'files/%s' % (accession))
    if dryrun:
        logger.info(
            "Dry run.  Would PATCH: %s with %s" %
            (accession, pprint.pformat(payload)))
        logger.info("Dry run.  Returning unchanged file object")
        new_file_object = common.encoded_get(
            urlparse.urljoin(server, '/files/%s' % (accession)), keypair)
    else:
        # r = requests.patch(url, auth=keypair, headers={'content-type': 'application/json'}, data=json.dumps(payload))
        logger.info("PATCH file %s" % (url))
        r = common.encoded_patch(url, keypair, payload, return_response=True)
        try:
            r.raise_for_status()
        except:
            logger.error(
                'PATCH file object failed: %s %s' % (r.status_code, r.reason))
            logger.error(r.text)
            new_file_object = None
        else:
            new_file_object = r.json()['@graph'][0]
            logger.info("Patched: %s" % (new_file_object.get('accession')))

    return new_file_object


def post_file(payload, keypair, server, dryrun):
    logger.debug('in post_file with %s' % (pprint.pformat(payload)))
    url = urlparse.urljoin(server, 'files/')
    if dryrun:
        logger.info("Dry run.  Would post: %s" % (pprint.pformat(payload)))
        new_file_object = None
    else:
        # r = requests.post(url, auth=keypair, headers={'content-type': 'application/json'}, data=json.dumps(payload))
        logger.info('POST new file to %s' % (url))
        r = common.encoded_post(url, keypair, payload, return_response=True)
        try:
            r.raise_for_status()
        except:
            logger.error(
                'POST file object failed: %s %s' % (r.status_code, r.reason))
            logger.error(r.text)
            new_file_object = None
        else:
            new_file_object = r.json()['@graph'][0]
            logger.info(
                "New accession: %s" % (new_file_object.get('accession')))

    return new_file_object


def add_tag(dx_fh, tag):
    try:
        dx_fh.add_tags([tag])
    except dxpy.exceptions.ResourceNotFound as e:
        logger.warning(
            '%s adding tag to %s.  Will try in current project.'
            % (e, dx_fh.name))
        try:
            current_project_fh = dxpy.DXFile(
                dx_fh.get_id(), project=dxpy.PROJECT_CONTEXT_ID)
            current_project_fh.add_tags([tag])
        except dxpy.exceptions.ResourceNotFound as e2:
            # give up
            logger.warning('%s.  Skipping saving tag.' % (e2))
            pass
        except:
            raise
    except:
        raise


def set_property(dx_fh, prop):
    try:
        dx_fh.set_properties(prop)
    except dxpy.exceptions.ResourceNotFound as e:
        logger.warning(
            '%s adding property %s to %s.  Will try in current project.'
            % (e, prop, dx_fh.name))
        try:
            current_project_fh = dxpy.DXFile(
                dx_fh.get_id(), project=dxpy.PROJECT_CONTEXT_ID)
            current_project_fh.set_properties(prop)
        except dxpy.exceptions.ResourceNotFound as e2:
            # give up
            logger.warning('%s.  Skipping saving property.' % (e2))
            pass
        except:
            raise
    except dxpy.exceptions.PermissionDenied as e:
        logger.warning(
            '%s adding property %s to %s.'
            % (e, prop, dx_fh.name))
        raise
    except:
        raise


def qckiller(f, server, keypair):
    QC_OBJECS_TO_KILL = [
        'chipseq_filter_quality_metric',
        'samtools_flagstats_quality_metric',
        'idr_quality_metric',
        'histone_chipseq_quality_metric']

    for object_type in QC_OBJECS_TO_KILL:
        url = \
            server + \
            '/search/?type=%s&quality_metric_of=%s&status!=deleted' % (object_type, f.get('@id'))
        logger.info(
            "%s: Searching for existing %s QC objects"
            % (f.get('accession'), url))
        objs = common.encoded_get(url, keypair)['@graph']
        if not objs:
            logger.info(
                "%s: No %s found"
                % (f.get('accession'), url))
        for o in objs:
            url = server + o.get('@id')
            logger.info(
                "%s: Remove existing qc object %s by setting status:deleted"
                % (f.get('accession'), url))
            logger.info("PATCH qc object %s" % (url))
            common.encoded_patch(url, keypair, {'status': 'deleted'})


def accession_file(f, server, keypair, dryrun, force_patch, force_upload, use_content_md5sum, accessioned_file=None):
    # check for duplication
    # - if it has ENCFF or TSTFF number in it's tag, or
    # - if there exists an accessioned file with the same submitted_file_name
    #   that is not deleted, replaced, revoked and has the same size
    # - then there should be a file with the same md5.  If not, warn of a
    #   mismatch between what's at DNAnexus and ENCODEd.
    # - If same md5, return the existing object.
    # - Next, check if there's already a file with the same md5.  If it's
    #   deleted, replaced, revoked, then remodel it if --force_patch=true,
    # - Else warn and return None
    # download
    # calculate md5 and add to f.md5sum
    # post file and get accession, upload credentials
    # upload to S3
    # remove the local file (to save space)
    # return the ENCODEd file object
    logger.debug(
        'in accession_file with f %s'
        % (pprint.pformat(f['submitted_file_name'])))
    dx_fh = f.pop('dx')
    local_fname = dx_fh.name
    logger.info("Downloading %s" % (local_fname))
    dxpy.download_dxfile(dx_fh.get_id(), local_fname)
    f.update({'md5sum': common.md5(local_fname)})
    f['notes'] = json.dumps(f.get('notes'))

    # check to see if md5 already in the database
    url = \
        server + \
        '/md5:%s?format=json&frame=object' % (f.get('md5sum'))
    r = common.encoded_get(url, keypair, return_response=True)
    try:
        r.raise_for_status()
    except:
        if r.status_code == 404:
            logger.info('No md5 matches %s' % (f.get('md5sum')))
            md5_exists = False
        else:
            logger.error(
                'MD5 duplicate check. GET failed: %s %s'
                % (r.status_code, r.reason))
            logger.error(r.text)
            md5_exists = None
    else:
        md5_exists = r.json()

    if use_content_md5sum and accessioned_file and not md5_exists:
        logger.info('Using file with matching content_md5sum')
        md5_exists = accessioned_file
        using_content_md5sum = True
    else:
        using_content_md5sum = False

    # check if an ENCODE accession number in in the list of tags, as it would
    # be if accessioned by this script or similar scripts
    for tag in dx_fh.tags:
        m = re.findall(r'ENCFF\d{3}\D{3}', tag)  # Deliberately ignore TSTFF accessions
        if m:
            logger.info(
                '%s appears to contain ENCODE accession number in tag %s.'
                % (dx_fh.get_id(), m))
            accession_in_tag = True
            # if not force_patch:
            #   return
        else:
            accession_in_tag = False

    # if the same file has been revoked or replaced that requires wrangler
    # intervention
    if md5_exists:
        existing_file_status = md5_exists['status']
        if existing_file_status in ['revoked', 'replaced']:
            raise AccessioningError(
                "File %s with matching MD5 %s exists but has status %s"
                % (md5_exists.get('accession'), md5_exists.get('md5sum'),
                   md5_exists.get('status')))

        if not (force_patch or force_upload):
            logger.info("accession_file: Returning duplicate file unchanged")
            return md5_exists

        if force_patch or force_upload:
            logger.info(
                "accession_file: MD5 exisits, but force_patch, so patching file metatdata")
            f['accession'] = md5_exists['accession']
            if using_content_md5sum and not force_upload:
                logger.info('Using content_md5sum match and not force_upload '
                            'so patching original md5sum and file_size')
                f.pop('md5sum', None)
                f.pop('file_size', None)
            # if the same file has been deleted then we "undelete" it by
            # resetting its status to uploading
            if existing_file_status in ['deleted'] or (existing_file_status in ['upload failed'] and force_upload):
                logger.info(
                    "File %s with matching MD5 %s will be reset to status uploading"
                    % (md5_exists.get('accession'), md5_exists.get('md5sum')))
                f['status'] = 'uploading'
            # blow away existing qc metrics from the existing file object
            qckiller(md5_exists, server, keypair)
            new_file_object = patch_file(f, keypair, server, dryrun)

        if force_upload:
            force_statuses = ["uploading", "upload failed"]
            if new_file_object['status'] not in force_statuses:
                logger.warning(
                    '%s: status is %s, so force_upload is not allowed. Skipping.'
                    % (new_file_object.get('accession'),
                       new_file_object.get('status')))
            else:
                logger.info(
                    "accession_file: MD5 exisits, but force_upload, so uploading and patching file metatdata")
                return_code = common.s3_cp(
                    new_file_object, local_fname, server, keypair)
                if not return_code:
                    logger.info('Upload succeeded')
                logger.debug('s3_cp returned %s' % (return_code))
                assert not return_code, '%s: s3_cp failed. Returned non-zero return code %s' % (new_file_object.get('accession'), return_code)
                add_tag(dx_fh, new_file_object.get('accession'))
    else:
        logger.info("accession_file: New MD5")
        logger.info('posting new file %s' % (f.get('submitted_file_name')))
        logger.debug('%s' % (f))
        new_file_object = post_file(f, keypair, server, dryrun)
        return_code = common.s3_cp(
            new_file_object, local_fname, server, keypair)
        logger.debug('s3_cp returned %s' % (return_code))
        assert not return_code, '%s: s3_cp failed. Returned non-zero return code %s' % (new_file_object.get('accession'), return_code)
        add_tag(dx_fh, new_file_object.get('accession'))
    try:
        os.remove(local_fname)
    except:
        pass

    return new_file_object


def accession_analysis_step_run(analysis_step_run_metadata, keypair, server,
                                dryrun, force_patch, force_upload, use_content_md5sum):
    logger.debug(
        "in accession_analysis_step_run with analysis_step_run_metadata\n%s"
        % (pprint.pformat(analysis_step_run_metadata)))

    if dryrun:
        logger.info("Dry run.  Would POST %s" % (analysis_step_run_metadata))
        return {}

    post_url = urlparse.urljoin(server, '/analysis-step-runs/')
    post_response = common.encoded_post(
        post_url, keypair, analysis_step_run_metadata, return_response=True)
    try:
        post_response.raise_for_status()
    except:  # POST failed
        if post_response.status_code == 409:
            existing_url = urlparse.urljoin(
                server,
                "/%s" % (analysis_step_run_metadata['aliases'][0]))  # assumes there's only one alias
            if force_patch:
                try:
                    put_response = common.encoded_put(
                        existing_url, keypair, analysis_step_run_metadata, return_response=True)
                except:  # PUT failed
                    logger.error(
                        'PUT analysis_step_run object failed: %s %s'
                        % (put_response.status_code, put_response.reason))
                    logger.error(put_response.text)
                    return None
                else:  # PUT succeeded
                    return_object = put_response.json()['@graph'][0]
            else:  # POST failed but not force_patch
                logger.warning(post_response.text)
                return_object = common.encoded_get(existing_url, keypair)
                logger.info(
                    'Using existing analysis_step_run object %s'
                    % (return_object.get('@id')))
        else:  # POST failed with something other than 409
            logger.error(
                'POST analysis_step_run object failed: %s %s'
                % (post_response.status_code, post_response.reason))
            logger.error(post_response.text)
            return_object = None
    else:  # POST succeeded
        logger.info('POST new analysis_step_run to %s' % (post_url))
        return_object = post_response.json()['@graph'][0]
        logger.info(
            "New analysis_step_run uuid: %s" % (return_object.get('uuid')))
    if return_object is None:
        raise AccessioningError(
            "accession_analysis_step_run failed to accession:\n%s"
            % (pprint.pformat(analysis_step_run_metadata)))
    return return_object


def encode_file(keypair, server, field, value):
    logger.info('Searching ENCODE for file with %s=%s' % (field, value))
    search_result = common.encoded_get(
        server + '/search/?type=File&%s=%s' % (field, value),
        keypair=keypair)
    filtered_result = [r for r in search_result.get('@graph', [])
                       if r.get('status') != 'replaced']
    if filtered_result:
        try:
            # Preferentially return relased files.
            released = [f for f in filtered_result
                        if f.get('status') == 'released']
            assert len(released) <= 1, 'More than one released file with {}={} found on portal'.format(
                field, value
            )
            return released[0]
        except IndexError:
            return filtered_result[0]
    else:
        None


def dx_file_at_encode(dx_fh, keypair, server, use_content_md5sum):
    match = encode_file(keypair, server, field='md5sum', value=dxf_md5(dx_fh))
    if not match and use_content_md5sum:
        match = encode_file(
            keypair,
            server,
            field='content_md5sum',
            value=dxf_content_md5(dx_fh)
        )
    return match


def accessioned_outputs(stages, keypair, server, use_content_md5sum):
    files = []
    for (stage_name, outputs) in stages.iteritems():
        stage_metadata = outputs['stage_metadata']
        for i, file_metadata in enumerate(outputs['output_files']):
            file_id = stage_metadata['output'][file_metadata['name']]
            project = stage_metadata['project']
            logger.debug(
                'in accessioned_outputs getting handler for file %s in %s'
                % (file_id, project))
            dx = dxpy.DXFile(file_id, project=project)
            accessioned_file = dx_file_at_encode(dx, keypair, server, use_content_md5sum)
            if accessioned_file:
                logger.info(
                    "Found dx file %s named %s accessioned at ENCODE as %s"
                    % (file_id, dx.name, accessioned_file.get('accession')))
                stages[stage_name]['output_files'][i].update(
                    {'encode_object': accessioned_file})
                files.append(accessioned_file)
    return files


def accession_outputs(stages, keypair, server,
                      dryrun, force_patch, force_upload, use_content_md5sum):
    files = []
    for (stage_name, outputs) in stages.iteritems():
        stage_metadata = outputs['stage_metadata']
        for i, file_metadata in enumerate(outputs['output_files']):
            file_id = stage_metadata['output'][file_metadata['name']]
            project = stage_metadata['project']
            logger.debug(
                'in accession_outputs getting handler for file %s in %s'
                % (file_id, project))
            dx = dxpy.DXFile(file_id, project=project)
            accessioned_file = dx_file_at_encode(dx, keypair, server, use_content_md5sum)

            if accessioned_file:
                logger.info(
                    "Found dx file %s named %s already accessioned at ENCODE as %s"
                    % (file_id, dx.name, accessioned_file.get('accession')))
            if accessioned_file and not force_patch and not force_upload:
                stages[stage_name]['output_files'][i].update(
                    {'encode_object': accessioned_file})
                files.append(accessioned_file)
            else:
                if accessioned_file:
                    if force_patch:
                        logger.info("File already accessioned, but force_patch so patching new metadata")
                    if force_upload:
                        logger.info("File already accessioned, but force_upload so patching new metadata")
                logger.info(
                    "Accessioning dx file %s named %s"
                    % (file_id, dx.name))
                analysis = stage_metadata['parentAnalysis']
                dataset_accession = get_experiment_accession(analysis)
                dx_desc = dx.describe()
                surfaced_outputs = \
                    [o for o in outputs['qc'] if isinstance(o, str)]  # this will be a list of strings
                calculated_outputs = \
                    [o for o in outputs['qc'] if not isinstance(o, str)]  # this will be a list of functions/methods
                logger.debug(
                    'in accession_outputs with stage metadata\n%s'
                    % (stage_metadata.get('name')))
                logger.debug(
                    'in accession_ouputs with surfaced_ouputs %s and calculated_outputs %s from project %s'
                    % (surfaced_outputs, calculated_outputs, project))
                notes_qc = dict(zip(
                    surfaced_outputs,
                    [stage_metadata['output'][metric]
                     for metric in surfaced_outputs]))
                notes_qc.update(dict(zip(
                    [f.__name__ for f in calculated_outputs],
                    [f(stages) for f in calculated_outputs])))
                post_metadata = {
                    'dx': dx,
                    'notes': {
                        'dx-id': dx.get_id(),
                        'dx-createdBy': dx_desc.get('createdBy'),
                        'dx-parentAnalysis': analysis,
                        'qc': notes_qc
                    },
                    # 'aliases': ['ENCODE:%s-%s' %(experiment.get('accession'), static_metadata.pop('name'))],
                    'dataset': dataset_accession,
                    'file_size': dx_desc.get('size'),
                    'submitted_file_name':
                        dx.get_proj_id() + ':' + '/'.join([dx.folder, dx.name])}
                post_metadata.update(file_metadata['metadata'])
                new_file = accession_file(
                    post_metadata, server, keypair,
                    dryrun, force_patch, force_upload,
                    use_content_md5sum, accessioned_file
                )
                stages[stage_name]['output_files'][i].update(
                    {'encode_object': new_file})
                files.append(new_file)
    return files


def new_metadata(old_obj, new_obj):
    logger.debug(
        'in new_metdata with old_object\n%s\nnew_object\n%s',
        pprint.pformat(old_obj),
        pprint.pformat(new_obj))
    for key in new_obj:
        if key not in old_obj:
            return True
        elif key == 'derived_from':
            logger.debug("%s" % (set([re.search("(ENC|TST)FF.{6}", s).group(0) for s in old_obj[key]])))
            logger.debug("%s" % (set([re.search("(ENC|TST)FF.{6}", s).group(0) for s in new_obj[key]])))
            try:
                if set([re.search("(ENC|TST)FF.{6}", s).group(0) for s in old_obj[key]]) \
                   != \
                   set([re.search("(ENC|TST)FF.{6}", s).group(0) for s in new_obj[key]]):
                    return True
            except:
                logger.warning(
                    "In new_metadata: failed to compare derived_from properties.  Old: %s, New: %s"
                    % (old_obj[key], new_obj))
                return True
        elif isinstance(new_obj[key], list):
            if set(new_obj[key]) != set(old_obj.get(key)):
                return True
        else:
            if new_obj[key] != old_obj[key]:
                return True
    return False


def patch_outputs(stages, keypair, server, dryrun):
    logger.debug('in patch_outputs')
    files = []
    for stage_name in [s for s in stages if s]:
        for n, file_metadata in enumerate(stages[stage_name]['output_files']):
            if file_metadata.get('encode_object'):
                logger.info(
                    'patch outputs stage_name %s n %s file_metadata[name] %s encode_object[accession] %s'
                    % (stage_name, n, file_metadata['name'],
                       file_metadata['encode_object'].get('accession')))
                logger.debug(
                    "encode_object %s"
                    % (pprint.pformat(file_metadata['encode_object']['@id'])))
                logger.debug("patch_outputs file_metadata")
                # logger.debug("%s" %(pprint.pformat(file_metadata)))
                accession = file_metadata['encode_object'].get('accession')
                derived_from_accessions = set()  # no duplicates allowed
                for derived_from in file_metadata['derived_from']:
                    # derived from can be a tuple specifying a name from
                    # different set of stages
                    if isinstance(derived_from, tuple):
                        logger.debug(
                            'different stage file_metadata[derived_from] =')
                        logger.debug('%s' % (pprint.pformat(derived_from[1])))
                        stages_to_use = derived_from[0]
                        name_to_use = derived_from[1]
                    else:
                        logger.debug('same stage file_metadata[derived_from]')
                        logger.debug('%s' % (pprint.pformat(derived_from)))
                        stages_to_use = stages
                        name_to_use = derived_from
                    # May see the same accession twice.  If, for example, a
                    # single control is reused, there
                    # will be two paths back to it (one via rep1 one via rep2)
                    # and so it will come out of this loop twice.
                    encode_accessions = resolve_name_to_accessions(stages_to_use, name_to_use)
                    if not encode_accessions:
                        raise AccessioningError(
                            "Expected but found no accessioned file for %s"
                            % (name_to_use))
                    for acc in encode_accessions:
                        if acc:
                            derived_from_accessions.add(acc)
                logger.debug(
                    'derived_from_accessions = %s'
                    % (pprint.pformat(derived_from_accessions)))
                patch_metadata = {
                    'accession': accession,
                    'derived_from': list(derived_from_accessions)
                }
                logger.debug(
                    'patch_metadata = %s'
                    % (pprint.pformat(patch_metadata)))
                if new_metadata(file_metadata['encode_object'], patch_metadata):
                    patched_file = patch_file(
                        patch_metadata, keypair, server, dryrun)
                    if patched_file:
                        stages[stage_name]['output_files'][n]['encode_object'] = patched_file
                        files.append(patched_file)
                    else:
                        logger.error("%s PATCH failed ... skipping" % (accession))
                else:
                    logger.info('Nothing new to patch to %s' % (accession))
            else:
                logger.warning(
                    '%s,%s: No encode object found ... skipping'
                    % (stage_name, file_metadata['name']))
                continue
    return files


def accession_qc_object(obj_type, obj, keypair, server,
                        dryrun, force_patch, force_upload, use_content_md5sum):

    logger.debug(
        'in accession_qc_object with obj_type %s obj.keys() %s'
        % (obj_type, obj.keys()))

    logger.debug(
        'obj[step_run] %s'
        % (obj.get('step_run')))

    # To avoid duplicating qc objects check for same analysis_step_run, and if
    # there is a qc object of this type for that analysis_step run, then PUT,
    # else POST.
    # Result is that the QC objects from the analysis being accessioned will
    # replace those that already exist.

    # get the objects of this object type already in this step run
    url = urlparse.urljoin(
        server,
        '/search/?type=%s&step_run=%s&datastore=database'
        % (obj_type, obj.get('step_run')))

    logger.debug(
        'get qc objects url %s'
        % (url))
    logger.info('Waiting 2 minutes to check if ENCODE server is indexing')
    time.sleep(120)
    while encode_unready(server):
        logger.info(
            'In accession_qc_object but ENCODE server is not ready.  Checking again in 60s.')
        time.sleep(60)

    r = common.encoded_get(url, keypair)
    # need the process_stage special case for the samtools_flagstat objects,
    # which are diffrentiated by processing_stage
    # Also obsolete an older object without processing_stage with a new one
    # with processing_stage
    # Problem:  if the object was just modified, the search above will not
    # pick up the new object.
    existing_objects = \
        [o for o in r['@graph']
         if o['status'] not in DEPRECATED and
         ((o.get('processing_stage') == obj.get('processing_stage')) or
          (obj.get('processing_stage') and o.get('processing_stage') is None))]

    logger.debug(
        'found %d qc objects of type %s'
        % (len(existing_objects), obj_type))

    # pick one to replace
    if existing_objects:
        object_to_replace = existing_objects.pop()
    else:
        object_to_replace = None
    # and delete the rest
    for object_to_delete in existing_objects:
        logger.debug('object_to_delete %s' % (object_to_delete))
        logger.debug(
            'new processing_stage %s old processing stage %s'
            % (obj.get('processing_stage'),
                object_to_delete.get('processing_stage')))
        logger.info(
            'Deleting obsolete qc metric object %s'
            % (object_to_delete['@id']))
        url = urlparse.urljoin(server, object_to_delete['@id'])
        logger.info("PATCH qc object %s to status: 'deleted'" % (url))
        common.encoded_patch(url, keypair, {'status': 'deleted'})
        existing_objects.remove(object_to_delete)

    if object_to_replace:
        # retain any existing links from this metric to existing files
        logger.debug('object_to_replace.quality_metric_of: %s' % (object_to_replace.get('quality_metric_of')))
        logger.debug('obj.quality_metric_of: %s' % (obj.get('quality_metric_of')))
        logger.debug('combined: %s' % (
            (object_to_replace.get('quality_metric_of') or []) +
            (obj.get('quality_metric_of') or [])))
        obj['quality_metric_of'] = list(set(
            [uri.split('/')[2] if '/' in uri else uri for uri in
                (object_to_replace.get('quality_metric_of') or []) +
                (obj.get('quality_metric_of') or [])]
            ))
        url = urlparse.urljoin(server, object_to_replace['@id'])
        logger.info('PUT to %s' % (url))
        logger.debug('PUT %s with %s' % (url, json.dumps(obj)))
        r = common.encoded_put(url, keypair, obj, return_response=True)
    else:
        url = urlparse.urljoin(server, '/%s/' % (obj_type))
        logger.info('POST qc object to %s' % (url))
        logger.debug('POST to %s with %s' % (url, json.dumps(obj)))
        r = common.encoded_post(url, keypair, obj, return_response=True)
    try:
        r.raise_for_status()
    except:
        logger.error('PUT or POST failed: %s %s' % (r.status_code, r.reason))
        logger.error('url was %s' % (url))
        logger.error(r.text)
        new_qc_object = None
    else:
        new_qc_object = r.json()['@graph'][0]

    return new_qc_object


def accession_pipeline(analysis_step_versions, keypair, server,
                       dryrun, force_patch, force_upload, use_content_md5sum):
    patched_files = []
    for (analysis_step_version_id, steps) in analysis_step_versions.iteritems():
        for step in steps:
            if not (step['stages']
                    and step['stage_name']
                    and step['file_names']):
                logger.warning(
                    '%s missing stage metadata (files or stage_name) ... skipping'
                    % (analysis_step_version_id))
                continue
            stage_name = step['stage_name']
            logger.info(stage_name)
            logger.debug(step.keys())
            logger.debug(step['stages'].keys())
            jobid = step['stages'][stage_name]['stage_metadata']['id']
            # analysis_step_version = \
            #     'versionof:%s' % (analysis_step_version_name)
            # versionof namespace has been removed, now need to find versions
            # by uuid
            alias = 'dnanexus:%s' % (jobid)
            if step.pop('virtual', None):
                alias += '-virtual-file-conversion-step'
            analysis_step_run_metadata = {
                'aliases': [alias],
                'analysis_step_version':
                    '/analysis-step-versions/%s/' % (analysis_step_version_id),
                # this used to be taken from the step definition,
                # but the Portal was changed to mean something different
                # so now steps are always set to released
                'status': 'released',
                'dx_applet_details': [{
                    'dx_status': 'finished',
                    'dx_job_id': 'dnanexus:%s' % (jobid),
                }]
            }
            analysis_step_run = accession_analysis_step_run(
                analysis_step_run_metadata, keypair, server,
                dryrun, force_patch, force_upload, use_content_md5sum)

            logger.debug(
                'in accession_pipeline analysis_step_run %s'
                % (pprint.pformat(analysis_step_run)))
            for qc in step['qc_objects']:
                qc_object_name, files_to_associate = next(qc.iteritems())
                qc_objects = \
                    globals()[qc_object_name](
                        analysis_step_run.get('@id'),
                        step['stages'],
                        files_to_associate)
                for qc_object in qc_objects:
                    new_object = accession_qc_object(
                        qc_object_name,
                        qc_object,
                        keypair,
                        server,
                        dryrun,
                        force_patch, force_upload, use_content_md5sum)
                    logger.info(
                        'New %s qc object %s aliases %s'
                        % (qc_object_name, new_object.get('uuid'),
                           new_object.get('aliases')))
                    logger.debug('%s' % (pprint.pformat(new_object)))

            for file_name in step['file_names']:
                for file_accession in resolve_name_to_accessions(step['stages'], file_name):
                    patch_metadata = {
                        'accession': file_accession,
                        'step_run': analysis_step_run.get('@id')
                    }
                    patched_file = \
                        patch_file(patch_metadata, keypair, server, dryrun)
                    patched_files.append(patched_file)

    return patched_files


def accession_mapping_analysis_files(
        mapping_analysis, keypair, server, dryrun, force_patch, force_upload,
        fqcheck, accession_raw, pipeline_version, use_content_md5sum):

    experiment_accession = get_experiment_accession(mapping_analysis)
    if not experiment_accession:
        logger.info(
            "Missing experiment accession or rep in %s, skipping."
            % (mapping_analysis['name']))
        return []

    m = re.match(
        '^Map (ENCSR[0-9]{3}[A-Z]{3}) rep(\d+)', mapping_analysis['name'])
    if m:
        repn = int(m.group(2))
    else:
        logger.error(
            "Missing rep in %s, skipping." % (mapping_analysis['name']))
        return []

    logger.info(
        "%s rep %d: accessioning mapping." % (experiment_accession, repn))

    # experiment = common.encoded_get(
    #     urlparse.urljoin(
    #         server, '/experiments/%s' % (experiment_accession)), keypair)

    analysis_stages = \
        [stage['execution'] for stage in mapping_analysis.get('stages')]
    mapping_stages = get_mapping_stages(
        mapping_analysis, keypair, server, fqcheck, repn)
    if not mapping_stages:
        logger.error('%s: failed to find mapping stages for rep%d'
                     % (mapping_analysis['id'], repn))
        return []

    if accession_raw:
        raw_mapping_stages = get_raw_mapping_stages(
            mapping_analysis, keypair, server, fqcheck, repn)
        if not raw_mapping_stages:
            logger.error('%s: failed to find raw mapping stages for rep%d'
                         % (mapping_analysis['id'], repn))
            return []
    else:
        raw_mapping_stages = None

    for stages in [i for i in [raw_mapping_stages, mapping_stages] if i]:
        output_files = \
            accession_outputs(stages, keypair, server, dryrun,
                              force_patch, force_upload, use_content_md5sum)
        if not output_files:
            logger.error(
                'in accession_mapping_analysis_files, accession_outputs failed')
        files_with_derived = patch_outputs(
            stages, keypair, server, dryrun)

    scrubbed = any([scrubbed_stage(stage) for stage in analysis_stages])
    filtered_bam = \
        'scrubbed_filtered_bam' if scrubbed else 'filtered_bam'
    unfiltered_bam = \
        'scrubbed_unfiltered_bam' if scrubbed else 'mapped_reads'
    mapping_analysis_step_versions = {
        STEP_VERSION_ALIASES[pipeline_version]['bwa-indexing-step']: [
            {
                'stages': "",
                'stage_name': "",
                'file_names': [],
                'status': 'released',
                'qc_objects': []
            }
        ],
        STEP_VERSION_ALIASES[pipeline_version]['bwa-alignment-step']: [
            {
                'stages': mapping_stages,
                'stage_name': get_stage_name(
                    'Filter and QC.*', analysis_stages),
                'file_names': [filtered_bam],
                'status': 'released',
                'qc_objects': [
                    {'chipseq_filter_quality_metric': [filtered_bam]},
                    {'samtools_flagstats_quality_metric': [filtered_bam]}
                ]
            }
        ]
    }

    if accession_raw:
        mapping_analysis_step_versions[STEP_VERSION_ALIASES[pipeline_version]['bwa-alignment-step']].append(
            {
                'stages': raw_mapping_stages,
                'stage_name': get_stage_name('Filter and QC.*' if scrubbed else 'Map ENCSR.*', analysis_stages),
                'file_names': [unfiltered_bam],
                'status': 'released',
                'qc_objects': [
                    {'samtools_flagstats_quality_metric': [unfiltered_bam]}
                ]
            }
        )

    patched_files = accession_pipeline(
        mapping_analysis_step_versions, keypair, server,
        dryrun, force_patch, force_upload, use_content_md5sum)
    return patched_files


def accession_raw_mapping_analysis_files(
        mapping_analysis, keypair, server, dryrun, force_patch, force_upload,
        fqcheck, pipeline_version, use_content_md5sum):

    experiment_accession = get_experiment_accession(mapping_analysis)
    if not experiment_accession:
        logger.info(
            "Missing experiment accession or rep in %s, skipping."
            % (mapping_analysis['name']))
        return []

    m = re.match(
        '^Map (ENCSR[0-9]{3}[A-Z]{3}) rep(\d+)', mapping_analysis['name'])
    if m:
        repn = int(m.group(2))
    else:
        logger.error(
            "Missing rep in %s, skipping." % (mapping_analysis['name']))
        return []

    logger.info(
        "%s rep %d: accessioning mapping." % (experiment_accession, repn))

    # experiment = common.encoded_get(urlparse.urljoin(server,'/experiments/%s' %(experiment_accession)), keypair)

    analysis_stages = \
        [stage['execution'] for stage in mapping_analysis.get('stages')]
    raw_mapping_stages = get_raw_mapping_stages(
        mapping_analysis, keypair, server, fqcheck, repn)

    scrubbed = any([scrubbed_stage(stage) for stage in analysis_stages])
    unfiltered_bam = \
        'scrubbed_unfiltered_bam' if scrubbed else 'mapped_reads'

    output_files = \
        accession_outputs(raw_mapping_stages, keypair, server, dryrun,
                          force_patch, force_upload, use_content_md5sum)
    files_with_derived = patch_outputs(raw_mapping_stages, keypair, server, dryrun)

    raw_mapping_analysis_step_versions = {
        STEP_VERSION_ALIASES[pipeline_version]['bwa-indexing-step']: [
            {
                'stages': "",
                'stage_name': "",
                'file_names': [],
                'status': 'released',
                'qc_objects': []
            }
        ],
        STEP_VERSION_ALIASES[pipeline_version]['bwa-raw-alignment-step']: [
            {
                'stages': raw_mapping_stages,
                'stage_name': get_stage_name('Filter and QC.*' if scrubbed else 'Map ENCSR.*', analysis_stages),
                'file_names': [unfiltered_bam],
                'status': 'released',
                'qc_objects': [
                    {'samtools_flagstats_quality_metric': [unfiltered_bam]}
                ]
            }
        ]
    }

    patched_files = accession_pipeline(
        raw_mapping_analysis_step_versions, keypair, server,
        dryrun, force_patch, force_upload, use_content_md5sum)
    return patched_files


def filtered_bam_output_name(mapping_stages):
    stage_name = next(name
                      for name in mapping_stages.keys()
                      if name.startswith('Filter and QC'))
    filter_stage_metadata = \
        mapping_stages.get(stage_name).get('stage_metadata')
    if 'scrubbed_filtered_bam' in filter_stage_metadata.get('output'):
        return 'scrubbed_filtered_bam'
    elif 'filtered_bam' in filter_stage_metadata.get('output'):
        return 'filtered_bam'
    else:
        return None


def accession_histone_analysis_files(peaks_analysis, keypair, server, dryrun,
                                     force_patch, force_upload, fqcheck,
                                     skip_control, pipeline_version, use_content_md5sum):

    experiment_accession = get_experiment_accession(peaks_analysis)
    if experiment_accession:
        logger.info('%s: accession histone peaks' % (experiment_accession))
    else:
        logger.error(
            "No experiment accession in %s, skipping."
            % (peaks_analysis['executableName']))
        return None
    experiment = common.encoded_get(
        urlparse.urljoin(server, '/experiments/%s' % (experiment_accession)),
        keypair)
    logger.debug('got experiment %s' % (experiment.get('accession')))
    unreplicated_analysis = is_unreplicated_analysis(peaks_analysis)

    # returns a list with two elements:  the mapping stages for [rep1,rep2]
    # in this context rep1,rep2 are the first and second replicates in the
    # pipeline.  They may have been accessioned on the portal with any
    # arbitrary biological_replicate_numbers.
    mapping_stages = \
        get_peak_mapping_stages(peaks_analysis, keypair, server, fqcheck)
    if not mapping_stages:
        logger.error("Failed to find peak mapping stages")
        return None
    # returns a list with three elements: the mapping stages for the controls
    # for [rep1, rep2, pooled], the control stages for rep1 and rep2 might be
    # the same as the pool if the experiment used pooled controls
    if skip_control:
        control_stages = []
        logger.info("skip_control, so ignoring control mapping stages")
    else:
        control_stages = get_control_mapping_stages(
            peaks_analysis, keypair, server, fqcheck)
        if not control_stages:
            logger.error("Failed to find control mapping stages")
            return None

    # returns the stages for peak calling
    peak_stages = get_histone_peak_stages(
        peaks_analysis,
        mapping_stages,
        control_stages,
        experiment,
        keypair,
        server)
    if not peak_stages:
        logger.error("Failed to find peak stages")
        return None

    # accession all the output files
    output_files = []
    # retrieve file metadata from ENCODE Portal for
    # exected-to-be-already-accessioned mapping files
    # fail if they are not accessioned
    for stages in control_stages:
        if stages:
            logger.info('Retrieving accessioned outputs for control mappings')
            output_files.extend(accessioned_outputs(
                stages, keypair, server, use_content_md5sum))
    for stages in mapping_stages:
        if stages:
            logger.info('Retrieving accessioned outputs for experiment mappings')
            output_files.extend(accessioned_outputs(
                stages, keypair, server, use_content_md5sum))

    # accession all the output files
    for stages in peak_stages:
        if stages:
            logger.info('accessioning output')
            output_files.extend(accession_outputs(
                stages, keypair, server, dryrun, force_patch, force_upload, use_content_md5sum))

    # now that we have file accessions, loop again and patch derived_from
    files_with_derived = []
    for stages in peak_stages:
        if stages:
            files_with_derived.extend(
                patch_outputs(stages, keypair, server, dryrun))

    replicated_peaks_filename = 'overlapping_peaks'
    replicated_peaks_bb_filename = 'overlapping_peaks_bb'

    full_analysis_step_versions = {

        STEP_VERSION_ALIASES[pipeline_version][
            'histone-unreplicated-peak-calling-step'
            if unreplicated_analysis else
            'histone-peak-calling-step']: [
            {
                'stages': peak_stage,
                'stage_name': 'ENCODE Peaks',
                'file_names':
                    [
                     'rep1_fc_signal', 'rep1_pvalue_signal', 'rep1_narrowpeaks'
                    ] if unreplicated_analysis else [
                     'rep1_fc_signal', 'rep2_fc_signal', 'pooled_fc_signal',
                     'rep1_pvalue_signal', 'rep2_pvalue_signal',
                     'pooled_pvalue_signal', 'rep1_narrowpeaks',
                     'rep2_narrowpeaks', 'pooled_narrowpeaks'],
                'status': 'released',
                'qc_objects': []
            } for peak_stage in peak_stages
        ],
        STEP_VERSION_ALIASES[pipeline_version][
            'histone-unreplicated-partition-concordance-step'
            if unreplicated_analysis else
            'histone-overlap-peaks-step']: [
            {
                'stages': peak_stage,
                'stage_name': next(
                    stage_name
                    for stage_name in peak_stage.keys()
                    if re.match('(Overlap|Final) narrowpeaks', stage_name)),
                'file_names': [replicated_peaks_filename],
                'status': 'released',
                'qc_objects': [{'histone_chipseq_quality_metric': [replicated_peaks_filename]}]
            } for peak_stage in peak_stages
        ],
        STEP_VERSION_ALIASES[pipeline_version][
            'histone-unreplicated-peaks-to-bigbed-step'
            if unreplicated_analysis else
            'histone-peaks-to-bigbed-step']: [
            {
                'stages': peak_stage,
                'stage_name': 'ENCODE Peaks',
                'file_names':
                    [
                     'rep1_narrowpeaks_bb'
                    ] if unreplicated_analysis else [
                     'rep1_narrowpeaks_bb', 'rep2_narrowpeaks_bb',
                     'pooled_narrowpeaks_bb'],
                'status': 'released',
                'virtual': True,
                'qc_objects': []
            } for peak_stage in peak_stages
        ],
        STEP_VERSION_ALIASES[pipeline_version][
            'histone-unreplicated-partition-concordance-peaks-to-bigbed-step'
            if unreplicated_analysis else
            'histone-replicated-peaks-to-bigbed-step']: [
            {
                'stages': peak_stage,
                'stage_name': next(
                    stage_name
                    for stage_name in peak_stage.keys()
                    if re.match('(Overlap|Final) narrowpeaks', stage_name)),
                'file_names': [replicated_peaks_bb_filename],
                'status': 'released',
                'virtual': True,
                'qc_objects': [{'histone_chipseq_quality_metric': [replicated_peaks_bb_filename]}]
            } for peak_stage in peak_stages
        ]
    }

    patched_files = accession_pipeline(
        full_analysis_step_versions, keypair, server,
        dryrun, force_patch, force_upload, use_content_md5sum)
    return patched_files


def stage_output_names(stages, stage_name):
    output_names = []
    for stage in stages:
        for output_file in stage[stage_name]['output_files']:
            output_names.append(output_file['name'])
    return output_names


def accession_tf_analysis_files(peaks_analysis, keypair, server, dryrun,
                                force_patch, force_upload, fqcheck,
                                signal_only, skip_control, pipeline_version, use_content_md5sum):

    experiment_accession = get_experiment_accession(peaks_analysis)
    if experiment_accession:
        logger.info('%s: accession TF peaks' % (experiment_accession))
    else:
        logger.error(
            "No experiment accession in %s, skipping."
            % (peaks_analysis['executableName']))
        return None
    experiment = common.encoded_get(
        urlparse.urljoin(server, '/experiments/%s' % (experiment_accession)),
        keypair)
    logger.debug('got experiment %s' % (experiment.get('accession')))
    unreplicated_analysis = is_unreplicated_analysis(peaks_analysis)

    # returns a list with elements:  the mapping stages for [rep1,...,repn]
    # in this context rep1,rep2 are the first and second replicates in the
    # pipeline.  They may have been accessioned on the portal with any
    # arbitrary biological_replicate_numbers.
    mapping_stages = \
        get_peak_mapping_stages(peaks_analysis, keypair, server, fqcheck)
    if not mapping_stages:
        logger.error("Failed to find peak mapping stages")
        return None

    # returns a list with three elements: the mapping stages for the controls
    # for [rep1, rep2, pooled], the control stages for rep1 and rep2 might be
    # the same as the pool if the experiment used pooled controls
    if skip_control:
        control_stages = []
        logger.info("skip_control, so ignoring control mapping stages")
    else:
        control_stages = get_control_mapping_stages(
            peaks_analysis, keypair, server, fqcheck)
        if not control_stages:
            logger.error("Failed to find control mapping stages")
            return None

    # returns the stages for peak calling
    peak_stages = get_tf_peak_stages(
        peaks_analysis, mapping_stages, control_stages,
        experiment, keypair, server, signal_only)
    if not peak_stages:
        logger.error("Failed to find peak stages")
        return None

    output_files = []

    # retrieve file metadata from ENCODE Portal for
    # exected-to-be-already-accessioned mapping files
    # fail if they are not accessioned
    for stages in control_stages:
        if stages:
            logger.info('Retrieving accessioned outputs for control mappings')
            output_files.extend(accessioned_outputs(
                stages, keypair, server, use_content_md5sum))
    for stages in mapping_stages:
        if stages:
            logger.info('Retrieving accessioned outputs for experiment mappings')
            output_files.extend(accessioned_outputs(
                stages, keypair, server, use_content_md5sum))

    # accession all the output files
    for stages in peak_stages:
        if stages:
            logger.info('accessioning output')
            output_files.extend(accession_outputs(
                stages, keypair, server, dryrun, force_patch, force_upload, use_content_md5sum))

    # now that we have file accessions, loop again and patch derived_from
    files_with_derived = []
    for stages in peak_stages:
        if stages:
            files_with_derived.extend(
                patch_outputs(stages, keypair, server, dryrun))

    signal_filenames = \
        ['rep1_fc_signal', 'rep1_pvalue_signal'] if unreplicated_analysis else \
        ['rep1_fc_signal', 'rep2_fc_signal', 'pooled_fc_signal',
         'rep1_pvalue_signal', 'rep2_pvalue_signal',
         'pooled_pvalue_signal']
    peaks_filenames = \
        ['rep1_peaks'] if unreplicated_analysis else \
        ['rep1_peaks', 'rep2_peaks', 'pooled_peaks']
    idr_filenames = \
        ['stable_set'] if unreplicated_analysis else \
        [fn for fn in ['optimal_set', 'conservative_set'] if fn in stage_output_names(peak_stages, 'Final IDR peak calls')]
    peaks_bb_filenames = \
        ['rep1_peaks_bb'] if unreplicated_analysis else \
        ['rep1_peaks_bb', 'rep2_peaks_bb', 'pooled_peaks_bb']
    idr_bb_filenames = \
        ['stable_set_bb'] if unreplicated_analysis else \
        [fn for fn in ['optimal_set_bb', 'conservative_set_bb'] if fn in stage_output_names(peak_stages, 'Final IDR peak calls')]

    full_analysis_step_versions = {
       STEP_VERSION_ALIASES[pipeline_version][
            'tf-unreplicated-macs2-signal-calling-step'
            if unreplicated_analysis else
            'tf-macs2-signal-calling-step']: [
            {
                'stages': peak_stage,
                'stage_name': 'ENCODE Peaks',
                'file_names': signal_filenames,
                'status': 'released',
                'qc_objects': []
            } for peak_stage in peak_stages
        ],
    }

    full_analysis_step_versions.update({
        STEP_VERSION_ALIASES[pipeline_version][
            'tf-unreplicated-spp-peak-calling-step'
            if unreplicated_analysis else
            'tf-spp-peak-calling-step']: [
            {
                'stages': peak_stage,
                'stage_name': 'SPP Peaks',
                'file_names': peaks_filenames,
                'status': 'released',
                'qc_objects': []
            } for peak_stage in peak_stages
        ],
        STEP_VERSION_ALIASES[pipeline_version][
            'tf-unreplicated-idr-step'
            if unreplicated_analysis else
            'tf-idr-step']: [
            {
                'stages': peak_stage,
                'stage_name': 'Final IDR peak calls',
                'file_names': idr_filenames,
                'status': 'released',
                'qc_objects': [{'idr_quality_metric': idr_filenames}]
            } for peak_stage in peak_stages
        ],
        STEP_VERSION_ALIASES[pipeline_version][
            'tf-unreplicated-peaks-to-bigbed-step'
            if unreplicated_analysis else
            'tf-peaks-to-bigbed-step']: [
            {
                'stages': peak_stage,
                'stage_name': 'SPP Peaks',
                'file_names': peaks_bb_filenames,
                'status': 'released',
                'virtual': True,
                'qc_objects': []
            } for peak_stage in peak_stages
        ],
        STEP_VERSION_ALIASES[pipeline_version][
            'tf-unreplicated-idr-peaks-to-bigbed-step'
            if unreplicated_analysis else
            'tf-idr-peaks-to-bigbed-step']: [
            {
                'stages': peak_stage,
                'stage_name': 'Final IDR peak calls',
                'file_names': idr_bb_filenames,
                'status': 'released',
                'virtual': True,
                'qc_objects': [{'idr_quality_metric': idr_bb_filenames}]
            } for peak_stage in peak_stages
        ]
    } if not signal_only else {})

    patched_files = accession_pipeline(
        full_analysis_step_versions,
        keypair, server, dryrun, force_patch, force_upload, use_content_md5sum)
    return patched_files


def infer_pipeline(analysis):
    if any(name.startswith('histone_chip_seq') for name in
           [analysis.get('executableName'), analysis.get('name')]):
        return "histone"
    elif any(name.startswith('tf_chip_seq') for name in
             [analysis.get('executableName'), analysis.get('name')]):
        return "tf"
    elif (analysis.get('executableName').startswith('ENCODE mapping pipeline') or
          (any([re.match("Map", stage['name'])
                for stage in analysis['workflow']['stages']]) and
           any([re.match("Filter", stage['name'])
                for stage in analysis['workflow']['stages']]))):
        return "mapping"
    elif (any([re.match("Map", stage['name'])
               for stage in analysis['workflow']['stages']]) and
          not any([re.match("Filter", stage['name'])
                   for stage in analysis['workflow']['stages']])):
        return "raw"
    else:
        return None


def pipeline_version_by_date(analysis):
    analysis_date = analysis.get('created')
    # get the largest version number that was activated on a date before
    # this analysis was created
    pipeline_version = str(max([
        float(version) for version in VERSION_TIMES
        if VERSION_TIMES[version] < analysis_date])) or None
    return pipeline_version


def infer_pipeline_version(analysis):
    try:
        workflow = dxpy.describe(
            analysis['workflow']['id'], fields={'properties': True})
    except dxpy.exceptions.ResourceNotFound:
        pipeline_version = pipeline_version_by_date(analysis)
        logger.warning(
            "Workflow for %s is missing.  Inferred version %s"
            % (analysis.get('id'), pipeline_version))
    else:
        pipeline_version = \
            workflow['properties'].get('pipeline_version') or pipeline_version_by_date(analysis)

    return pipeline_version


@dxpy.entry_point('accession_analysis_id')
def accession_analysis_id(debug, key, keyfile, dryrun, force_patch,
                          force_upload, fqcheck, analysis_id, pipeline,
                          project, accession_raw, signal_only, skip_control, use_content_md5sum):

    if debug:
        logger.setLevel(logging.DEBUG)
        logger.debug('In accession_anlysis_id with logging level DEBUG')
    else:
        logger.setLevel(logging.INFO)
        logger.info('In accession_anlysis_id with logging level INFO')

    # fetch the credentials from the DCC Credentials project
    dxpy.download_folder(
        DCC_CREDENTIALS_PROJECT, '.', folder=DCC_CREDENTIALS_FOLDER)

    authid, authpw, server = common.processkey(key, keyfile)
    keypair = (authid, authpw)

    logger.debug(analysis_id)
    analysis = dxpy.describe(analysis_id.strip())
    experiment = get_experiment_accession(analysis)

    if not pipeline:
        inferred_pipeline = infer_pipeline(analysis)
    else:
        inferred_pipeline = pipeline

    pipeline_version = infer_pipeline_version(analysis)
    logger.info(
        "Accessioning as %s version %s"
        % (inferred_pipeline, pipeline_version))

    output = {
        'analysis': analysis_id,
        'experiment': experiment,
    }
    logger.info(
        'Accessioning %s name %s executableName %s'
        % (analysis.get('id'),
           analysis.get('name'),
           analysis.get('executableName')))

    try:
        if inferred_pipeline == "histone":
            logger.info('accession histone analysis started')
            output.update(
                {'dx_pipeline': 'histone_chip_seq'})
            accessioned_files = \
                accession_histone_analysis_files(
                    analysis, keypair, server, dryrun, force_patch,
                    force_upload, fqcheck, skip_control, pipeline_version,
                    use_content_md5sum)
            logger.info('accession histone analysis completed')
        elif inferred_pipeline == "mapping":
            logger.info('accession mapping analysis started')
            output.update(
                {'dx_pipeline': 'ENCODE mapping pipeline'})
            accessioned_files = \
                accession_mapping_analysis_files(
                    analysis, keypair, server, dryrun, force_patch,
                    force_upload, fqcheck, accession_raw, pipeline_version,
                    use_content_md5sum)
            logger.info('accession mapping analysis completed')
        elif inferred_pipeline == "tf":
            logger.info('accession tf_chip_seq analysis started')
            output.update(
                {'dx_pipeline': 'tf_chip_seq'})
            accessioned_files = \
                accession_tf_analysis_files(
                    analysis, keypair, server, dryrun, force_patch,
                    force_upload, fqcheck, signal_only, skip_control,
                    pipeline_version, use_content_md5sum)
            logger.info('accession tf_chip_seq analysis completed')
        elif inferred_pipeline == "raw":
            logger.info('accession raw mapping analysis started')
            output.update(
                {'dx_pipeline': 'ENCODE raw mapping pipeline'})
            accessioned_files = \
                accession_raw_mapping_analysis_files(
                    analysis, keypair, server, dryrun, force_patch,
                    force_upload, fqcheck, pipeline_version,
                    use_content_md5sum)
            logger.info('accession raw mapping analysis completed')
        else:
            logger.error(
                'unrecognized analysis pattern %s %s ... skipping.'
                % (analysis.get('name'),
                   analysis.get('executableName')))
            output.update(
                {'dx_pipeline': 'unrecognized'})
            accessioned_files = None
    except:
        raise
        # # if only accessioning one job, throw an error
        # # otherwise print the traceback and try to accession the
        # # other jobs
        # if len(ids) == 1:
        #     raise
        # else:
        #     traceback.print_exc()
        #     accessioned_files = None
        #     file_accessions = None

    else:
        file_accessions = \
            [f.get('accession') for f in (accessioned_files or [])]
        if file_accessions:
            url = server+"/experiments/%s" % (experiment)
            experiment_obj = common.encoded_get(
                url, keypair=keypair, frame='page')
            target = experiment_obj.get('target')
            # mapping non-controls is not a complete pipeline run
            if target and ('control' not in target['investigated_as']) and (inferred_pipeline in ['mapping', 'raw']):
                internal_status = 'processing'
            else:  # everything else is assumed to be complete
                internal_status = 'pipeline completed'
            logger.info("PATCH status %s to experiment %s" % (internal_status, url))
            r = common.encoded_patch(
                url,
                keypair,
                {"internal_status": internal_status},
                return_response=True)
            try:
                r.raise_for_status()
            except:
                logger.error(
                    "Tried but failed to update experiment "
                    "internal_status to pipeline completed")
                logger.error(r.text)

    logger.info("Accessioned: %s" % (file_accessions))
    output.update({'files': file_accessions})
    output_row = copy.copy(output)
    output.update({'output_row': output_row})
    return output


@dxpy.entry_point('postprocess')
def postprocess(outfn, output_rows):
    with open(outfn, 'w') as fh:
        fieldnames = [
            'analysis',
            'experiment',
            'dx_pipeline',
            'files',
            'error'
        ]
        output_writer = csv.DictWriter(fh, fieldnames, delimiter='\t')
        output_writer.writeheader()
        for output_row in output_rows:
            output_writer.writerow(output_row)

    for line in open(outfn, 'r'):
        print(line)

    output = {"outfile": dxpy.upload_local_file(outfn)}

    return output


def encode_unready(server):
    url = server + "/_indexer"
    try:
        indexing_report = common.encoded_get(url)
    except Exception as e:
        logger.error(
            '%s: Could not get indexing report from %s' % (e, url))
        return True
    try:
        indexing_status = indexing_report.get('status')
    except:
        logger.error(
            'Could not interpret _indexer report: %s' % (indexing_report))
        return True
    logger.debug('encode_indexing: found status %s' % (indexing_status))
    if not indexing_status:
        logger.error(
            'Could not get indexing status from %s' % (indexing_report))
        return True
    elif indexing_status == 'waiting':
        return False
    elif indexing_status == 'indexing':
        return True
    else:
        return True


@dxpy.entry_point('main')
def main(outfn, debug, dryrun,
         force_patch, force_upload, fqcheck, use_content_md5sum,
         key=None, keyfile=None, pipeline=None, analysis_ids=None, infile=None, project=None,
         accession_raw=False, signal_only=False, skip_control=False,
         wait_on_files=None, encoded_check=True):

    # wait_on_files is never used here, it is just a place-holder input field
    # to block the platform from running accession_analysis until all the
    # previous stages are complete

    if debug:
        logger.setLevel(logging.DEBUG)
        logger.debug('Set logger level to logging.DEBUG')
    else:
        logger.setLevel(logging.INFO)
        logger.info('Set logger level to logging.INFO')

    # fetch the credentials from the DCC Credentials project
    dxpy.download_folder(
        DCC_CREDENTIALS_PROJECT, '.', folder=DCC_CREDENTIALS_FOLDER)

    if infile is not None:
        infile = dxpy.DXFile(infile)
        dxpy.download_dxfile(infile.get_id(), "infile")
        ids = open("infile", 'r')
    elif analysis_ids is not None:
        ids = analysis_ids
    else:
        logger.error(
            "Must supply one of --infile or a list of analysis-ids")
        return

    if not key or key in ['www', 'submit', 'production']:
        key = dxpy.api.system_whoami()['id']
    elif key == 'test':
        key = dxpy.api.system_whoami()['id'] + "-test"

    key_tuple = common.processkey(key, keyfile)
    assert key_tuple, "ERROR: Key %s is not found in the keyfile %s" % (key, keyfile)
    authid, authpw, server = key_tuple
    keypair = (authid, authpw)

    accession_subjobs = []
    for (i, analysis_id) in enumerate(ids):

        if analysis_id == 'self':
            self_analysis_id = dxpy.describe(dxpy.JOB_ID)['analysis']
            analysis_id = self_analysis_id

        try:
            subjob_name = "Accession %s" % ((dxpy.DXAnalysis(analysis_id)).name)
        except:
            logger.warning(
                "Failed to construct subjob name from analysis %s. Using default"
                % (analysis_id))
            subjob_name = "accession_%s" % (analysis_id)

        accession_subjob_input = {
            "debug": debug,
            "key": key,
            "keyfile": keyfile,
            "dryrun": dryrun,
            "force_patch": force_patch,
            "force_upload": force_upload,
            "use_content_md5sum": use_content_md5sum,
            "fqcheck": fqcheck,
            "pipeline": pipeline,
            "analysis_id": analysis_id,
            "project": project,
            "accession_raw": accession_raw,
            "signal_only": signal_only,
            "skip_control": skip_control
        }

        logger.info("Accession job input: %s" % (accession_subjob_input))
        while encode_unready(server):
            if not encoded_check:
                logger.info(
                    'ENCODE server is not ready but encoded_check=False so continuing.')
                break
            else:
                logger.info(
                    'ENCODE server is not ready.  Checking again in 60s.')
                time.sleep(60)
        accession_subjobs.append(
            dxpy.new_dxjob(
                fn_input=accession_subjob_input,
                fn_name="accession_analysis_id",
                name=subjob_name,
                properties={'id': analysis_id}))

    postprocess_subjob = dxpy.new_dxjob(
        fn_input={
            "outfn": outfn,
            "output_rows":
                [accession_subjob.get_output_ref("output_row")
                 for accession_subjob in accession_subjobs]
        },
        fn_name="postprocess",
        depends_on=accession_subjobs,
        name="output report"
    )

    output = {"outfile": postprocess_subjob.get_output_ref("outfile")}
    return output


dxpy.run()
