# Path: nanopore/analyses/abstractAnalysis.py
# class AbstractAnalysis(Target):
#     """Base class to for analysis targets. Inherit this class to create an analysis.
#     """
#     def __init__(self, readFastqFile, readType, referenceFastaFile, samFile, outputDir):
#         Target.__init__(self)
#         self.readFastqFile = readFastqFile
#         self.referenceFastaFile = referenceFastaFile
#         self.samFile = samFile
#         self.outputDir = outputDir
#         self.readType = readType
#         print(str(self.samFile))
#         
#     def run(self):
#         """Base method that does some logging
#         """
#         logger.info("This analysis target has read fastq file: %s, reference fasta file: %s, sam file: %s and will output to the directory: %s" % \
#                     (self.readFastqFile, self.referenceFastaFile, self.samFile, self.outputDir))
#         
#     def finish(self):
#         """Method called when a analysis has finished successfully to indicate that it should not be repeated.
#         """
#         open(os.path.join(self.outputDir, "DONE"), 'w').close()
# 
#     @staticmethod
#     def reset(outputDir):
#         if AbstractAnalysis.isFinished(outputDir):
#             os.remove(os.path.join(outputDir, "DONE"))
#     
#     @staticmethod
#     def isFinished(outputDir):
#         return os.path.exists(os.path.join(outputDir, "DONE"))
#     
#     @staticmethod
#     def formatRatio(numerator, denominator):
#         if denominator == 0:
#             return float("nan")
#         return float(numerator)/denominator
# 
# Path: nanopore/analyses/utils.py
# def samIterator(sam):
#     """Creates an iterator over the aligned reads in a sam file, filtering out
#     any reads that have no reference alignment.
#     """
#     for aR in sam:
#         if aR.rname != -1:
#             yield aR


from nanopore.analyses.abstractAnalysis import AbstractAnalysis
from jobTree.src.bioio import fastqRead, fastaRead, system, reverseComplement
from nanopore.analyses.utils import samIterator
from collections import Counter
from math import log
import pysam, os, itertools

        for name, seq, qual in fastqRead(self.readFastqFile):
            for i in xrange(self.kmerSize, len(seq)):
                s = seq[ i - self.kmerSize : i ]
                if "N" not in s:
                    readKmers[s] += 1
                    readKmers[reverseComplement(s)] += 1

        return (refKmers, readKmers)

    def analyzeCounts(self, refKmers, readKmers, name):
        refSize, readSize = sum(refKmers.values()), sum(readKmers.values())
        outf = open(os.path.join(self.outputDir, name + "kmer_counts.txt"), "w")
        outf.write("kmer\trefCount\trefFraction\treadCount\treadFraction\tlogFoldChange\n")
        
        for kmer in itertools.product("ATGC", repeat=5):
            kmer = "".join(kmer)
            refFraction, readFraction = 1.0 * refKmers[kmer] / refSize, 1.0 * readKmers[kmer] / readSize
            if refFraction == 0:
                foldChange = "-Inf"
            elif readFraction == 0:
                foldChange = "Inf"
            else:
                foldChange = -log(readFraction / refFraction)
            outf.write("\t".join(map(str,[kmer, refKmers[kmer], refFraction, readKmers[kmer], readFraction, foldChange]))+"\n")
        outf.close()
        
        system("Rscript nanopore/analyses/kmer_analysis.R {} {} {} {} {}".format(os.path.join(self.outputDir, name + "kmer_counts.txt"), os.path.join(self.outputDir, name + "pval_kmer_counts.txt"), os.path.join(self.outputDir, name + "top_bot_sigkmer_counts.txt"), os.path.join(self.outputDir, name + "volcano_plot.pdf"), "Kmer"))

    def run(self, kmerSize=5):

# Path: nanopore/analyses/abstractAnalysis.py
# class AbstractAnalysis(Target):
#     """Base class to for analysis targets. Inherit this class to create an analysis.
#     """
#     def __init__(self, readFastqFile, readType, referenceFastaFile, samFile, outputDir):
#         Target.__init__(self)
#         self.readFastqFile = readFastqFile
#         self.referenceFastaFile = referenceFastaFile
#         self.samFile = samFile
#         self.outputDir = outputDir
#         self.readType = readType
#         print(str(self.samFile))
#         
#     def run(self):
#         """Base method that does some logging
#         """
#         logger.info("This analysis target has read fastq file: %s, reference fasta file: %s, sam file: %s and will output to the directory: %s" % \
#                     (self.readFastqFile, self.referenceFastaFile, self.samFile, self.outputDir))
#         
#     def finish(self):
#         """Method called when a analysis has finished successfully to indicate that it should not be repeated.
#         """
#         open(os.path.join(self.outputDir, "DONE"), 'w').close()
# 
#     @staticmethod
#     def reset(outputDir):
#         if AbstractAnalysis.isFinished(outputDir):
#             os.remove(os.path.join(outputDir, "DONE"))
#     
#     @staticmethod
#     def isFinished(outputDir):
#         return os.path.exists(os.path.join(outputDir, "DONE"))
#     
#     @staticmethod
#     def formatRatio(numerator, denominator):
#         if denominator == 0:
#             return float("nan")
#         return float(numerator)/denominator
# 
# Path: nanopore/analyses/utils.py
# def samIterator(sam):
#     """Creates an iterator over the aligned reads in a sam file, filtering out
#     any reads that have no reference alignment.
#     """
#     for aR in sam:
#         if aR.rname != -1:
#             yield aR


# Path: nanopore/analyses/kmerAnalysis.py
from nanopore.analyses.abstractAnalysis import AbstractAnalysis
from jobTree.src.bioio import fastqRead, fastaRead, system, reverseComplement
from nanopore.analyses.utils import samIterator
from collections import Counter
from math import log
import pysam, os, itertools


        for name, seq, qual in fastqRead(self.readFastqFile):
            for i in xrange(self.kmerSize, len(seq)):
                s = seq[ i - self.kmerSize : i ]
                if "N" not in s:
                    readKmers[s] += 1
                    readKmers[reverseComplement(s)] += 1

        return (refKmers, readKmers)

    def analyzeCounts(self, refKmers, readKmers, name):
        refSize, readSize = sum(refKmers.values()), sum(readKmers.values())
        outf = open(os.path.join(self.outputDir, name + "kmer_counts.txt"), "w")
        outf.write("kmer\trefCount\trefFraction\treadCount\treadFraction\tlogFoldChange\n")
        
        for kmer in itertools.product("ATGC", repeat=5):
            kmer = "".join(kmer)
            refFraction, readFraction = 1.0 * refKmers[kmer] / refSize, 1.0 * readKmers[kmer] / readSize
            if refFraction == 0:
                foldChange = "-Inf"
            elif readFraction == 0:
                foldChange = "Inf"
            else:
                foldChange = -log(readFraction / refFraction)
            outf.write("\t".join(map(str,[kmer, refKmers[kmer], refFraction, readKmers[kmer], readFraction, foldChange]))+"\n")
        outf.close()
        
        system("Rscript nanopore/analyses/kmer_analysis.R {} {} {} {} {}".format(os.path.join(self.outputDir, name + "kmer_counts.txt"), os.path.join(self.outputDir, name + "pval_kmer_counts.txt"), os.path.join(self.outputDir, name + "top_bot_sigkmer_counts.txt"), os.path.join(self.outputDir, name + "volcano_plot.pdf"), "Kmer"))

    def run(self, kmerSize=5):

        AbstractAnalysis.run(self)