Commit 079e4381 by Eric Coissac

Patch bug in circular unfolding

parent ffd6006e
...@@ -23,7 +23,7 @@ import sys ...@@ -23,7 +23,7 @@ import sys
from orgasm.multialign import multiAlignReads, consensus # @UnresolvedImport from orgasm.multialign import multiAlignReads, consensus # @UnresolvedImport
from orgasm.assembler import Assembler # @UnresolvedImport from orgasm.assembler import Assembler # @UnresolvedImport
from orgasm.assembler import buildstem # @UnresolvedImport from orgasm.assembler import buildstem # @UnresolvedImport
from orgasm.assembler import tango # @UnresolvedImport from orgasm.assembler import tango, getusedreads # @UnresolvedImport
from orgasm.utils import tags2str from orgasm.utils import tags2str
from time import time from time import time
...@@ -552,6 +552,7 @@ def fillGaps(self,minlink=5, ...@@ -552,6 +552,7 @@ def fillGaps(self,minlink=5,
snp=False, snp=False,
nodeLimit=1000000, nodeLimit=1000000,
onlyLinking=False, onlyLinking=False,
useonce=True,
logger=None): logger=None):
''' '''
...@@ -577,6 +578,7 @@ def fillGaps(self,minlink=5, ...@@ -577,6 +578,7 @@ def fillGaps(self,minlink=5,
__cacheAli = __cacheAli2 __cacheAli = __cacheAli2
__cacheAli2 = set() __cacheAli2 = set()
def isInitial(n): def isInitial(n):
return len(list(assgraph.parentIterator(n[0],edgePredicate=lambda e: 'stemid' in assgraph.getEdgeAttr(*e))))==0 return len(list(assgraph.parentIterator(n[0],edgePredicate=lambda e: 'stemid' in assgraph.getEdgeAttr(*e))))==0
def isTerminal(n): def isTerminal(n):
...@@ -755,11 +757,13 @@ def fillGaps(self,minlink=5, ...@@ -755,11 +757,13 @@ def fillGaps(self,minlink=5,
def insertFragment(self,seq,cycle=1): def insertFragment(self,seq,cycle=1,useonce=True):
index = self.index index = self.index
rsize = index.getReadSize() rsize = index.getReadSize()
readmax=len(index)+1 readmax=len(index)+1
seeds = set() seeds = set()
usedreads = getusedreads()
ireadidE=None ireadidE=None
if len(seq) >= rsize: if len(seq) >= rsize:
...@@ -770,6 +774,7 @@ def insertFragment(self,seq,cycle=1): ...@@ -770,6 +774,7 @@ def insertFragment(self,seq,cycle=1):
for i in range(1,len(seq)-rsize+1): for i in range(1,len(seq)-rsize+1):
coverage = readid[1] coverage = readid[1]
ireadid = readid[0] ireadid = readid[0]
if not useonce or ireadid not in usedreads:
seeds.add(ireadid) seeds.add(ireadid)
if ireadid not in graph: if ireadid not in graph:
...@@ -809,7 +814,7 @@ def insertFragment(self,seq,cycle=1): ...@@ -809,7 +814,7 @@ def insertFragment(self,seq,cycle=1):
readid = readidE readid = readidE
if ireadidE is not None: if ireadidE is not None and (not useonce or ireadid not in usedreads):
seeds.add(ireadidE) seeds.add(ireadidE)
return [(-abs(i),("Consensus",),0) for i in seeds] return [(-abs(i),("Consensus",),0) for i in seeds]
...@@ -1823,7 +1828,7 @@ def unfoldAssembling(self,assgraph,constraints=None, ...@@ -1823,7 +1828,7 @@ def unfoldAssembling(self,assgraph,constraints=None,
# We try to close the graph # We try to close the graph
c = assgraph.minpath(end, c = assgraph.minpath(end,
[first], [first],
edgePredicate=lambda e: 'stemid' in assgraph.getEdgeAttr(*e) edgePredicate=lambda e: assgraph.getEdgeAttr(*e)['class']=="sequence"
) )
i = first i = first
ep = [i] ep = [i]
...@@ -2521,7 +2526,7 @@ def path2sam(self,assgraph,path,identifier="contig",minlink=10,nlength=20,back=2 ...@@ -2521,7 +2526,7 @@ def path2sam(self,assgraph,path,identifier="contig",minlink=10,nlength=20,back=2
reads.append(rl) reads.append(rl)
if rid !=0 and not self.index.isFake(rid): if rid !=0 and not self.index.isFake(rid):
allreads = self.index.getIds(rid) allreads = self.index.getIds(rid)
for rid in allreads: for rid in allreads[2]:
prid = self.index.getPairedRead(rid) prid = self.index.getPairedRead(rid)
rl[rid]=(prid,min(abs(rid),abs(prid))) rl[rid]=(prid,min(abs(rid),abs(prid)))
if rid not in location: if rid not in location:
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment