# Path: checkm/util/seqUtils.py
# def readFasta(fastaFile, trimHeader=True):
#     '''Read sequences from FASTA file.'''
#     try:
#         if fastaFile.endswith('.gz'):
#             openFile = gzip.open
#         else:
#             openFile = open
# 
#         seqs = {}
#         for line in openFile(fastaFile):
#             # skip blank lines
#             if not line.strip():
#                 continue
# 
#             if line[0] == '>':
#                 if trimHeader:
#                     seqId = line[1:].split(None, 1)[0]
#                 else:
#                     seqId = line[1:].rstrip()
#                 seqs[seqId] = []
#             else:
#                 seqs[seqId].append(line[0:-1])
# 
#         for seqId, seq in seqs.items():
#             seqs[seqId] = ''.join(seq)
#     except:
#         logger = logging.getLogger('timestamp')
#         logger.error("Failed to process sequence file: {}".format(fastaFile))
#         sys.exit(1)
# 
#     return seqs
# 
# def baseCount(seq):
#     testSeq = seq.upper()
#     a = testSeq.count('A')
#     c = testSeq.count('C')
#     g = testSeq.count('G')
#     t = testSeq.count('T') + testSeq.count('U')
# 
#     return a, c, g, t


from checkm.plot.AbstractPlot import AbstractPlot
from checkm.util.seqUtils import readFasta, baseCount
from numpy import mean, array, log, poly1d, polyfit

        gc = []
        coverage = []
        for seqId, gcInfo in gcProfile.items():
            gc += gcInfo[1]
            coverage += coverageProfile[seqId][1]

        windowAxes.scatter(gc, coverage, c=abs(array(coverage)), s=10, lw=0.5, cmap='gray_r')
        windowAxes.set_xlabel('GC (mean = %.1f%%)' % (mean(gc) * 100))
        windowAxes.set_ylabel('Coverage (mean = %.1f)' % mean(coverage))

        # plot linear regression line
        if len(gc) > 1:
            slope, inter = polyfit(gc, coverage, 1)
            fit_fn = poly1d([slope, inter])  # fit_fn is now a function which takes in x and returns an estimate for y
            windowAxes.plot([min(gc), max(gc)], fit_fn([min(gc), max(gc)]), '--r', lw=0.5)
            windowAxes.set_title('GC vs. Coverage\n(window size = %d bp, slope = %.2f)' % (self.options.window_size, slope))
        else:
            # not possible to calculate best fit line
            windowAxes.set_title('GC vs. Coverage\n(window size = %d bp, no best fit line)' % self.options.window_size)

        # Prettify plot
        for a in windowAxes.yaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for a in windowAxes.xaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for line in windowAxes.yaxis.get_ticklines():

# Path: checkm/util/seqUtils.py
# def readFasta(fastaFile, trimHeader=True):
#     '''Read sequences from FASTA file.'''
#     try:
#         if fastaFile.endswith('.gz'):
#             openFile = gzip.open
#         else:
#             openFile = open
# 
#         seqs = {}
#         for line in openFile(fastaFile):
#             # skip blank lines
#             if not line.strip():
#                 continue
# 
#             if line[0] == '>':
#                 if trimHeader:
#                     seqId = line[1:].split(None, 1)[0]
#                 else:
#                     seqId = line[1:].rstrip()
#                 seqs[seqId] = []
#             else:
#                 seqs[seqId].append(line[0:-1])
# 
#         for seqId, seq in seqs.items():
#             seqs[seqId] = ''.join(seq)
#     except:
#         logger = logging.getLogger('timestamp')
#         logger.error("Failed to process sequence file: {}".format(fastaFile))
#         sys.exit(1)
# 
#     return seqs
# 
# def baseCount(seq):
#     testSeq = seq.upper()
#     a = testSeq.count('A')
#     c = testSeq.count('C')
#     g = testSeq.count('G')
#     t = testSeq.count('T') + testSeq.count('U')
# 
#     return a, c, g, t


# Path: checkm/plot/gcBiasPlots.py
from checkm.plot.AbstractPlot import AbstractPlot
from checkm.util.seqUtils import readFasta, baseCount
from numpy import mean, array, log, poly1d, polyfit

        gc = []
        coverage = []
        for seqId, gcInfo in gcProfile.items():
            gc += gcInfo[1]
            coverage += coverageProfile[seqId][1]

        windowAxes.scatter(gc, coverage, c=abs(array(coverage)), s=10, lw=0.5, cmap='gray_r')
        windowAxes.set_xlabel('GC (mean = %.1f%%)' % (mean(gc) * 100))
        windowAxes.set_ylabel('Coverage (mean = %.1f)' % mean(coverage))

        # plot linear regression line
        if len(gc) > 1:
            slope, inter = polyfit(gc, coverage, 1)
            fit_fn = poly1d([slope, inter])  # fit_fn is now a function which takes in x and returns an estimate for y
            windowAxes.plot([min(gc), max(gc)], fit_fn([min(gc), max(gc)]), '--r', lw=0.5)
            windowAxes.set_title('GC vs. Coverage\n(window size = %d bp, slope = %.2f)' % (self.options.window_size, slope))
        else:
            # not possible to calculate best fit line
            windowAxes.set_title('GC vs. Coverage\n(window size = %d bp, no best fit line)' % self.options.window_size)

        # Prettify plot
        for a in windowAxes.yaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for a in windowAxes.xaxis.majorTicks:
            a.tick1On = True
            a.tick2On = False

        for line in windowAxes.yaxis.get_ticklines():

            line.set_color(self.axesColour)
