import os,sys
from optparse import OptionParser
import time
import fastq

maxLen = 200
allbases = ("A", "T", "C", "G");

########################### QualityControl
class QualityControl:
    readLen = 0
    readCount = 0
    counts = {}
    percents = {}
    qualities = [0 for x in xrange(maxLen)]

    def __init__(self):
        for base in allbases:
            self.counts[base] = [0 for x in xrange(maxLen)]
            self.percents[base] = [0.0 for x in xrange(maxLen)]

    def statRead(self, read):
        seq = read[1]
        for i in xrange(len(seq)):
            b = seq[i]
            if b in allbases:
                self.counts[b][i] += 1

    def calcReadLen(self):
        for pos in xrange(maxLen):
            hasData = False
            for base in allbases:
                if self.counts[base][pos]>0:
                    hasData = True
            if hasData == False:
                self.readLen = pos
                break

    def calcPercents(self):
        #calc percents of each base
        for pos in xrange(self.readLen):
            total = 0
            for base in allbases:
                total += self.counts[base][pos]
            for base in allbases:
                self.percents[base][pos] = float(self.counts[base][pos])/float(total)
        
    def statFile(self, filename):
        reader = fastq.Reader(filename)
        #sample up to maxSample reads for stat
        maxSample = 5000
        while True:
            read = reader.nextRead()
            self.readCount += 1
            if read==None or self.readCount>maxSample:
                break
            self.statRead(read)
                
        self.calcReadLen()
        self.calcPercents()

    def autoTrim(self):
        #use (center-5, center+5) as initial good segment        
        center = int(self.readLen/2)
        left = center-5
        right = center+5
        
        threshold = 0.05
        lastStepIsLeft = False
        leftFinished = False
        rightFinished = False
        current = -1
        
        #expand the good segment
        meanPercents = {}
        while not (leftFinished and rightFinished):
            for base in allbases:
                meanPercents[base] = 0.0
                for pos in xrange(left, right):
                    meanPercents[base] += self.percents[base][pos]
                meanPercents[base] /= (right-left)
            
            if leftFinished:
                current = right + 1
                lastStepIsLeft = False
            elif rightFinished:
                current = left - 1
                lastStepIsLeft = True
            elif lastStepIsLeft:
                current = right + 1
                lastStepIsLeft = False
            else:
                current = left - 1
                lastStepIsLeft = True
                                
            percentBias = 0.0
            for base in allbases:
                percentBias += abs(meanPercents[base] - self.percents[base][current])
            
            if percentBias > threshold:
                if lastStepIsLeft:
                    leftFinished = True
                else:
                    rightFinished = True
            else:
                if lastStepIsLeft:
                    left = current
                    if left == 0: leftFinished = True
                else:
                    right = current
                    if right == self.readLen-1: rightFinished = True
                    
        #find the bad segment from front, considering a small window
        #if any in the window is bad, it is bad
        trimFront = left
        window = 3
        for pos in xrange(0, left):
            isGood = True
            for posInWindow in xrange(pos, min(pos+3, self.readLen)):
                percentBias = 0.0
                for base in allbases:
                    percentBias += abs(meanPercents[base] - self.percents[base][posInWindow])
                if percentBias > threshold:
                    isGood = False
            if isGood:    
                trimFront = pos
                break
        #find the bad segment from tail, considering a small window
        #if any in the window is bad, it is bad
        trimTail = right
        for pos in xrange(self.readLen-1, right, -1):
            isGood = True
            for posInWindow in xrange(pos, max(pos-3, 0), -1):
                percentBias = 0.0
                for base in allbases:
                    percentBias += abs(meanPercents[base] - self.percents[base][posInWindow])
                if percentBias > threshold:
                    isGood = False
            if isGood: 
                trimTail = pos
                break
        
        trimFront = min(self.readLen*0.1,trimFront)
        trimTail = min(self.readLen*0.05,self.readLen-1-trimTail)
        # the last base should be definitely trimmed for illumina sequencer output
        trimTail = max(1, trimTail)
        
        return (int(trimFront), int(trimTail))

if __name__  == "__main__":
    qc = QualityControl()
    qc.statFile("R1.fq")
    print(qc.autoTrim())

@@ -65,16 +65,28 @@
 rt fastq
+%0Aimport util
 %0A%0AmaxLen
@@ -88,17 +88,18 @@
 axLen = 
-2
+10
 00%0Aallba
@@ -268,22 +268,110 @@
     
-qualities
+totalQual = %5B0 for x in xrange(maxLen)%5D%0A    totalNum = %5B0 for x in xrange(maxLen)%5D%0A    meanQual
  = %5B0
+.0
  for
@@ -624,16 +624,39 @@
 read%5B1%5D%0A
+        qual = read%5B3%5D%0A
         
@@ -678,24 +678,113 @@
 (len(seq)):%0A
+            self.totalNum%5Bi%5D += 1%0A            self.totalQual%5Bi%5D += util.qualNum(qual%5Bi%5D)%0A
             
@@ -1487,16 +1487,309 @@
 t(total)
+%0A%0A    def calcQualities(self):%0A        for pos in xrange(self.readLen):%0A            self.meanQual%5Bpos%5D = float(self.totalQual%5Bpos%5D)/float(self.totalNum%5Bpos%5D)%0A%0A    def qc(self):    %0A        self.calcReadLen()%0A        self.calcPercents()%0A        self.calcQualities()%0A        print(self.meanQual)
 %0A       
@@ -1935,16 +1935,20 @@
 e = 5000
+0000
 %0A       
@@ -2142,77 +2142,24 @@
 ad)%0A
-                %0A        self.calcReadLen()%0A        self.calcPercents
+%0A        self.qc
 ()%0A%0A

add qual stat