update

Author: Bill Majoros

Codon.py

@@ -35,7 +35,7 @@
 
 class Codon:
     def __init__(self,exon,triplet,relative,absolute,isInterrupted):
-        if(not relative): raise Exception("relative is not set")
+        if(relative is None): raise Exception("relative is not set")
         self.triplet=triplet
         self.exon=exon
         self.relativeCoord=relative

CodonIterator.py

@@ -34,15 +34,15 @@
 class CodonIterator:
     def __init__(self,transcript,axisSequenceRef,stopCodons):
         self.transcript=transcript
-        self.stopCodon=selCodons
+        self.stopCodons=stopCodons
 
         # Advance to the exon containing the start codon
         exons=transcript.exons
         strand=transcript.strand
         startCodon=transcript.startCodonAbsolute
         if(startCodon):
             if(len(axisSequenceRef)==0): raise Exception("empty sequence")
-            transcript.loadExonSequences(axisSequenceRef)
+            transcript.loadExonSequences(axisSequenceRef,transcript.exons)
             numExons=len(exons)
             exon=None
             for i in range(0,numExons):
@@ -120,7 +120,7 @@ def nextCodon(self):
                 self.exon=nextExon
         else: # codon was not interrupted by end of exon
             triplet=exonSeq[relative:relative+3]
-            if(isStopCodon[triplet]): self.exon=None
+            if(isStopCodon.get(triplet,False)): self.exon=None
             else:
                 self.relative+=3
                 self.absolute+=3 if strand=="+" else -3

Exon.py

@@ -7,6 +7,7 @@
    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
 
 ######################################################################
 #
@@ -53,6 +54,7 @@
 #   gff=exon.toGff()
 #   begin=exon.getBegin()
 #   end=exon.getEnd()
+#   interval=exon.asInterval()
 #   frame=exon.getFrame()
 #   exon.setFrame(frame)
 #   type=exon.getType()
@@ -87,6 +89,9 @@ def setEnd(self,end):
     def containsCoordinate(self,x):
         return x>=self.begin and x<self.end
 
+    def asInterval(self):
+        return Interval(self.begin,self.end)
+
     def getLength(self):
         return self.end-self.begin
 

FastaReader.py

@@ -23,6 +23,7 @@
 #   reader.doUppercase()
 # Class Methods:
 #   size=FastaReader.getSize(filename)
+#   num=FastaReader.countEntries(filename)
 #   FastaReader.readAll(filename) # returns hash : id->sequence
 #   FastaReader.readAllAndKeepDefs(filename) # returns hash : id->[def,seq]
 #   (defline,seq)=FastaReader.firstSequence(filename)
@@ -94,13 +95,23 @@ def firstSequence(cls,filename):
         reader.close()
         return [defline,seq]
 
+    @classmethod
+    def countEntries(cls,filename):
+        n=0
+        reader=FastaReader(filename)
+        while(True):
+            (defline,seq)=reader.nextSequence()
+            if(not defline): break
+            n+=1
+        return n
+
     @classmethod
     def readAll(cls,filename):
         hash={}
         reader=FastaReader(filename)
         while(True):
             [defline,seq]=reader.nextSequence()
-            if(not defline): break;
+            if(not defline): break
             match=re.search("^\s*>(\S+)",defline)
             if(not match): raise Exception("can't parse defline: "+defline)
             id=match.group(1)
@@ -114,7 +125,7 @@ def readAllAndKeepDefs(cls,filename):
         reader=FastaReader(filename)
         while(True):
             [defline,seq]=reader.nextSequence()
-            if(not defline): break;
+            if(not defline): break
             match=re.search("^\s*>(\S+)",defline)
             if(not match): raise Exception("can't parse defline: "+defline)
             id=match.group(1)

GFF3Parser.py

@@ -150,7 +150,7 @@ def labelStructure(self,root):
         t=root["type"]
         if(t=="gene"):
             obj=self.makeGene(root)
-        elif(t=="transcript"):
+        elif(t=="transcript" or t=="mRNA"):
             obj=self.makeTranscript(root)
         elif(t=="exon" or t=="CDS"):
             obj=self.makeExon(root)
@@ -335,4 +335,4 @@ def test_parser7(filename):
     for gene in genes:
         print(gene.toGff())
 
-test_parser7("/Users/bmajoros/python/test/data/subset.gff3")
+#test_parser7("/Users/bmajoros/python/test/data/subset.gff3")

GffTranscriptReader.py

@@ -197,6 +197,7 @@ def loadGFF_transcript(self,fields,line,transcriptBeginEnd,GFF,
             transcriptId=rex[1]
             transcriptBeginEnd[transcriptId]=[begin,end]
             strand=fields[6]
+            score=fields[5]
             transcriptExtraFields=""
             for i in range(8,len(fields)):
                 transcriptExtraFields+=fields[i]+" "
@@ -211,6 +212,8 @@ def loadGFF_transcript(self,fields,line,transcriptBeginEnd,GFF,
                 transcript.source=fields[1]
                 transcript.setBegin(begin)
                 transcript.setEnd(end)
+            if(transcript.score is None and
+               score!="."): transcript.score=float(score)
             geneId=None
             if(rex.find("genegrp=(\S+)",line)): geneId=rex[1]
             elif(rex.find('gene_id[:=]?\s*\"?([^\s\;"]+)\"?',line)):
@@ -389,26 +392,21 @@ def loadGFF(self,gffFilename):
         transcriptBeginEnd={}
         while(True):
             line=GFF.readline()
-            #print("LINE="+line)
             if(not line): break
             if(not re.search("\S+",line)): continue
             if(re.search("^\s*\#",line)): continue
             fields=line.split("\t") ### \t added 3/24/2017
             if(len(fields)<8): raise Exception("can't parse GTF:"+line)
             if(fields[2]=="transcript"):
-                #print("loading transcript line")
                 self.loadGFF_transcript(fields,line,transcriptBeginEnd,GFF,
                                    transcripts,readOrder,genes)
             elif("UTR" in fields[2] or "utr" in fields[2]):
-                #print("loading UTR")
                 self.loadGFF_UTR(fields,line,transcriptBeginEnd,GFF,
                             transcripts,readOrder,genes)
             elif(fields[2]=="exon"):
-                #print("loading exon: "+line)
                 self.loadGFF_exon(fields,line,transcriptBeginEnd,GFF,
                              transcripts,readOrder,genes)
             elif("CDS" in fields[2] or "-exon" in fields[2]):
-                #print("loading CDS")
                 self.loadGFF_CDS(fields,line,transcriptBeginEnd,GFF,
                             transcripts,readOrder,genes)
         GFF.close()

Interval.py

@@ -24,6 +24,8 @@
 #   union=interval.union(other) # returns an array of intervals
 #   diff=interval.minus(other)  # returns an array of intervals
 #   length=interval.getLength()
+#   begin=interval.getBegin()
+#   end=interval.getEnd()
 #   length=interval.length()
 #   bool=interval.equals($other)
 #   other=interval.clone()
@@ -82,6 +84,12 @@ def length(self):
    def getLength(self):
       return self.length()
 
+   def getBegin(self):
+      return self.begin
+
+   def getEnd(self):
+      return self.end
+
    def equals(self,other):
       return self.begin==other.begin and self.end==other.end
 

SlurmWriter.py

@@ -23,15 +23,15 @@
 #   writer.mem(1500)
 #   writer.threads(16)
 #   writer.setQueue("new,all")
-#   writer.writeArrayScript(slurmDir,jobName,runDir,maxParallel,
+#   writer.writeArrayScript(slurmDir,jobName,maxParallel,
 #                           additional_SBATCH_lines)
 #=========================================================================
 class SlurmWriter:
     """SlurmWriter"""
     def __init__(self):
         self.commands=[]
         self.niceValue=0
-        self.MemValue=None
+        self.memValue=0
         self.threadsValue=0
         self.queue=None
 
@@ -50,8 +50,7 @@ def threads(self,value):
     def setQueue(self,value):
         self.queue=value
 
-    def writeArrayScript(self,slurmDir,jobName,runDir,maxParallel,
-                         moreSBATCH=""):
+    def writeArrayScript(self,slurmDir,jobName,maxParallel,moreSBATCH=""):
         if(moreSBATCH is None): moreSBATCH=""
         if(maxParallel<1): raise Exception("specify maxParallel parameter")
         moreSBATCH=moreSBATCH.rstrip()

Transcript.py

@@ -25,6 +25,7 @@
 #                regardless of strand
 #   startCodonAbsolute : absolute coordinates of start codon, 
 #                        relative to genomic axis
+#   score : float
 #   strand : + or -
 #   exons : pointer to array of Exons (which are actually CDS segments)
 #   UTR : pointer to array of UTR segments
@@ -124,13 +125,16 @@ class Transcript:
     def __init__(self,id,strand=None):
         if(type(id)!=EssexNode): # not an EssexNode
             self.transcriptId=id
+            self.score=None
             self.strand=strand
             self.exons=[]
             self.UTR=[]
             self.rawExons=None
             self.stopCodons={"TAG":1,"TGA":1,"TAA":1}
             self.startCodon=None
             self.extraFields=None
+            self.structureChange=None
+            self.fate=None
         else: # EssexNode
             essex=id
             self.transcriptId=essex.getAttribute("ID")
@@ -140,12 +144,24 @@ def __init__(self,id,strand=None):
             self.end=essex.getAttribute("end")
             self.geneId=essex.getAttribute("gene")
             self.substrate=essex.getAttribute("substrate")
+            score=essex.getAttribute("score")
+            self.score=float(score) if score!="." else None
             self.exons=[]
             self.UTR=[]
             self.rawExons=None
             self.startCodon=None
             self.extraFields=None
             self.stopCodons={"TAG":1,"TGA":1,"TAA":1}
+            changeNode=essex.findChild("structure-change")
+            if(changeNode is not None):
+                changeString=""
+                numElem=changeNode.numElements()
+                for i in range(0,numElem):
+                    elem=changeNode.getIthElem(i)
+                    if(EssexNode.isaNode(elem)): continue
+                    if(changeString!=""): changeString+=" "
+                    changeString+=elem
+                if(len(changeString)>0): self.structureChange=changeString
             exons=self.exons
             UTR=self.UTR
             exonsElem=essex.findChild("exons")
@@ -442,9 +458,14 @@ def toGff(self):
             strand=self.strand
             extra=""
             if(re.search("\S",extraFields)): extra="; "+extraFields
+            change=self.structureChange
+            if(change): extra+="; structure_change \""+change+"\""
+            score=self.score
+            if(score is None): score="."
             gff+=substrate+"\t"+source+"\ttranscript\t"+str(begin)+"\t" \
-                  +str(end)+"\t.\t"+strand+"\t.\ttranscript_id \""+transID+ \
-                  "\"; gene_id \""+geneID+"\""+extra+"\n"
+                  +str(end)+"\t"+str(score)+"\t"+strand+ \
+                  "\t.\ttranscript_id \""+ \
+                  transID+ "\"; gene_id \""+geneID+"\""+extra+"\n"
         for exon in exons: gff+=exon.toGff()
         UTR=self.UTR
         for exon in UTR:
@@ -596,12 +617,13 @@ def trimUTR(self,axisSequenceRef):
         self.recomputeBoundaries()
 
     def getScore(self):
-        exons=self.exons
-        score=0;
-        for exon in exons:
-            exonScore=exon.getScore()
-            if(exonScore!="."): score+=exonScore
-        return score
+        return self.score
+        #exons=self.exons
+        #score=0;
+        #for exon in exons:
+        #    exonScore=exon.getScore()
+        #    if(exonScore!="."): score+=exonScore
+        #return score
 
     def getIntrons(self):
         exons=self.getRawExons()