<gh_stars>10-100
from ChipseqReport import *


class FoldChangeTracker(TrackerSQL):

    '''the fold change tracker ignores unstimulated tracks.'''

    def __init__(self, *args, **kwargs):
        TrackerSQL.__init__(self, *args, backend=DATABASE, **kwargs)

    def getTracks(self, subset=None):
        tracks = TrackerSQL.getTracks(self, subset=subset)
        return [x for x in tracks if TAG_UNSTIM not in x]

##########################################################################
##########################################################################
##########################################################################
# Annotation of bases with SNPs
##########################################################################


class IntervalsSummary(DefaultTracker):

    """Summary stats of intervals called by the peak finder.
    """

    mPattern = "_intervals$"

    def __call__(self, track, slice=None):
        data = self.getFirstRow(
            "SELECT COUNT(*), AVG(length), SUM(nprobes)  FROM %(track)s_intervals" % locals())
        return odict(list(zip(("nintervals", "<length>", "nprobes"), data)))

##########################################################################
##########################################################################
##########################################################################
# Distribution of interval lengths
##########################################################################


class IntervalLengths(DefaultTracker):

    """Distribution of interval sizes.
    """

    mPattern = "_intervals$"

    def __call__(self, track, slice=None):
        data = self.getValues(
            "SELECT length FROM %(track)s_intervals" % locals())
        return {"length": data}

##########################################################################
##########################################################################
##########################################################################
# Distribution of peak values
##########################################################################


class IntervalPeakValues(DefaultTracker):

    """Distribution of peak values (the number of reads at peak).
    """

    mPattern = "_intervals$"

    def __call__(self, track, slice=None):
        data = self.getValues(
            "SELECT peakval FROM %(track)s_intervals" % locals())
        return {"peakval": data}

##########################################################################
##########################################################################
##########################################################################
# Distribution of peak values
##########################################################################


class EnrichmentOverUnstim(DefaultTracker):

    """For each peakval, present the fold enrichment of a track
    compared to the Unstim set.

    Useful for diagnosing cutoffs.
    """

    mPattern = "_intervals$"

    def __call__(self, track, slice=None):

        cellline, condition, replicate = splitTrack(track)
        if condition == "Unstim":
            return []

        if replicate is None:
            replicate = ""
        other_track = "run" + cellline + "Unstim" + replicate

        data_fg = self.getValues(
            "SELECT peakval FROM %(track)s_intervals ORDER BY peakval DESC" % locals())
        data_bg = self.getValues(
            "SELECT peakval FROM %(other_track)s_intervals ORDER BY peakval DESC" % locals())

        mi = min(data_bg + data_fg)
        ma = max(data_bg + data_fg)

        n_fg = len(data_fg)
        n_bg = len(data_bg)

        hist_fg, bin_edges = numpy.histogram(
            data_fg, bins=numpy.arange(mi, ma + 1.0, 1.0))
        hist_bg, bin_edges = numpy.histogram(
            data_bg, bins=numpy.arange(mi, ma + 1.0, 1.0))

        hist_fg = hist_fg[::-1].cumsum()[::-1]
        hist_bg = hist_bg[::-1].cumsum()[::-1]

        fold = []
        for fg, bg in zip(hist_fg, hist_bg):
            fold.append(float(bg) / fg)

        return odict((("peakval", bin_edges[:-1]),
                      ("fold", fold)))

##########################################################################
##########################################################################
##########################################################################
# Distribution of average values
##########################################################################


class IntervalAverageValues(DefaultTracker):

    """Distribution of average values (the average number of reads within the interval)
    """

    mPattern = "_intervals$"

    def __call__(self, track, slice=None):
        data = self.getValues(
            "SELECT avgval FROM %(track)s_intervals" % locals())
        return {"avgval": data}

##########################################################################
##########################################################################
##########################################################################
# Distribution of average values
##########################################################################


class IntervalLengthVsAverageValue(DefaultTracker):

    """Length vs average value.
    """

    mPattern = "_intervals$"

    def __call__(self, track, slice=None):
        data = self.get(
            "SELECT length, avgval FROM %(track)s_intervals" % locals())
        return odict(list(zip(("length", "avgval"), list(zip(*data)))))

##########################################################################
##########################################################################
##########################################################################
# Distribution of average values
##########################################################################


class IntervalLengthVsPeakValue(DefaultTracker):

    """Length vs peak value
    """

    mPattern = "_intervals$"

    def __call__(self, track, slice=None):
        data = self.get(
            "SELECT length, peakval FROM %(track)s_intervals" % locals())
        return odict(list(zip(("length", "peakval"), list(zip(*data)))))

##########################################################################
##########################################################################
##########################################################################
# peak location within interval
##########################################################################


class PeakLocation(DefaultTracker):
    mPattern = "_intervals$"

    def __call__(self, track, slice=None):
        data1 = self.getValues(
            "SELECT (PeakCenter - start) / CAST( Length as FLOAT) - 0.5 FROM %(track)s_intervals" % locals())
        data2 = self.getValues(
            "SELECT (end - PeakCenter) / CAST( Length as FLOAT) - 0.5 FROM %(track)s_intervals" % locals())
        return {"distance": data1 + data2}

##########################################################################
##########################################################################
##########################################################################
# distance of peak to end of interval
##########################################################################


class PeakDistance(DefaultTracker):
    mPattern = "_intervals$"

    def __call__(self, track, slice=None):
        data1 = self.getValues(
            "SELECT PeakCenter - start FROM %(track)s_intervals" % locals())
        data2 = self.getValues(
            "SELECT end - PeakCenter FROM %(track)s_intervals" % locals())
        return {"distance": data1 + data2}

##########################################################################
##########################################################################
##########################################################################
##
##########################################################################


class IntervalList(DefaultTracker):

    '''list of intervals.'''

    nresults = 10
    mColumnsFixed = ("pos", "length")
    mColumnsVariable = ("peakval", "avgval")
    mPattern = "_intervals$"

    def getSQLStatement(self, track, slice=None):

        statement = '''
          SELECT 
               i.interval_id, i.contig, i.start, i.end, i.peakval, i.avgval
          FROM 
               %(track)s_intervals AS i
          ORDER BY i.peakval DESC''' % locals()

        if self.nresults:
            statement += " LIMIT %i" % self.nresults

        return statement

    def __call__(self, track, slice=None):

        statement = self.getSQLStatement(track, slice)
        data = self.get(statement)
        ucsc_genome = UCSC_GENOME
        n = odict()
        for d in data:
            id, contig, start, end = d[:4]
            pos = "`%(contig)s:%(start)i..%(end)i <http://genome.ucsc.edu/cgi-bin/hgTracks?db=%(ucsc_genome)s&position=%(contig)s:%(start)i..%(end)i>`_" \
                % locals()
            n[str(id)] = odict(
                list(zip(self.mColumnsFixed + self.mColumnsVariable, (pos, end - start,) + d[4:])))

        return n

##########################################################################
##########################################################################
##########################################################################
##
##########################################################################


class IntervalListFull(DefaultTracker):

    '''list of all intervals.

    Table for export.
    '''

    nresults = None
    mPattern = "_intervals$"

    def __call__(self, track, slice=None):

        statement = '''
          SELECT 
               i.contig, i.start, i.end, i.peakval, i.avgval
          FROM 
               %(track)s_intervals AS i
          ORDER BY i.peakval DESC''' % locals()

        data = self.get(statement)
        return odict(list(zip(
            ("contig", "start", "end", "peakval", "avgval"),
            list(zip(*data)))))

##########################################################################
##########################################################################
##########################################################################
##
##########################################################################


class IntervalListPeakval(IntervalList):

    '''list of intervals.'''

    def getSQLStatement(self, track, slice=None):

        nresults = self.nresults

        statement = '''
          SELECT 
               i.interval_id, i.contig, i.start, i.end, i.peakval, i.avgval, i.length
          FROM 
               %(track)s_intervals AS i
          ORDER BY i.peakval DESC
          LIMIT %(nresults)s''' % locals()

        return statement

##########################################################################
##########################################################################
##########################################################################
##
##########################################################################


class IntervalListFoldChange(FoldChangeTracker, IntervalList):

    '''list of intervals.'''

    mColumnsVariable = ("fg_anysense_mean", "bg_anysense_mean", "fold_mean",
                        "fg_anysense_max", "bg_anysense_max", "fold_max")
    mMinFold = 1.5
    mPattern = "_readcounts$"

    def getSQLStatement(self, track, slice=None):
        nresults = self.nresults
        minfold = self.mMinFold

        statement = '''
          SELECT 
               i.interval_id, i.contig, i.start, i.end, 
               fg_anysense_mean, bg_anysense_mean,
               fg_anysense_mean / bg_anysense_mean AS fold_mean,
               fg_anysense_max, bg_anysense_max,
               cast(fg_anysense_max as float)/ bg_anysense_max AS fold_max
          FROM 
               %(track)s_intervals AS i,
               %(track)s_readcounts AS c
          WHERE
               i.interval_id = c.gene_id AND 
               cast(fg_anysense_max as float)/ bg_anysense_max > %(minfold)f
          ORDER BY fold_max DESC
          LIMIT %(nresults)s''' % locals()

        return statement

##########################################################################
##########################################################################
##########################################################################
# correlations
##########################################################################


class Correlations(ChipseqTracker):

    """Correlation between all sets.
    """

    pattern = "(.*)_annotations$"

    def __call__(self, track, slice=None):
        table = "%s_correlation" % self.field
        return self.getValues("SELECT %s AS %s FROM %s ORDER BY id" % (track,
                                                                       self.field,
                                                                       table))


class CorrelationsPeakval(Correlations):
    field = "peakval"


class CorrelationsAvgval(Correlations):
    field = "avgval"


class CorrelationsLength(Correlations):
    field = "length"

##########################################################################
##########################################################################
##########################################################################
# fold change
##########################################################################


class FoldChange(FoldChangeTracker):

    """return fold changes for all intervals.
    """
    pattern = "(.*)_readcounts$"

    def __call__(self, track):
        return self.getValues("SELECT CAST(fg_anysense_max AS float)/ bg_anysense_max AS fold FROM %(track)s_readcounts" % locals())

##########################################################################
##########################################################################
##########################################################################
# fold change
##########################################################################


class FoldChangeCounts(FoldChangeTracker):

    """Correlation between all sets.
    """
    pattern = "(.*)_readcounts$"
    mMinFoldChange = 2.0

    def __call__(self, track):
        data = []
        upfold = self.mMinFoldChange
        downfold = 1.0 / upfold
        data.append(("> %5.2f fold" % upfold, self.getValue(
            "SELECT COUNT(*) FROM %(track)s_readcounts WHERE CAST(fg_anysense_max AS float)/ bg_anysense_max > %(upfold)f " % locals())))
        data.append(("unchanged", self.getValue(
            "SELECT COUNT(*) FROM %(track)s_readcounts WHERE CAST(fg_anysense_max AS float)/ bg_anysense_max between %(downfold)f and %(upfold)f" % locals())))
        data.append(("< %5.2f fold" % downfold, self.getValue(
            "SELECT COUNT(*) FROM %(track)s_readcounts WHERE CAST(fg_anysense_max AS float)/ bg_anysense_max < %(downfold)f " % locals())))

        return odict(data)

##########################################################################
##########################################################################
##########################################################################
# peak shape
##########################################################################


class PeakShapeTracker(Tracker):

    '''return peakshape data.

    Only 1000 rows are returned.
    '''

    tracks = [os.path.basename(x)[:-len(".peakshape.tsv.gz")]
              for x in glob.glob(os.path.join(DATADIR, "*.peakshape.tsv.gz"))]
    slices = ["peak_height", "peak_width"]

    def __call__(self, track, slice=None):
        fn = os.path.join(
            DATADIR, "%(track)s.peakshape.tsv.gz.matrix_%(slice)s.gz" % locals())
        if not os.path.exists(fn):
            return

        matrix, rownames, colnames = iotools.readMatrix(iotools.openFile(fn))
        nrows = len(rownames)
        if nrows == 0:
            return

        if nrows > 1000:
            take = numpy.array(
                numpy.floor(numpy.arange(0, nrows, nrows / 1000)), dtype=int)
            rownames = [rownames[x] for x in take]
            matrix = matrix[take]

        return odict((('matrix', matrix),
                      ('rows', rownames),
                      ('columns', colnames)))
