# Path: seqio.py
# def fastaIterator(fh):
# 
#     l = fh.readline()
#     if(not l or not l.startswith(">")):
#         raise Exception("No \">\" at start of Fasta File")
#     name = l.strip()[1:]
#     seq = ""
#     while True:
#         l = fh.readline()
#         if not l or l.startswith(">"):
#             yield FastaRecord(name,seq)
#             if not l:
#                 break
#             name = l.strip()[1:]
#             seq = ""
#         else:
#             seq += l.strip()
# 
# Path: nucio.py
# def lineRecordIterator(fh, nt, nt_types, filter_func=trueFunc, delim=None, cleaner_func=None):
# def lineItemIterator(fh, filter_func=trueFunc):
# def getNucmerAlignmentIterator(fh):
# def fileIterator(filename, itemIterator, open_func=open):
# def recordToString( record, delim="\t" ):
# def nucRecordToString(nuc_record):
# def typeify(mylist, mytypes):
# def prepareDeltaRecord(buf):
# def getDeltaAlignmentIterator(fh):
# def deltaRecordHeaderToString(dr):
# def deltaAlignmentHeaderToString(alignment):
# def deltaRecordToOriginalFormat(dr):
#     def __init__(self, s, *args ):
#     def __enter__(self):
#     def __exit__(self, ttype ,value, traceback):
# def openerFromExtension(filename, default=None):
# class FileOrStream:
# 
# Path: cov.py
# def fillc(l,s,e):
#     '''l is an array
#     s and e are start and end of
#     a region the coverage should be incremented in'''
# 
#     if(len(l) < e):
#         l += list(repeat(0, e-len(l)))
#     for cr in range(s,e):
#         l[cr] += 1
# 
# def getMarkedRanges(v):
#     '''Takes an array v with wanted elements marked as 1
#     all other elements are 0. This function returns
#     the index ranges of these 1's'''
#     
#     #subtract adjacent elements to find -1
#     #flip these to 1's and change everything else to 0
#     #so that we can use it to find the index
#     breaks = map(lambda x: 1 if x == -1 else 0 , pairwise(v + [0],operator.sub))
# 
#     #create a cumsum of the inverse to know how many elements
#     #are a part of this region
#     lengths = accumulate_mod(v)
# 
#     #zip up the indexes with the lengths
#     z = zip(count(),list(lengths))
# 
#     #use the breaks to select only the indexes we want
#     #and use the cumsum to know how many elements came
#     #previously
#     endAndLength = compress(z,breaks)
# 
#     return map(lambda (e,l): (e-l+1,e), endAndLength)
# 
# def getCoverageFromNucAlignments(alignments):
#     '''Gets coverage from nucmer alignments'''
#     it = iter(alignments)
#     n = it.next()
# 
#     def nucAGetter():
#         yield (n.sstart-1, n.send-1)
#         for g in alignments:
#             yield (n.sstart-1, n.send-1)
# 
#     return getCoverage( nucAGetter(), n.slen)
# 
# Path: gccontent.py
# def getGCSlidingWindow(sequence, window_size):
#     '''Calculates %GC in sequence of sliding window window_size'''
#     return map(lambda x: (x.count("G")+x.count("C")) / float(window_size), 
#                map(lambda r: sequence[r:r+window_size], 
#                    range(len(sequence)-window_size+1)))


import sys
import operator
from seqio import fastaIterator
from itertools import groupby,izip,imap
from nucio import NucRecord, NucRecordTypes, getNucmerAlignmentIterator
from cov import fillc, getMarkedRanges, getCoverageFromNucAlignments
from gccontent import getGCSlidingWindow
#!/usr/bin/env python




GC_WINDOW_SIZE = 300
GC_THRESHOLD = 0.7
MIN_COV_GAP = 100

if not len(sys.argv) == 4:
    print "gc_count.py reads.fa alignments.sc outprefix"
    sys.exit(1)


rfh = open(sys.argv[1])
afh = open(sys.argv[2])
ofh = open(sys.argv[3]+".uncov.gc.bases","w")


# Path: seqio.py
# def fastaIterator(fh):
# 
#     l = fh.readline()
#     if(not l or not l.startswith(">")):
#         raise Exception("No \">\" at start of Fasta File")
#     name = l.strip()[1:]
#     seq = ""
#     while True:
#         l = fh.readline()
#         if not l or l.startswith(">"):
#             yield FastaRecord(name,seq)
#             if not l:
#                 break
#             name = l.strip()[1:]
#             seq = ""
#         else:
#             seq += l.strip()
# 
# Path: nucio.py
# def lineRecordIterator(fh, nt, nt_types, filter_func=trueFunc, delim=None, cleaner_func=None):
# def lineItemIterator(fh, filter_func=trueFunc):
# def getNucmerAlignmentIterator(fh):
# def fileIterator(filename, itemIterator, open_func=open):
# def recordToString( record, delim="\t" ):
# def nucRecordToString(nuc_record):
# def typeify(mylist, mytypes):
# def prepareDeltaRecord(buf):
# def getDeltaAlignmentIterator(fh):
# def deltaRecordHeaderToString(dr):
# def deltaAlignmentHeaderToString(alignment):
# def deltaRecordToOriginalFormat(dr):
#     def __init__(self, s, *args ):
#     def __enter__(self):
#     def __exit__(self, ttype ,value, traceback):
# def openerFromExtension(filename, default=None):
# class FileOrStream:
# 
# Path: cov.py
# def fillc(l,s,e):
#     '''l is an array
#     s and e are start and end of
#     a region the coverage should be incremented in'''
# 
#     if(len(l) < e):
#         l += list(repeat(0, e-len(l)))
#     for cr in range(s,e):
#         l[cr] += 1
# 
# def getMarkedRanges(v):
#     '''Takes an array v with wanted elements marked as 1
#     all other elements are 0. This function returns
#     the index ranges of these 1's'''
#     
#     #subtract adjacent elements to find -1
#     #flip these to 1's and change everything else to 0
#     #so that we can use it to find the index
#     breaks = map(lambda x: 1 if x == -1 else 0 , pairwise(v + [0],operator.sub))
# 
#     #create a cumsum of the inverse to know how many elements
#     #are a part of this region
#     lengths = accumulate_mod(v)
# 
#     #zip up the indexes with the lengths
#     z = zip(count(),list(lengths))
# 
#     #use the breaks to select only the indexes we want
#     #and use the cumsum to know how many elements came
#     #previously
#     endAndLength = compress(z,breaks)
# 
#     return map(lambda (e,l): (e-l+1,e), endAndLength)
# 
# def getCoverageFromNucAlignments(alignments):
#     '''Gets coverage from nucmer alignments'''
#     it = iter(alignments)
#     n = it.next()
# 
#     def nucAGetter():
#         yield (n.sstart-1, n.send-1)
#         for g in alignments:
#             yield (n.sstart-1, n.send-1)
# 
#     return getCoverage( nucAGetter(), n.slen)
# 
# Path: gccontent.py
# def getGCSlidingWindow(sequence, window_size):
#     '''Calculates %GC in sequence of sliding window window_size'''
#     return map(lambda x: (x.count("G")+x.count("C")) / float(window_size), 
#                map(lambda r: sequence[r:r+window_size], 
#                    range(len(sequence)-window_size+1)))


# Path: gc_count.py
import sys
import operator
from seqio import fastaIterator
from itertools import groupby,izip,imap
from nucio import NucRecord, NucRecordTypes, getNucmerAlignmentIterator
from cov import fillc, getMarkedRanges, getCoverageFromNucAlignments
from gccontent import getGCSlidingWindow

#!/usr/bin/env python




GC_WINDOW_SIZE = 300
GC_THRESHOLD = 0.7
MIN_COV_GAP = 100

if not len(sys.argv) == 4:
    print "gc_count.py reads.fa alignments.sc outprefix"
    sys.exit(1)


rfh = open(sys.argv[1])
afh = open(sys.argv[2])
ofh = open(sys.argv[3]+".uncov.gc.bases","w")


reads = {}