Commit b330de3e by Eric Coissac

Adds options to control minread estimate

parent 359c026f
......@@ -28,11 +28,13 @@ default_config = { 'reformat' : False,
'smallbranches' : None,
'lowcomplexity' : False,
'snp' : False,
'assmax' : 500000,
'testrun' : 15000,
'assmax' : 500000, # Maximum size of the assembling
'testrun' : 15000, # length of to assembling to estimate coverage
'clean' : False,
'forceseeds' : False,
'maxfillgaps' : -1
'maxfillgaps' : -1,
'covratio' : 3,
'covratiofill' : 12
}
def addOptions(parser):
......@@ -65,6 +67,18 @@ def addOptions(parser):
default=None,
help='the expected sequencing coverage [default: <estimated>]')
parser.add_argument('--coverage-ratio', dest='buildgraph:covratio',
type=int,
action='store',
default=None,
help='the ratio between the expected and the minimum coverage [default 3]')
parser.add_argument('--fillgaps-ratio', dest='buildgraph:covratio',
type=int,
action='store',
default=None,
help='the ratio between the expected and the minimum coverage during gap filling [default 12]')
parser.add_argument('--minratio', dest='buildgraph:minratio',
type=float, action='store',
default=None,
......@@ -200,13 +214,12 @@ def addOptions(parser):
def estimateMinRead(index,minoverlap,coverage):
def estimateMinRead(index,minoverlap,coverage,minreadcor):
MINREAD=10
MINREADCOR=3
MINOVERLAP=50
minread = (index.getReadSize() - minoverlap) * coverage / index.getReadSize() / MINREADCOR
minread = (index.getReadSize() - minoverlap) * coverage / index.getReadSize() / minreadcor
if minread < MINREAD:
minoverlap = index.getReadSize() - (MINREAD * MINREADCOR * index.getReadSize() / coverage)
minoverlap = index.getReadSize() - (MINREAD * minreadcor * index.getReadSize() / coverage)
minread = MINREAD
if minoverlap< MINOVERLAP:
minread = MINREAD
......@@ -223,6 +236,8 @@ def run(config):
coverageset=config['buildgraph']['coverage'] is not None
snp=config['buildgraph']['snp']
assmax = config['buildgraph']['assmax']*2
covratio = config['buildgraph']['covratio']
covratiofill = config['buildgraph']['covratiofill']
......@@ -269,7 +284,9 @@ def run(config):
if coverageset:
logger.info('Coverage forced by user at : %d' % coverage)
minread,minoverlap = estimateMinRead(r,config['buildgraph']['minoverlap'],coverage)
minread,minoverlap = estimateMinRead(r,config['buildgraph']['minoverlap'],
coverage,
covratio)
else:
##########################
......@@ -281,7 +298,9 @@ def run(config):
##########################
if config['buildgraph']['minread'] is None:
minread,minoverlap = estimateMinRead(r,config['buildgraph']['minoverlap'],coverage)
minread,minoverlap = estimateMinRead(r,config['buildgraph']['minoverlap'],
coverage,
covratio)
minread//=4
if minread<5:
minread=5
......@@ -328,7 +347,9 @@ def run(config):
score,length,ecoverage = coverageEstimate(asm,x,r) # @UnusedVariable
if not coverageset:
coverage = ecoverage
minread,minoverlap = estimateMinRead(r,config['buildgraph']['minoverlap'],coverage)
minread,minoverlap = estimateMinRead(r,config['buildgraph']['minoverlap'],
coverage,
covratio)
logger.info("coverage estimated : %dx based on %d bp (minread: %d)" %(coverage,length/2,minread))
else:
......@@ -397,12 +418,16 @@ def run(config):
# according to the minread option estimate it from coverage or use the specified value
if config['buildgraph']['minread'] is None:
minread,minoverlap = estimateMinRead(r,config['buildgraph']['minoverlap'],coverage)
minread,minoverlap = estimateMinRead(r,config['buildgraph']['minoverlap'],
coverage,
covratio)
if minread<5:
minread=5
else:
minread,minoverlap = estimateMinRead(r,config['buildgraph']['minoverlap'],coverage)
minread,minoverlap = estimateMinRead(r,config['buildgraph']['minoverlap'],
coverage,
covratio)
minread=config['buildgraph']['minread']
......@@ -494,13 +519,8 @@ def run(config):
# if coverage < 30:
# sys.exit()
if config['buildgraph']['minread'] is None:
minread,minoverlap = estimateMinRead(r,config['buildgraph']['minoverlap'],coverage)
minread/=4
if minread<5:
minread=5
logger.info("coverage estimated : %d based on %d bp (minread: %d)" %(coverage,length/2,minread))
logger.info("coverage estimated : %d based on %d bp" %(coverage,length/2))
meanlength,sdlength = estimateFragmentLength(asm)
......@@ -539,6 +559,16 @@ def run(config):
#
###################################################
if config['buildgraph']['minread'] is None:
logger.info("Estimating minread")
minread,minoverlap = estimateMinRead(r,config['buildgraph']['minoverlap'],
coverage,
covratiofill)
if minread<5:
minread=5
logger.info("minread set to: %d" % minread)
delta = 1
maxfillgaps=config['buildgraph']['maxfillgaps']
fillcycle=0
......
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