annotate microsatellite.py @ 0:70f8259b0b30 draft

Uploaded

#!/usr/bin/env python
"""
Snoop thru a fasta file looking for microsatellite repeats of given periods
Output format: length_of_repeat left_flank_length right_flank_length  repeat_motif  hamming_distance  read_name read_sequence read_quality  (additional columns)

If --r option turned on, output format will have additional columns behind:
read_name read_chr  pre_s pre_e tr_s  tr_e  suf_s suf_e tr_len  tr_ref_seq

pre_s           where the read start
pre_e           the last position before microsatellite
tr_s            where microsatellite start
tr_e            where microsatellite end
suf_s           first base after microsatellite
tr_ref_seq      reference sequence corresponding to microsatellite

* output positions are 0 based

:Author: Chen Sun ([e-mail]); Bob Harris ([e-mail])
          
modifing log:

09/27/2013
replace function dense_intervals with function non_negative_intervals, which do not need to import such file.

10/18/2013
modify function find_repeat_element to get a quick speed, under the condition that hamming_distance = 0, which means do not allowed any mutation/indel

02/25/2014
add function that can deal with mapped reads
with additional output

02/28/2014
modify the 0-based end point, as in 0-base area, it is half-open [ )
so the 0-based site, should always be added by 1

03/05/2014
deal with multi-fasta
"""
from sys          import argv,stdin,stderr,exit
from string       import maketrans
from md5          import new as md5_new
import re
#from pyfracluster import dense_intervals

def usage(s=None):
    message = """
usage: microsat_snoop [fasta_file] [options]
  <fasta_file>                Name of file to read sequences from;  if absent,
                              sequences are read from stdin
  --fasta                     Input file is in fasta format
                              (this is the default)
  --fastq                     Input file is in fastq format
                              (default is fasta unless filename is .fastq)
  --fastq:noquals             Input file is in fastq format, but discard quals
  --sam                       Input file is SAM file 
  --r                         Indicate additional output information, if indicated,
                              --ref option is mendatory
  --ref=<filepath>            Reference file (absolute) path
  --period=<length>           (mandatory,cumulative) repeat length(s) to be
                              searched for
                              <length> is expected to be small, less than 10
                              <length> can also be a comma-separated list, or
                              a range <low>..<high>
  --rate=<fraction>           control the candidate repeat interval detector;
                              it will consider intervals with at least
                              <fraction> of matches when shifted by the period;
                              <fraction> is between 0 and 1 and can be either a
                              real number or <n>/<d>
                              (default is 6/7)
  --minlength=<length>        minimum length of intervals reported, in bp
                              (default is 20)
  --progress=<count>          how often to report the sequence we're searching
                              (default is no progress report)
  --allowduplicates           process all input sequences
                              (this is the default)
  --noduplicates              ignore any input sequence that's the same as an
                              earlier sequence
  --nonearduplicates          ignore any input sequence that has the same first
                              100 bp as an earlier sequence
  --nonearduplicate=<length>  ignore any input sequence that has the same first
                              <length> bp as an earlier sequence
  --hamming=<count>           Don't report candidate repeat intervals that have
                              more than <count> mismatches
                              (default is to do no such filtering)
  --prefix=<length>           Don't report candidate repeat intervals that
                              start within <length> of the sequence start
                              (default is to do no such filtering)
  --suffix=<length>           Don't report candidate repeat intervals that
                              end within <length> of the sequence end
                              (default is to do no such filtering)
  --subsample=<k>/<n>         Process only the <k>th sequence of every group of
                              <n> sequences;  <k> ranges from 1 to <n>
  --multipleruns              Consider all candidate intervals in a sequence
                              (default is to consider only the longest)
  --partialmotifs             Consider microatelites with a partial motif
                              (default is to consider only whole motifs)
  --splitbyvalidity           Preprocess sequences, splitting at Ns;  this
                              prevents candidates from including Ns
                              (default is not to split)
  --noflankdisplay            Show entire sequence as flanking regions
                              (this is the default)
  --flankdisplay=<length>     Limit length of flanking regions shown
  --readnamesuffix=<string>   Root of suffix to append to read names;  e.g. 1
                              for forward, 2 for reverse;  this triggers other
                              info to be included in the suffix
                              (default is "1" for fastq;  no suffix for fasta)
  --head=<number>             limit the number of sequences processed
  --markend                   Write a marker line upon completion
                              (default is not to write a marker)
  --help=details              Describe the process, and quit"""

    if (s == None): exit (message)
    else:           exit ("%s\n%s" % (s,message))


detailedDescription = """In broad terms, the process works as follows:

(1) Identify intervals that are highly correlated with the interval shifted by
    P (the repeat period).  These intervals are called "runs" or "candidates".
    The level of correlation required is controlled by rateThreshold. 
    Depending on whether we want to look for more than one microsat, we either
    find the longest such run (simple algorithm) or many runs (more complicated
    algorithm). The following steps are then performed on each run.

(2) Find the most likely repeat motif in the run.  This is done by counting
    all kmers (of length P) and choosing the most frequent.  If that kmer is
    itself covered by a sub-repeat we discard this run.  The idea is that we
    can ignore a 6-mer like ACGACG because we will find it when we are looking
    for 3-mers.

(3) Once we identify the most likely repeat motif, we then modify the
    interval, adjusting start and end to find the interval that has the fewest
    mismatches vs. a sequence of the motif repeated (hamming distance).  Only
    whole copies of the motif are considered.

(4) At this point we have a valid microsat interval (in the eyes of the
    program). It is subjected to some filtering stages (hamming distance or too
    close to an end), and if it satisfies those conditions, it's reported to
    the user."""

def main():
    global debug

    #=== parse the command line ===

    inputFilename         = None
    referenceFileName     = None #add by Chen Sun on 02/25
    inputFormat           = None
    repeatPeriods         = []
    rateThreshold         = 6 / 7.0
    lengthThreshold       = 20
    reportProgress        = None
    discardDuplicates     = False
    discardNearDuplicates = False
    nearDuplicatePrefix   = 100
    hammingThreshold      = 0
    prefixThreshold       = None
    suffixThreshold       = None
    subsampleK            = None
    subsampleN            = None
    reportMultipleRuns    = False
    allowPartialMotifs    = False
    splitByValidity       = False
    flankDisplayLimit     = None
    readNameSuffix        = None
    headLimit             = None
    markEndOfFile         = False
    additionalInfo        = False
    debug                 = []

    for arg in argv[1:]:
        if (arg == "--fasta"):
            inputFormat = "fasta"
        elif (arg == "--fastq"):
            inputFormat = "fastq"
        elif (arg == "--fastq:noquals"):
            inputFormat = "fastq:noquals"
        elif (arg == "--sam"):
            inputFormat = "sam"
        elif (arg == "--r"):
            additionalInfo = True
        elif (arg.startswith("--ref=")):
            referenceFileName = arg.split("=",1)[1]
        elif (arg.startswith("--period=")):
            val = arg.split("=",1)[1]
            for period in val.split(","):
                if (".." in period):
                    (lowPeriod,highPeriod) = period.split("..",1)
                    lowPeriod  = int(lowPeriod)
                    highPeriod = int(highPeriod)
                    for period in xrange(lowPeriod,highPeriod+1):
                        repeatPeriods += [period]
                else:
                    repeatPeriods += [int(period)]
        elif (arg.startswith("--rate=")):
            val = arg.split("=",1)[1]
            rateThreshold = float_or_fraction(val)
            assert (0.0 < rateThreshold <= 1.0), "%s not a valid rate" % val
        elif (arg.startswith("--minlength=")):
            val = arg.split("=",1)[1]
            lengthThreshold = int(val)
            assert (lengthThreshold >= 0)
        elif (arg.startswith("--progress=")):
            val = arg.split("=",1)[1]
            reportProgress = int(val)
        elif (arg == "--allowduplicates"):
            discardDuplicates     = False
            discardNearDuplicates = False
        elif (arg == "--noduplicates"):
            discardDuplicates     = True
            discardNearDuplicates = False
        elif (arg == "--nonearduplicates"):
            discardDuplicates     = False
            discardNearDuplicates = True
        elif (arg.startswith("--nonearduplicate=")):
            val = arg.split("=",1)[1]
            discardDuplicates     = False
            discardNearDuplicates = True
            nearDuplicatePrefix   = int(val)
            assert (nearDuplicatePrefix > 0)
        elif (arg.startswith("--hamming=")):
            val = arg.split("=",1)[1]
            hammingThreshold = int(val)
            assert (hammingThreshold >= 0)
        elif (arg.startswith("--prefix=")):
            val = arg.split("=",1)[1]
            prefixThreshold = int(val)
            assert (prefixThreshold >= 0)
        elif (arg.startswith("--suffix=")):
            val = arg.split("=",1)[1]
            suffixThreshold = int(val)
            assert (suffixThreshold >= 0)
        elif (arg.startswith("--subsample=")):
            val = arg.split("=",1)[1]
            (k,n) = val.split("/",2)
            subsampleK = int(k)
            subsampleN = int(n)
            assert (0 < subsampleK <= subsampleN)
        elif (arg == "--multipleruns"):
            reportMultipleRuns = True
        elif (arg == "--partialmotifs"):
            allowPartialMotifs = True
        elif (arg == "--splitbyvalidity"):
            splitByValidity = True
        elif (arg == "--noflankdisplay"):
            flankDisplayLimit = None
        elif (arg.startswith("--flankdisplay=")):
            val = arg.split("=",1)[1]
            flankDisplayLimit = int(val)
            assert (flankDisplayLimit >= 0)
        elif (arg.startswith("--readnamesuffix")):
            readNameSuffix = arg.split("=",1)[1]
        elif (arg.startswith("--head=")):
            headLimit = int_with_unit(arg.split("=",1)[1])
        elif (arg == "--markend"):
            markEndOfFile = True
        elif (arg == "--help=details"):
            exit (detailedDescription)
        elif (arg.startswith("--debug=")):
            debug += (arg.split("=",1)[1]).split(",")
        elif (arg.startswith("--")):
            usage("unrecognized option: %s" % arg)
        elif (inputFilename == None):
            inputFilename = arg
        else:
            usage("unrecognized option: %s" % arg)

    #=== determine periods of interest ===

    if (repeatPeriods == []):
        usage("you gotta give me a repeat period")

    if (additionalInfo == True):
        if (referenceFileName == None):
            usage("reference file path needed. use --ref=<reference> to indicate")

    periodSeed = {}
    for period in repeatPeriods:
        if (period < 1): usage("period %d is not valid" % period)
        periodSeed[period] = True

    repeatPeriods = [period for period in periodSeed]
    repeatPeriods.sort()

    #=== determine input format ===

    if   (inputFormat == "fasta"):           sequence_reader = fasta_sequences
    elif (inputFormat == "fastq"):           sequence_reader = fastq_sequences
    elif (inputFormat == "fastq:noquals"):   sequence_reader = fastq_sequences
    elif (inputFormat == "sam"):             sequence_reader = sam_sequences
    elif (inputFilename == None):            sequence_reader = fasta_sequences
    elif (inputFilename.endswith(".fastq")): sequence_reader = fastq_sequences
    elif (inputFilename.endswith(".fq")):    sequence_reader = fastq_sequences
    elif (inputFilename.endswith(".sam")):   sequence_reader = sam_sequences
    else:                                    sequence_reader = fasta_sequences

    if (inputFilename != None): inputF = file(inputFilename,"rt")
    else:                       inputF = stdin

    if   (readNameSuffix == None) \
     and (sequence_reader == fastq_sequences) \
     and (inputFormat != "fastq:noquals"):
        readNameSuffix = "1"

    #=== process the sequences ===
    
    refSequence = {}
    rightName = ""
    sequence = ""
    if additionalInfo:
        firstFasta = True
        originalRefF = open(referenceFileName)
        for line in originalRefF.readlines():
            line = line.replace('\r','')
            line = line.replace('\n','')
            if line.startswith(">"):
                if firstFasta:
                    firstFasta = False
                else:
                    refSequence[rightName] = sequence
                rightName = line[1:]
                sequence = ""
                continue
            sequence += line
        originalRefF.close()
        refSequence[rightName] = sequence

    sequenceSeen = {}

    numSequences = 0
    for seqInfo in sequence_reader(inputF):
        numSequences += 1
        if (headLimit != None) and (numSequences > headLimit):
            print >>stderr, "limit of %d sequences reached" % headLimit
            break

        if (sequence_reader == sam_sequences):
            #seqName,"".join(seqNucs).upper().translate(nonDnaMap), refName, pre_s, cigar
            (name, sequence, refName, pre_s, cigar) = seqInfo
            quals = None
        elif (sequence_reader == fastq_sequences):
            (name,sequence,quals) = seqInfo
            if (inputFormat == "fastq:noquals"): quals = None
        else:
            (name,sequence) = seqInfo
            quals = None

        if (reportProgress != None) and (numSequences % reportProgress == 0):
            print >>stderr, "%s %d" % (name,numSequences)

        # if we're subsampling and not interested in this sequence, skip it

        if (subsampleN != None):
            if ((numSequences-1) % subsampleN != (subsampleK-1)):
                continue

        # if this sequence is shorter than the length of interest, skip it

        seqLen = len(sequence)
        if (seqLen < period) or (seqLen < lengthThreshold): continue

        # if we're not interested in duplicates and this is one, skip it;
        # note that we assume no hash collisions occur, i.e. that all hash
        # matches are truly sequence matches

        if (discardDuplicates):
            h = hash108(sequence)
            if (h in sequenceSeen): continue
            sequenceSeen[h] = True
        elif (discardNearDuplicates):
            h = hash108(sequence[:nearDuplicatePrefix])
            if (h in sequenceSeen): continue
            sequenceSeen[h] = True

        # split the sequence into chunks of valid nucleotides

        if (splitByValidity):
            chunks = [(start,end) for (start,end) in nucleotide_runs(sequence)]
        else:
            chunks = [(0,len(sequence))]

        # evaluate for each period of interest

        for period in repeatPeriods:

            # operate on each chunk

            for (chunkStart,chunkEnd) in chunks:
                chunkLen = chunkEnd - chunkStart
                if (chunkLen < period) or (chunkLen < lengthThreshold): continue

                if ("validity" in debug) or ("correlation" in debug) or ("runs" in debug):
                    print >>stderr, ">%s_%d_%d" % (name,chunkStart,chunkEnd)

                # compute correlation sequence

                corr = correlation_sequence(sequence,period,chunkStart,chunkEnd)

                if ("correlation" in debug) or ("runs" in debug):
                    print >>stderr, sequence[chunkStart:chunkEnd]
                    print >>stderr, corr

                # find runs (candidates for being a microsat) 

                if (reportMultipleRuns):
                    runs = all_suitable_runs(corr,lengthThreshold-period,rateThreshold, hammingThreshold)
                else:
                    runs = longest_suitable_run(corr,lengthThreshold,rateThreshold)
                if (runs == []): continue


                if ("runs" in debug):
                    for (start,end) in runs:
                        run = [" "] * seqLen
                        for ix in xrange(start-period,end):
                            run[ix] = "*"
                        print >>stderr, "".join(run)

                if ("candidates" in debug):
                    for (start,end) in runs:
                        print >>stderr, "%s %d %d" % (name,start,end)

                # process runs and report those that pass muster

                runCount = 0
                for (start,end) in runs:
                    runCount += 1

                    start = chunkStart + start - period
                    end   = chunkStart + end

                    (kmer,d,start,end) = find_repeat_element(hammingThreshold, period,sequence,start,end,allowPartials=allowPartialMotifs)
                    if (kmer == None): continue    # (no useful repeat kmer was found)

                    rptExtent = end - start
                    prefixLen = start
                    suffixLen = seqLen - end
                    if (rptExtent <= period): continue
                    if (hammingThreshold != None) and (d         > hammingThreshold): continue
                    if (prefixThreshold  != None) and (prefixLen < prefixThreshold):  continue
                    if (suffixThreshold  != None) and (suffixLen < suffixThreshold):  continue

                    if (flankDisplayLimit == None):
                        seq = sequence[:start] \
                            + sequence[start:end].lower() \
                            + sequence[end:]
                    else:
                        seq = sequence[max(chunkStart,start-flankDisplayLimit):start] \
                            + sequence[start:end].lower() \
                            + sequence[end:min(chunkEnd,end+flankDisplayLimit)]
                    reportName = name
                    if (readNameSuffix != None):
                        reportName += "_"+readNameSuffix+"_per"+str(period)+"_"+str(runCount)
                    if (quals == None or quals == "." or quals == "\t."): quals = "\t."
                    else:               quals = "\t" + quals
                    if not additionalInfo:
                        print "%d\t%d\t%d\t%s\t%d\t%s\t%s%s" \
                            % (rptExtent,prefixLen,suffixLen,kmer,d,reportName,seq,quals)
                    else:
                        #pre_e = pre_s + prefixLen - 1
                        refPoint = pre_s
                        donorPoint = 0
                        
                        donorBeforeStart = prefixLen - 1 #pre_e
                        donorMicroStart = prefixLen     #tr_s
                        donorMicroEnd = donorMicroStart + rptExtent - 1 #tr_e
                        donorAfterMicro = donorMicroEnd + 1 #suf_s
                        donorEnd = len(seq) - 1    #suf_e
                        
                        set_pre_e = False
                        set_tr_s = False
                        set_tr_e = False
                        set_suf_s = False
                        set_suf_e = False
                        
                        pre_e = 0
                        tr_s = 0
                        tr_e = 0
                        suf_s = 0
                        suf_e = 0
                        
                        matchList = re.findall('(\d+)([IDM])', cigar)
                        unCognitiveCigar = False
                        for matchN, matchType in matchList:
                            matchNum = int(matchN)
                            if matchType == "M":
                                donorPoint = donorPoint + matchNum
                                refPoint = refPoint + matchNum
                            elif matchType == "D":
                                refPoint = refPoint + matchNum
                                continue
                            elif matchType == "I":
                                donorPoint = donorPoint + matchNum
                            else:
                                unCognitiveCigar = True
                                break

                            if not set_pre_e:
                                if donorPoint >= donorBeforeStart:
                                    pre_e = refPoint - (donorPoint - donorBeforeStart)
                                    set_pre_e = True
                                else:
                                    continue
                                    
                            if not set_tr_s:
                                if donorPoint >= donorMicroStart:
                                    tr_s = refPoint - (donorPoint - donorMicroStart)
                                    set_tr_s = True
                                else:
                                    continue
                                    
                            if not set_tr_e:
                                if donorPoint >= donorMicroEnd:
                                    tr_e = refPoint - (donorPoint - donorMicroEnd)
                                    set_tr_e = True
                                else:
                                    continue
                                    
                            if not set_suf_s:
                                if donorPoint >= donorAfterMicro:
                                    suf_s = refPoint - (donorPoint - donorAfterMicro)
                                    set_suf_s = True
                                else:
                                    continue
                                    
                            if not set_suf_e:
                                if donorPoint >= donorEnd:
                                    suf_e = refPoint - (donorPoint - donorEnd)
                                    set_suf_e = True
                                else:
                                    continue
                                
                        if unCognitiveCigar:
                            break
                        tr_len = tr_e - tr_s + 1

                        if refName not in refSequence:
                            tr_ref_seq = "."
                        else:
                            if refSequence[refName] == "":
                                tr_ref_seq = "."
                            elif len(refSequence[refName]) <= tr_e:
                                tr_ref_seq = "."
                            else:
                                tr_ref_seq = refSequence[refName][tr_s:tr_e+1]

                        pre_e += 1
                        tr_e += 1
                        suf_e += 1
                        print "%d\t%d\t%d\t%s\t%d\t%s\t%s%s\t%s\t%s\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%s" \
                            % (rptExtent,prefixLen,suffixLen,kmer,d,reportName,seq,quals,reportName,refName,pre_s,pre_e,tr_s,tr_e,suf_s,suf_e,tr_len,tr_ref_seq)

    if (markEndOfFile):
        print "# microsat_snoop end-of-file"

    if (inputF != stdin):
        inputF.close()

# non_negative_intervals
#     find intervals with exactly + and no -
#     from string like this : +++++++++---+++++++++
def non_negative_intervals(seq, minLength=None):

    start = -1
    end = -1
    firstPlus = 1
    #print seq
    for ix in range(len(seq)): # for every char in seq    
        ch = seq[ix]
        if(ch == "+"):
            if(firstPlus):
                firstPlus = 0
                start = ix
            else:
                continue
        elif(ch == "-"):
            if(start >= 0):
                end = ix-1
                if((end - start + 1) >= minLength):
                    yield (start,end+1)
                start = -1
                firstPlus = 1
    if(start > 0):
        if((ix - start + 1) >= minLength):
            yield (start, ix+1)


###################################################################
# modified by Chen Sun on 7/11/2014
# We do not want other modules, so parse these functions inside
#
###################################################################

# parse a string of the form {positives}/{positives_and_neutrals}

def parse_spec(s):
    if ("/" not in s): raise ValueError
    (n,d) = s.split("/",1)
    if (not n.startswith("{")) or (not n.endswith("}")): raise ValueError
    if (not d.startswith("{")) or (not d.endswith("}")): raise ValueError

    positives = n[1:-1]
    d         = d[1:-1]

    for ch in positives:
        if (ch not in d): raise ValueError

    neutrals = [ch for ch in d if (ch not in positives)]
    return (positives,neutrals)


# convert a string to a number, allowing fractions

def float_or_fraction(s):
    if ("/" in s):
        (numer,denom) = s.split("/",1)
        return float(numer)/float(denom)
    else:
        return float(s)


# dense_intervals--
#    Find all non-overlapping runs with a good enough rate (of positives), and
#    which meet our length threshold.
#
#    The algorithm used is adapted from Zhang, Berman, Miller, "Post-processing
#    long pairwise alignments", Bioinformatics Vol. 15 no. 12 1999.
#
# $$$ we use the denominator as the threshold, but we really should use the
# $$$ .. numerator, comparing it to minLength*rate

def dense_intervals(seq,rate,positives,neutrals,blockers="",minLength=None):

    if (blockers == None):
        blockers = "".join([chr(n) for n in range(1,256)
                                   if  (chr(n) not in positives)
                                   and (chr(n) not in neutrals)])

    stackLeft       = [None]    # stack with each entry containing five
    stackRight      = [None]    # .. elements;  note that entry zero is not
    stackLeftScore  = [None]    # .. used
    stackRightScore = [None]
    stackLower      = [None]
    top   = 0
    score = 0

    for ix in range(len(seq)):
        ch = seq[ix]
        if (ch in blockers):
            # emit intervals

            for sp in range(1,top+1):
                left  = stackLeft [sp] + 1
                right = stackRight[sp]

                while (left < right) and (seq[left]  not in positives): left  += 1
                while (right > left) and (seq[right] not in positives): right -= 1

                right += 1
                if (minLength == None) or (right - left >= minLength):
                    yield (left,right)

            #empty stack

            stackLeft       = [None]
            stackRight      = [None]
            stackLeftScore  = [None]
            stackRightScore = [None]
            stackLower      = [None]
            top   = 0
            score = 0
            continue

        if   (ch in positives): weight = 1-rate
        elif (ch in neutrals):  weight = -rate
        else: raise ValueError

        score += weight
        #if ("algorithm" in debug):
        #    print >>sys.stderr, "%3d: %c %5.2f" % (ix, ch, score),

        if (weight < 0):
            #if ("algorithm" in debug):
            #    print >>sys.stderr
            continue

        if (top > 0) and (stackRight[top] == ix-1):
            # add this site to the interval on top of the stack

            stackRight     [top] = ix
            stackRightScore[top] = score

            #if ("algorithm" in debug):
            #    print >>sys.stderr, \
            #          " extending [%d] %d-%d %4.1f %4.1f" \
            #        % (top,
            #           stackLeft     [top], stackRight     [top],
            #           stackLeftScore[top], stackRightScore[top]),

        else:
            # create a one site interval

            top += 1
            if (top >= len(stackLeft)):
                stackLeft       += [None]
                stackRight      += [None]
                stackLeftScore  += [None]
                stackRightScore += [None]
                stackLower      += [None]

            stackLeft      [top] = ix - 1
            stackLeftScore [top] = score - weight
            stackRight     [top] = ix
            stackRightScore[top] = score
            stackLower     [top] = top - 1

            while (stackLower[top] > 0) \
              and (stackLeftScore[stackLower[top]] > stackLeftScore[top]):
                stackLower[top] = stackLower[stackLower[top]]

            #if ("algorithm" in debug):
            #    print >>sys.stderr, \
            #          " creating  [%d] %d-%d %4.1f %4.1f -> %d" \
            #        % (top,
            #           stackLeft     [top], stackRight     [top],
            #           stackLeftScore[top], stackRightScore[top],
            #           stackLower    [top]),

        # merge intervals;  if there is a previous interval with a no-higher
        # left score and no-higher right score, merge this interval (and all
        # intervening ones) into that one

        while (top > 1) \
          and (stackLower[top] > 0) \
          and (stackRightScore[stackLower[top]] <= stackRightScore[top]):
            stackRight     [stackLower[top]] = stackRight     [top]
            stackRightScore[stackLower[top]] = stackRightScore[top]
            top = stackLower[top]

            #if ("algorithm" in debug):
            #    print >>sys.stderr, \
            #          "\n%*s merging   [%d] %d-%d %4.1f %4.1f" \
            #        % (13, "", top,
            #           stackLeft[top],      stackRight     [top],
            #           stackLeftScore[top], stackRightScore[top]),

        #if ("algorithm" in debug):
        #    print >>sys.stderr

    # emit intervals

    for sp in range(1,top+1):
        left  = stackLeft [sp] + 1
        right = stackRight[sp]

        while (left < right) and (seq[left]  not in positives): left  += 1
        while (right > left) and (seq[right] not in positives): right -= 1

        right += 1
        if (minLength == None) or (right - left >= minLength):
            yield (left,right)
            
            
###################################################################
# modified by Chen Sun on 7/11/2014
#
###################################################################

# correlation_sequence--
#    Compute the correlation sequence for a given period.  This is a sequence
#    of + and - indicating whether the base at a given position matches the one
#    P positions earlier (where P is the period).  The first P positions are
#    blank.  Positions with single character runs longer than the period are
#    considered as non-matches, unless the period is 1.

def correlation_sequence(sequence,period,start=None,end=None):
    if (start == None): start = 0
    if (end   == None): end   = len(sequence)

    prevCh = sequence[start]
    run    = 1
    for ix in xrange(start+1,start+period):
        ch = sequence[ix]
        if (ch != prevCh): run =  1
        else:              run += 1
        prevCh = ch

    corr = [" "] * period
    for ix in xrange(start+period,end):
        rptCh = sequence[ix-period]
        ch    = sequence[ix]
        if (ch != prevCh): run =  1
        else:              run += 1
        if    (ch    in "ACGT") \
          and (ch == rptCh) \
          and ((period == 1) or (run < period)):
            corr += ["+"]
        else:
            corr += ["-"]
        prevCh = ch

    return "".join(corr)


# longest_suitable_run--
#    Find longest run with a good enough rate (of positives).
#
#    We score a "+" as 1-r and anything else as -r.  This is based on the fol-
#    lowing derivation (p is the number of "+"s, n is the number of non-"+"s):
#        p/(p+n) >= r
#        ==> p >= rp + rn
#        ==> (1-r)p - rn >= 0
#
#    We adapt an algorithm from "Programming Pearls", pg. 81 (2000 printing).
#
# $$$ we use the denominator as the threshold, but we really should use the
# $$$ .. numerator, comparing it to minLength*rate
#
# $$$ this needs to account for $$$ this situation:
# $$$   sequence: ACGACGACGACGTTATTATTATTA
# $$$   matches:     +++++++++---+++++++++
# $$$ this is currently considered to be one interval (if rate <= 6/7), but it
# $$$ ought to be two;  we can't just post-process, though, because some other
# $$$ interval might be longer than the longest half of this;  maybe what we
# $$$ need to do is consider matches at distances -P and -2P, or if we match
# $$$ -P but that itself was a mismatch, we should carry the mismatch forward

def longest_suitable_run(seq,minLength,rate):
    maxEndingHere = 0
    maxSoFar      = 0
    start         = None

    for ix in xrange(len(seq)):
        if (seq[ix] == "+"): s = 1-rate
        else:                s = -rate

        if (maxEndingHere+s < 0):
            maxEndingHere = 0
            block         = ix
        else:
            maxEndingHere += s
            if (maxEndingHere >= maxSoFar):
                maxSoFar = maxEndingHere
                start    = block + 1
                end      = ix + 1

    if (start == None) or (end - start < minLength):
        return []
    else:
        return [(start,end)]


# all_suitable_runs--
#    Find all non-overlapping runs with a good enough rate (of positives), and
#    which meet our length threshold.
# $$$ this needs to post-process the intervals, splitting them to account for
# $$$ this situation:
# $$$   sequence: ACGACGACGACGTTATTATTATTA
# $$$   matches:     +++++++++---+++++++++
# $$$ this is currently reported as one interval (if rate <= 6/7), but it
# $$$ ought to be two

def all_suitable_runs(seq,minCorrLength,rate, hammingThreshold):
    
    ################################################################
    # modified by Chen Sun on 07/11/2014
    #
    ################################################################
    
    if hammingThreshold > 0:
        return [(start,end) for (start,end) in dense_intervals(seq,rate,"+","-",blockers=None,minLength=minCorrLength)]
    elif hammingThreshold == 0:
        return [(start,end) for (start,end) in non_negative_intervals(seq, minLength=minCorrLength)]


# find_repeat_element--
#    Find the most plausible repeat element for a run, and nudge the ends of
#    the run if needed.  Note that we will not consider kmers that represent
#    shorter repeats.  For example, we won't report ACTACT as a 6-mer since we
#    consider this to have a shorter period than 6.

def find_repeat_element(hammingThreshold, period,seq,start,end,allowPartials=False):

    if hammingThreshold > 0:
        (kmer,bestD,bestStart,bestEnd) = find_hamming_repeat_element(period,seq,start,end,allowPartials)
        return (kmer,bestD,bestStart,bestEnd)
    # count the number of occurences of each k-mer;  note that we can't
    # reject kmers containing smaller repeats yet, since for a sequence like
    # ACACACACACAAACACACACACACACACAC we must first discover ACACAC as the best
    # 6-mer, and THEN reject it;  if we reject ACACAC while counting, we'd end
    # up reporting something like ACACAA as the best motif 

    if ("element" in debug):
        print >>stderr, "find_repeat_element(%d,%d,%d)" % (period,start,end)

    if ("partial" in debug):
        print period, seq, start, end, allowPartials;
        print seq[start:end]

    kmerToCount = {}
    kmerToFirst = {}
    for ix in xrange(start,end-(period-1)):
        kmer = seq[ix:ix+period]
        if ("N" in kmer): continue
        if (kmer not in kmerToCount):
            kmerToCount[kmer] = 1
            kmerToFirst[kmer] = ix
        else:
            kmerToCount[kmer] += 1
        #if ("element" in debug):
        #    print >>stderr, "    %d: %s" % (ix,kmer)

    # choose the best k-mer;  this is simply the most frequently occurring one,
    # with ties broken by whichever one came first

    kmers = [(-kmerToCount[kmer],kmerToFirst[kmer],kmer) for kmer in kmerToCount]
    if (kmers == []): return (None,None,start,end)
    kmers.sort()

    if ("element" in debug):
        for (count,first,kmer) in kmers:
            print >>stderr, "    %s: %d" % (kmer,-count)

    (count,first,kmer) = kmers[0]
    if (contains_repeat(kmer)): return (None,None,start,end)

    # determine the hamming distance between the run and a simple repeat, for
    # each "plausible" start and end;  we compute the distance for each such
    # interval, and choose the one with the lowest hamming distance;  ties are
    # broken in a deterministic-but-unspecified manner

    bestD = bestStart = bestEnd = None
    ###################################################################################
    # modified by Chen Sun([e-mail]) on 10/18/2013
    #     since we do not allow hamming_distance > 0, which means we do not allow mutation,
    # we do not need this section to produce bestStart and End
    ###################################################################################

    #for (s,e) in plausible_intervals(start,end,period,len(seq),allowPartials=allowPartials):
    #    d = hamming_distance(seq,s,e,kmer)
    #    if (d == None): continue
    #    if (bestD == None) or (d <= bestD):
    #        (bestD,bestStart,bestEnd) = (d,s,e)
    


    bestStart = start

    if(allowPartials):
        bestEnd = end
    elif(not allowPartials):
        bestEnd = start
        pattern = seq[start:start+period]
        if ("partial" in debug):
            print "kmer:", kmer
            if(pattern != kmer):
                print "pattern:", pattern

        while(bestEnd <= end-period):
            bestEnd += period

    # bestD will always be 0, as we do not allow mutation
    bestD = 0
    
    if ("partial" in debug):
        print bestD, bestStart, bestEnd

    ###################################################################################
    # modified by Chen Sun([e-mail]) on 10/10
    # 
    ###################################################################################
    return (kmer,bestD,bestStart,bestEnd)


def find_hamming_repeat_element(period,seq,start,end,allowPartials=False):
    
    # count the number of occurences of each k-mer;  note that we can't
    # reject kmers containing smaller repeats yet, since for a sequence like
    # ACACACACACAAACACACACACACACACAC we must first discover ACACAC as the best
    # 6-mer, and THEN reject it;  if we reject ACACAC while counting, we'd end
    # up reporting something like ACACAA as the best motif 

    if ("element" in debug):
        print >>stderr, "find_repeat_element(%d,%d,%d)" % (period,start,end)

    kmerToCount = {}
    kmerToFirst = {}
    for ix in xrange(start,end-(period-1)):
        kmer = seq[ix:ix+period]
        if ("N" in kmer): continue
        if (kmer not in kmerToCount):
            kmerToCount[kmer] = 1
            kmerToFirst[kmer] = ix
        else:
            kmerToCount[kmer] += 1
        #if ("element" in debug):
        #    print >>stderr, "    %d: %s" % (ix,kmer)

    # choose the best k-mer;  this is simply the most frequently occurring one,
    # with ties broken by whichever one came first

    kmers = [(-kmerToCount[kmer],kmerToFirst[kmer],kmer) for kmer in kmerToCount]
    if (kmers == []): return (None,None,start,end)
    kmers.sort()

    if ("element" in debug):
        for (count,first,kmer) in kmers:
            print >>stderr, "    %s: %d" % (kmer,-count)

    (count,first,kmer) = kmers[0]
    if (contains_repeat(kmer)): return (None,None,start,end)

    # determine the hamming distance between the run and a simple repeat, for
    # each "plausible" start and end;  we compute the distance for each such
    # interval, and choose the one with the lowest hamming distance;  ties are
    # broken in a deterministic-but-unspecified manner

    bestD = bestStart = bestEnd = None

    for (s,e) in plausible_intervals(start,end,period,len(seq),allowPartials=allowPartials):
        d = hamming_distance(seq,s,e,kmer)
        if (d == None): continue
        if (bestD == None) or (d <= bestD):
            (bestD,bestStart,bestEnd) = (d,s,e)

    return (kmer,bestD,bestStart,bestEnd)

# plausible_intervals--
#    Yield all plausible intervals intersecting with a run.  We generate all
#    starts within P bp of the run's start.  For each of these, we either (a) try
#    all ends within P bp of run's end, or (b) trim the new interval to a whole
#    multiple of the period, and report this short interval and the longer
#    interval with one more period appended.  Case (a) allows partial motifs,
#    while case (b) only allows whole motifs.

def plausible_intervals(start,end,period,seqLen,allowPartials=False):

    # generate intervals that allow a partial copy of the motif

    if (allowPartials):
        for candStart in xrange(start-(period-1),start+period):
            if (candStart < 0): continue
            for candEnd in xrange(end-(period-1),end+period):
                if (candEnd > seqLen): continue
                if (candEnd <= candStart+period): continue
                yield (candStart,candEnd)

    # -OR- generate intervals that allow only whole copies of the motif

    else:
        for candStart in xrange(start-(period-1),start+period):
            if (candStart < 0): continue
            candEnd = candStart + ((end-candStart)/period)*period
            yield (candStart,candEnd)
            candEnd += period
            if (candEnd <= seqLen): yield (candStart,candEnd)


# hamming_distance--
#    Determine the hamming distance between the run and a simple repeat.
# $$$ improve this by allowing gaps, and stopping when we reach a threshold

kmerToDiffs = {}  # (this is used for memo-ization)

def hamming_distance(seq,start,end,kmer):
    period = len(kmer)
    if (end < start + period): return None

    wholeEnd = start + ((end-start)/period)*period

    if (kmer not in kmerToDiffs):
        kmerToDiffs[kmer] = { kmer:0 }

    d = 0
    for ix in xrange(start,wholeEnd,period):
        qmer = seq[ix:ix+period]    # same size as the kmer motif
        if (qmer in kmerToDiffs[kmer]):
            d += kmerToDiffs[kmer][qmer]
            continue
        diffs = 0
        for iy in xrange(0,period):
            if (qmer[iy] != kmer[iy]): diffs += 1
        kmerToDiffs[kmer][qmer] = diffs
        d += diffs

    if (end > wholeEnd):
        qmer = seq[wholeEnd:end]    # shorter than the kmer motif
        if (qmer in kmerToDiffs[kmer]):
            d += kmerToDiffs[kmer][qmer]
        else:
            diffs = 0
            for iy in xrange(0,len(qmer)):
                if (qmer[iy] != kmer[iy]): diffs += 1
            kmerToDiffs[kmer][qmer] = diffs
            d += diffs

    return d


# fasta_sequences--
#    Read the fasta sequences from a file.  Note that we convert to upper case,
#    and convert any letter other than ACGT to N.

nonDnaMap = maketrans("BDEFHIJKLMOPQRSUVWXYZ","NNNNNNNNNNNNNNNNNNNNN")

def fasta_sequences(f):
    seqName = None
    seqNucs = None

    for line in f:
        line = line.strip()
        if (line.startswith(">")):
            if (seqName != None):
                yield (seqName,"".join(seqNucs))
            seqName = sequence_name(line)
            seqNucs = []
        elif (seqName == None):
            assert (False), "first sequence has no header"
        else:
            seqNucs += [line]

    if (seqName != None):
        yield (seqName,"".join(seqNucs).upper().translate(nonDnaMap))


# fastq_sequences--
#    Read the fastq sequences from a file.  Note that we convert to upper case,
#    and convert any letter other than ACGT to N.

def fastq_sequences(f):
    lineNum = 0
    for line in f:
        lineNum += 1
        line = line.strip()

        if (lineNum % 4 == 1):
            assert (line.startswith("@")), \
                   "bad read name at line %d" % lineNum
            seqName = line[1:]
            continue

        if (lineNum % 4 == 2):
            seqNucs = line
            continue

        if (lineNum % 4 == 3):
            assert (line.startswith("+")), \
                   "can't understand line %d:\n%s" % (lineNum,line)
            continue

        quals = line
        assert (len(quals) == len(seqNucs)), \
               "length mismatch read vs. qualities at line %d" % lineNum
        yield (seqName,"".join(seqNucs).upper().translate(nonDnaMap),quals)

    assert (lineNum % 4 == 0), \
           "incomplete read at end of file"

def sam_sequences(f):
    lineNum = 0
    for line in f:
        lineNum += 1
        line = line.strip()

        if line.startswith("@"):
            continue

        columns = line.split("\t")
        seqName = columns[0]
        refName = columns[2]
        pre_s = int(columns[3]) - 1
        cigar = columns[5]
        seqNucs = columns[9]
        
        yield (seqName,"".join(seqNucs).upper().translate(nonDnaMap), refName, pre_s, cigar)

# sequence_name--
#    Extract the sequence name from a fasta header.
#    $$$ this may need to be improved $$$

def sequence_name(s):
    s = s[1:].strip()
    if (s == ""): return ""
    else:         return s.split()[0]


# nucleotide_runs--
#    Yield (start,end) for all runs of valid nucleotides in a sequence.

def nucleotide_runs(s):
    runs  = []
    start = None
    for (ix,nuc) in enumerate(s):
        if (nuc in "ACGT"):
            if (start == None):
                start = ix
        else:
            if (start != None):
                yield (start,ix)
                start = None

    if (start != None): yield (start,len(s))


# contains_repeat--
#    Determine whether a short sequence contains a repeated element, such as a
#    6-mer containing a repeated 2-mer (ACACAC) or 3-mer (ACTACT).  The repeat
#    must cover the entire sequence, without mismatches.

def contains_repeat(kmer):
    kmerLength = len(kmer)
    hasRepeat = False
    rptLen = 1
    while (not hasRepeat) and (2 * rptLen <= kmerLength):
        if (kmerLength % rptLen != 0):
            rptLen += 1
            continue
        isRepeat = True
        for i in xrange(rptLen,kmerLength,rptLen):
            if (kmer[i:i+rptLen] != kmer[:rptLen]):
                isRepeat = False
                break
        if (isRepeat):
            hasRepeat = True
            break
        rptLen += 1
    return hasRepeat


# hash108--
#    Return a 108-bit hash "value" of a string

def hash108(s):
    m = md5_new()
    m.update(s)
    return m.hexdigest()[:27]


# float_or_fraction--
#    Convert a string to a number, allowing fractions

def float_or_fraction(s):
    if ("/" in s):
        (numer,denom) = s.split("/",1)
        return float(numer)/float(denom)
    else:
        return float(s)


# int_with_unit--
#    Parse a string as an integer, allowing unit suffixes

def int_with_unit(s):
    if (s.endswith("K")):
        multiplier = 1000
        s = s[:-1]
    elif (s.endswith("M")):
        multiplier = 1000 * 1000
        s = s[:-1]
    elif (s.endswith("G")):
        multiplier = 1000 * 1000 * 1000
        s = s[:-1]
    else:
        multiplier = 1

    try:               return               int(s)   * multiplier
    except ValueError: return int(math.ceil(float(s) * multiplier))


if __name__ == "__main__": main()