update

Author: Bill Majoros

CigarOp.py

@@ -7,16 +7,50 @@
    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 Interval import Interval
+
+ADVANCE_QUERY=set(["M","I","S","=","X"])
+ADVANCE_REF=set(["M","D","N","=","X"])
 
 #=========================================================================
 # Attributes:
-#   op : M/I/D/S
 #   length : integer
+#   interval1 : Interval (in sequence 1 = query)
+#   interval2 : Interval (in sequence 2 = reference)
+#   op : M/I/D/S:
+#                                                              consumes
+#                                                              query ref
+#     M 0 alignment match (can be a sequence match or mismatch) yes yes
+#     I 1 insertion to the reference                            yes  no
+#     D 2 deletion from the reference                            no yes
+#     N 3 skipped region from the reference                      no yes
+#     S 4 soft clipping (clipped sequences present in SEQ)      yes  no
+#     H 5 hard clipping (clipped sequences NOT present in SEQ)   no  no
+#     P 6 padding (silent deletion from padded reference)        no  no
+#     = 7 sequence match                                        yes yes
+#     X 8 sequence mismatch                                     yes yes
 # Instance Methods:
 #   op=CigarOp("M",135)
+#   bool=op.advanceInQuery()   # matches, insertions, etc.
+#   bool=op.advanceInRef() # matches, deletions, etc.
+#   op=op.getOp()
+#   L=op.getLength()
 #=========================================================================
 class CigarOp:
     def __init__(self,op,L):
         self.op=op
         self.length=L
+        self.interval=None
+
+    def getOp(self): 
+        return self.op
+
+    def getLength(self):
+        return self.length
+
+    def advanceInQuery(self):
+        return self.op in ADVANCE_QUERY
+
+    def advanceInRef(self):
+        return self.op in ADVANCE_REF
 

CigarString.py

@@ -10,6 +10,7 @@
 from Rex import Rex
 rex=Rex()
 from CigarOp import CigarOp
+from Interval import Interval
 
 #=========================================================================
 # Attributes:
@@ -19,12 +20,28 @@
 #   bool=cigar.completeMatch()
 #   numOps=cigar.length()
 #   cigarOp=cigar[i] # returns a CigarOp object
+#   str=cigar.toString()
+#   cigar.computeIntervals(refPos)
 #=========================================================================
 class CigarString:
     """CigarString parses CIGAR strings (alignments)"""
     def __init__(self,cigar):
         self.ops=self.parse(cigar)
 
+    def computeIntervals(self,refPos):
+        ops=self.ops
+        n=len(ops)
+        begin1=0; begin2=refPos
+        for i in range(n):
+            op=ops[i]
+            L=op.getLength()
+            end1=begin1; end2=begin2
+            if(op.advanceInQuery()): end1+=L
+            if(op.advanceInRef()): end2+=L
+            op.interval1=Interval(begin1,end1)
+            op.interval2=Interval(begin2,end2)
+            begin1=end1; begin2=end2
+
     def length(self):
         return len(self.ops)
 

DataFrame.py

@@ -20,6 +20,7 @@
 #   df=DataFrame()
 #   rowNames=df.getRowNames()
 #   colNames=df.getColumnNames()
+#   df.addRow(DataFrameRow)
 #   n=df.nrow()
 #   n=df.ncol()
 #   row=df[index]
@@ -29,12 +30,18 @@
 #   header=df.getHeader()
 #   df.hashRowNames()
 #   df.hashColNames()
+#   row=df.getRowI(i)
+#   col=df.getColI(i)
 #   row=df.getRow(rowName) # call hashRowNames() first!
 #   col=df.getColumn(columnName) # call hashColNames() first!
 #   bool=df.rowExists(rowName) # call hashRowNames() first!
 #   bool=df.columnExists(colName) # call hashColNames() first!
+#   newDataFrame=df.subsetColumns(colIndices)
+#   idx=df.addColumn(colName,defaultValue) # returns index of new column
+#   df.print(handle)
+#   array=df.toDataArray()
 # Class methods:
-#   df=DataFrame.readTable(filename,hasHeader=True,hasRowNames=True)
+#   df=DataFrame.readTable(filename,header=False,rowNames=False)
 #=========================================================================
 
 class DataFrame:
@@ -44,6 +51,39 @@ def __init__(self):
       self.rowHash=None
       self.colHash=None
 
+   def addRow(self,row):
+      self.matrix.append(row)
+
+   def toDataArray(self):
+      array=[]
+      for row in self.matrix:
+         array.append(row.values)
+      return array
+
+   def print(self,handle):
+      print("\t".join(self.header),file=handle)
+      for row in self.matrix: row.print(handle)
+
+   def addColumn(self,name,defaultValue):
+      colIndex=len(self.header)
+      self.header.append(name)
+      for row in self.matrix:
+         row.append(defaultValue)
+      return colIndex
+
+   def subsetColumns(self,colIndices):
+      newDF=DataFrame()
+      header=self.header
+      newHeader=newDF.header
+      for i in colIndices: newHeader.append(header[i])
+      for i in range(self.nrow()):
+         row=self[i]
+         newRow=DataFrameRow()
+         newRow.rename(row.getLabel())
+         for j in colIndices: newRow.values.append(row[j])
+         newDF.matrix.append(newRow)
+      return newDF
+
    def rowExists(self,rowName):
       if(self.rowHash is None): raise Exception("call hashRowNames() first")
       return self.rowHash.get(rowName,None) is not None
@@ -59,7 +99,16 @@ def getRowNames(self):
       return names
 
    def getColumnNames(self):
-      return header
+      return self.header
+
+   def getRowI(self,rowIndex):
+      return self.matrix[rowIndex]
+
+   def getColumnI(self,colIndex):
+      column=DataFrameRow()
+      for row in self.matrix:
+         column.values.append(row[colIndex])
+      return column
 
    def getRow(self,rowName):
       if(self.rowHash is None): raise Exception("call hashRowNames() first")
@@ -75,6 +124,7 @@ def getColumn(self,colName):
       column.label=colName
       for row in self.matrix:
          colum.values.append(row[colIndex])
+      return column
 
    def hashRowNames(self):
       h=self.rowHash={}
@@ -108,17 +158,17 @@ def toFloat(self):
       for row in self.matrix: row.toFloat()
 
    @classmethod
-   def readTable(cls,filename,hasHeader=True,hasRowNames=True):
+   def readTable(cls,filename,header=False,rowNames=False):
       df=DataFrame()
       with open(filename,"rt") as IN:
-         if(hasHeader):
+         if(header):
             df.header=IN.readline()
-            df.header=df.header.rstrip().split()
+            df.header=df.header.rstrip().split("\t")
          for line in IN:
-            fields=line.rstrip().split()
+            fields=line.rstrip().split("\t")
             if(len(fields)<1): continue
             label=""
-            if(hasRowNames):
+            if(rowNames):
                label=fields[0]
                fields=fields[1:]
             row=DataFrameRow()

DataFrameRow.py

@@ -21,6 +21,8 @@
 #   n=row.length()
 #   row.toInt()
 #   row.toFloat()
+#   row.append(value)
+#   row.print(handle)
 #=========================================================================
 
 class DataFrameRow:
@@ -31,6 +33,15 @@ def __init__(self):
    def __getitem__(self,i):
       return self.values[i]
 
+   def __setitem__(self,i,value):
+      self.values[i]=value
+
+   def print(self,handle):
+      print(self.label+"\t","\t".join([str(x) for x in self.values]),sep="")
+
+   def append(self,value):
+      self.values.append(value)
+
    def length(self):
       return len(self.values)
 

Interval.py

@@ -15,6 +15,7 @@
 #   end
 # Methods:
 #   i=Interval(begin,end)
+#   s=interval.toString()
 #   print(file=STDOUT)
 #   bool=interval.overlaps(other)
 #   bool=interval.contains(position)
@@ -49,6 +50,9 @@ def __init__(self,begin=0,end=0):
    def print(self,file=sys.stdout):
       print("(",self.begin,",",self.end,")",sep="",end="",file=file)
 
+   def toString(self):
+      return "("+str(self.begin)+","+str(self.end)+")"
+
    def overlaps(self,other):
       return self.begin<other.end and other.begin<self.end
 

SamReader.py

@@ -18,7 +18,7 @@
 #   fh : file handle
 # Instance Methods:
 #   reader=SamReader(filename)
-#   [ID,seq,qual]=reader.nextSequence() # returns None at EOF
+#   samRecord=reader.nextSequence() # returns None at EOF
 #   reader.close()
 # Class Methods:
 #=========================================================================
@@ -43,7 +43,7 @@ def nextSequence(self):
         if(len(fields)<11): raise Exception("can't parse sam line: "+line)
         (ID,flags,refName,refPos,mapQual,cigar,rnext,pnext,templateLen,
          seq,qual)=fields[:11]
-        refPos=int(refPos)
+        refPos=int(refPos)-1 # convert 1-based to 0-based
         flags=int(flags)
         CIGAR=CigarString(cigar)
         rec=SamRecord(ID,refName,refPos,CIGAR,seq,flags)

SamRecord.py

@@ -15,8 +15,26 @@
 #   refPos = position in reference where alignment begins
 #   CIGAR = CigarString
 #   seq = read sequence
+#   flags = bitfield
 # Instance Methods:
-#   rec=SamReader(ID,refName,refPos,cigar,seq)
+#   rec=SamReader(ID,refName,refPos,cigar,seq,flags)
+#   ID=rec.getID()
+#   cigar=rec.getCigar()
+#   seq=rec.getSequence()
+#   refName=rec.getRefName()
+#   refPos=rec.getRefPos()
+#   bool=rec.flag_hasMultipleSegments()
+#   bool=rec.flag_properlyAligned()
+#   bool=rec.flag_unmapped()
+#   bool=rec.flag_nextSegmentUnmapped()
+#   bool=rec.flag_revComp()
+#   bool=rec.flag_nextSegmentRevComp()
+#   bool=rec.flag_firstOfPair()
+#   bool=rec.flag_secondOfPair()
+#   bool=rec.flag_secondaryAlignment()
+#   bool=rec.flag_failedFilters()
+#   bool=rec.flag_PCRduplicate()
+#   bool=rec.flag_supplAlignment()
 # Class Methods:
 #=========================================================================
 class SamRecord:
@@ -29,6 +47,21 @@ def __init__(self,ID,refName,refPos,CIGAR,seq,flags):
         self.seq=seq
         self.flags=flags
 
+    def getRefName(self):
+        return self.refName
+
+    def getRefPos(self):
+        return self.refPos
+    
+    def getCigar(self):
+        return self.CIGAR
+
+    def getID(self):
+        return self.ID
+
+    def getSequence(self):
+        return self.seq
+
     def flag_hasMultipleSegments(self):
         return bool(self.flags & 0x1)