update

Author: Bill Majoros

CigarString.py

@@ -28,19 +28,30 @@
 #   (numMatches,numMismatches)=cigar.longestMatchStats(seq1,seq2) 
 #       # ^ Returns none if no match; must call computeIntervals() first!
 #   L=cigar.totalAlignmentLength()
+#   cigar.setOps(ops)
 #=========================================================================
 class CigarString:
     """CigarString parses CIGAR strings (alignments)"""
     def __init__(self,cigar):
         self.ops=self.parse(cigar)
 
+    def setOps(self,ops):
+        self.ops=ops
+
     def matchesByLength(self):
         matches=[]
         for op in self.ops:
             if(op.getOp() in ("M","=","X")): matches.append(op)
         matches.sort(key=lambda x: -x.getLength())
         return matches
 
+    def countIndelBases(self):
+        N=0
+        for op in self.ops:
+            c=op.getOp()
+            if(c=="I" or c=="D"): N+=1
+        return N
+
     def totalAlignmentLength(self):
         L=0
         for op in self.ops:

SamHSP.py

@@ -0,0 +1,85 @@
+#=========================================================================
+# This is OPEN SOURCE SOFTWARE governed by the Gnu General Public
+# License (GPL) version 3, as described at www.opensource.org.
+# Copyright (C)2016 William H. Majoros (martiandna@gmail.com).
+#=========================================================================
+from __future__ import (absolute_import, division, print_function,
+   unicode_literals, generators, nested_scopes, with_statement)
+from builtins import (bytes, dict, int, list, object, range, str, ascii,
+   chr, hex, input, next, oct, open, pow, round, super, filter, map, zip)
+from CigarString import CigarString
+from Interval import Interval
+
+#=========================================================================
+# Attributes:
+#   cigar : CigarString
+#   refName : string
+#   readInterval : Interval
+#   refInterval : Interval
+#   score : float
+#   rec : the SamRecord this HSP came from
+# Instance Methods:
+#   hsp=SamHSP(rec) # rec is a SamRecord
+#   cigar=hsp.getCigar() # returns CigarString object
+#   refName=hsp.getRefName()
+#   boolean=hsp.overlapsOnRead(otherHSP)
+#   boolean=hsp.overlapsOnRef(otherHSP)
+#   interval=hsp.getReadInterval()  # returns Interval object
+#   interval=hsp.getRefInterval() # returns Interval object
+#   hsp.computeScore()
+#   score=hsp.getScore()
+# Private Methods:
+#   self.computeIntervals()
+# Class Methods:
+#   none
+#=========================================================================
+class SamHSP:
+    """SamHSP"""
+    def __init__(self,rec):
+        self.cigar=rec.getCigar()
+        self.refName=rec.getRefName()
+        self.rec=rec
+        self.computeIntervals()
+        self.score=None
+
+    def getScore(self):
+        return self.score
+
+    def computeScore(self):
+        mismatches=self.rec.countMismatches()
+        matches=self.cigar.totalAlignmentLength()-mismatches
+        indelBases=self.cigar.countIndelBases()
+        numerator=matches
+        denominator=1+mismatches+indelBases
+        self.score=float(numerator)/float(denominator)
+
+    def overlapsOnRead(self,other):
+        return self.readInterval.overlaps(other.readInterval)
+
+    def overlapsOnRef(self,other):
+        return self.refInterval.overlaps(other.refInterval)
+
+    def getCigar(self):
+        return self.cigar
+
+    def getRefName(self):
+        return self.refName
+
+    def getReadInterval(self):
+        return self.readInterval
+
+    def getRefInterval(self):
+        return self.refInterval
+
+    def computeIntervals(self):
+        cigar=self.cigar
+        n=cigar.length()
+        firstOp=cigar[0]
+        lastOp=cigar[n-1]
+        self.readInterval=Interval(firstOp.getQueryInterval().getBegin(),
+                                   lastOp.getQueryInterval().getEnd())
+        self.refInterval=Interval(firstOp.getRefInterval().getBegin(),
+                                  lastOp.getRefInterval().getEnd())
+            
+
+

SamHspClusterer.py

@@ -0,0 +1,41 @@
+#=========================================================================
+# This is OPEN SOURCE SOFTWARE governed by the Gnu General Public
+# License (GPL) version 3, as described at www.opensource.org.
+# Copyright (C)2016 William H. Majoros (martiandna@gmail.com).
+#=========================================================================
+from __future__ import (absolute_import, division, print_function,
+   unicode_literals, generators, nested_scopes, with_statement)
+from builtins import (bytes, dict, int, list, object, range, str, ascii,
+   chr, hex, input, next, oct, open, pow, round, super, filter, map, zip)
+
+#=========================================================================
+# Attributes:
+#   
+# Instance Methods:
+#   clusterer=SamHspClusterer()
+# Class Methods:
+#   clustered=SamHspClusterer.cluster(HSPs)
+#=========================================================================
+class SamHspClusterer:
+    """SamHspClusterer"""
+    def __init__(self):
+        pass
+
+    @classmethod
+    def cluster(cls,raw):
+        HSPs=[x for x in raw]
+        HSPs.sort(key=lambda hsp: hsp.getReadInterval().getBegin())
+        for hsp in HSPs: hsp.computeScore()
+        n=len(HSPs)
+        i=0
+        while(i<n-1):
+            while(HSPs[i].overlapsOnRead(HSPs[i+1])):
+                if(HSPs[i].getScore()<HSPs[i+1].getScore()):
+                    del HSPs[i]
+                else:
+                    del HSPs[i+1]
+                n-=1
+                if(i>=n-1): break
+            i+=1
+        return HSPs
+

SamHspFactory.py

@@ -0,0 +1,52 @@
+#=========================================================================
+# This is OPEN SOURCE SOFTWARE governed by the Gnu General Public
+# License (GPL) version 3, as described at www.opensource.org.
+# Copyright (C)2016 William H. Majoros (martiandna@gmail.com).
+#=========================================================================
+from __future__ import (absolute_import, division, print_function,
+   unicode_literals, generators, nested_scopes, with_statement)
+from builtins import (bytes, dict, int, list, object, range, str, ascii,
+   chr, hex, input, next, oct, open, pow, round, super, filter, map, zip)
+from SamHSP import SamHSP
+from SamRecord import SamRecord
+from CigarString import CigarString
+
+#=========================================================================
+# Attributes:
+#   keepOps : set of string
+# Instance Methods:
+#   factory=SamHspFactory()
+#   HSPs=factory.makeHSPs(reads)
+# Private Methods:
+#   cigar=self.processCigar(cigar)
+# Class Methods:
+#   none
+#=========================================================================
+class SamHspFactory:
+    """SamHspFactory"""
+    def __init__(self):
+        self.keepOps=set(["M","I","D","=","X"])
+
+    def makeHSPs(self,reads):
+        HSPs=[]
+        for read in reads:
+            cigar=read.getCigar()
+            cigar.computeIntervals(read.getRefPos())
+            cigar=self.processCigar(cigar)
+            hsp=SamHSP(read)
+            HSPs.append(hsp)
+        return HSPs
+
+    def processCigar(self,cigar):
+        keepOps=self.keepOps
+        keep=[]
+        n=cigar.length()
+        for i in range(n):
+            op=cigar[i]
+            #if(op.getOp()=="S" and i>0 and i<n-1):
+            #    raise Exception("Nested soft mask detected in cigar string")
+            if(op.getOp() in keepOps):
+                keep.append(op)
+        new=CigarString("")
+        new.setOps(keep)
+        return new

SamPairedRead.py

@@ -35,6 +35,13 @@ def totalAlignedLength(self):
         L2=self.read2.getCigar().totalAlignmentLength()
         return L1+L2
 
+    def countIndelBases(self):
+        cigar1=self.read1.getCigar()
+        cigar2=self.read2.getCigar()
+        N1=cigar1.countIndelBases()
+        N2=cigar2.countIndelBases()
+        return N1+N2
+
     def numMismatches(self):
         mis1=self.read1.countMismatches()
         mis2=self.read2.countMismatches()
@@ -46,6 +53,14 @@ def matchProportion(self):
         score=float(L-mis)/float(L)
         return score
 
+    def computeScore(self):
+        mismatches=self.numMismatches()
+        matches=self.totalAlignedLength()-mismatches
+        numerator=matches
+        indelBases=self.countIndelBases()
+        denominator=1+mismatches+indelBases
+        score=float(numerator)/float(denominator)
+        return score
 
 
 

SamPairedReadStream.py

@@ -14,10 +14,12 @@
 # Attributes:
 #   reader : SamReader
 #   dedup : boolean
-#   buffer : SamRecord
+#   bufferedRec : SamRecord
+#   bufferedPair : SamPairedRead
 # Instance Methods:
 #   stream=SamPairedReadStream(filename,dedup=True)
 #   pair=stream.nextPair() # returns SamPairedRead
+#   readGroup=stream.nextGroup() # returns array of SamPairedRead
 # Class Methods:
 #   none
 #=========================================================================
@@ -26,14 +28,37 @@ class SamPairedReadStream:
     def __init__(self,filename,dedup=True):
         self.reader=SamReader(filename)
         self.dedup=dedup
-        self.buffer=None
+        self.bufferedRec=None
+        self.bufferedPair=None
+
+    def nextGroup(self):
+        group=[]
+        readID=None
+        while(True):
+            pair=None
+            if(self.bufferedPair is not None):
+                pair=self.bufferedPair
+                self.bufferedPair=None
+            else: pair=self.nextPair()
+            if(pair is None): return None
+            if(readID is None):
+                group.append(pair)
+                readID=pair.getID()
+                continue
+            thisID=pair.getID()
+            if(thisID==readID):
+                group.append(pair)
+                continue
+            self.bufferedPair=pair
+            break
+        return group
 
     def nextPair(self):
         while(True):
             rec=None
-            if(self.buffer is not None): 
-                rec=self.buffer
-                self.buffer=None
+            if(self.bufferedRec is not None): 
+                rec=self.bufferedRec
+                self.bufferedRec=None
             else: rec=self.reader.nextSequence()
             if(rec is None): return None
             if(rec.flag_unmapped()): continue
@@ -45,7 +70,7 @@ def nextPair(self):
             if(rec.flag_unmapped()): continue
             if(self.dedup and rec.flag_PCRduplicate()): continue
             if(rec.flag_firstOfPair()): 
-                self.buffer=rec
+                self.bufferedRec=rec
                 continue
             read2=rec
             return SamPairedRead(read1,read2)