diff --git a/ngs_mapper/runsample.py b/ngs_mapper/runsample.py index 32d00771..2ec053c6 100644 --- a/ngs_mapper/runsample.py +++ b/ngs_mapper/runsample.py @@ -427,8 +427,9 @@ def select_keys(d, keys): # Filter on index quality and Ns # Mapping + cmd_args['tmp_bam'] = os.path.join( tdir, 'tmp.bam' ) with open(bwalog, 'wb') as blog: - cmd = 'run_bwa_on_samplename {trim_outdir} {reference} -o {bamfile}' + cmd = 'run_bwa_on_samplename {trim_outdir} {reference} -o {tmp_bam}' if cmd_args['config']: cmd += ' -c {config}' p = run_cmd( cmd.format(**cmd_args), stdout=blog, stderr=subprocess.STDOUT ) @@ -440,6 +441,13 @@ def select_keys(d, keys): logger.critical( "{0} failed to complete sucessfully. Please check the log file {1} for more details".format(cmd,bwalog) ) sys.exit(1) + tcmd = "samtools view -F 2048 -b {tmp_bam} > {bamfile}" + p = run_cmd( tcmd.format(**cmd_args), stdout=lfile, stderr=subprocess.STDOUT ) + r = p.wait() + if r != 0: + logger.critical( "{0} did not exit sucessfully".format(cmd.format(**cmd_args)) ) + rets.append( r ) + # Tag Reads cmd = 'tagreads {bamfile} -CN {CN}' if cmd_args['config']: @@ -450,15 +458,19 @@ def select_keys(d, keys): logger.critical( "{0} did not exit sucessfully".format(cmd.format(**cmd_args)) ) rets.append( r ) # Variant Calling + from ngs_mapper.config import load_config + bc_cfg = load_config(cmd_args['config'])['base_caller'] + lofreq = bc_cfg['lofreq']['default'] + lofreq_options = bc_cfg['lofreq_options']['default'] or '' if lofreq: cmd1 = 'lofreq call -f {reference} {bamfile} -o {vcf} ' + lofreq_options r1 = run_cmd( cmd1.format(**cmd_args), stdout=lfile, stderr=subprocess.STDOUT ).wait() if r1 != 0: logger.critical( "{0} did not exit sucessfully".format(cmd1.format(**cmd_args)) ) - - cmd2 = 'lf_consensus --ref {reference} --vcf {vcf} --minbq 30 --mind 10 --bam {bamfile} > {consensus}' - r2 = run_cmd( cmd2.format(**cmd_args), stdout=lfile, stderr=subprocess.STDOUT ).wait() - if r2 != 0: logger.critical( "{0} did not exit sucessfully".format(cmd2.format(**cmd_args)) ) - rets += [r1, r2] + with open(cmd_args['consensus'], 'w') as cons_out: + cmd2 = 'lf_consensus --ref {reference} --vcf {vcf} --minbq 30 --mind 10 --bam {bamfile}' + r2 = run_cmd( cmd2.format(**cmd_args), stdout=cons_out, stderr=subprocess.STDOUT ).wait() + if r2 != 0: logger.critical( "{0} did not exit sucessfully".format(cmd2.format(**cmd_args)) ) + rets += [r1, r2] else: cmd = 'base_caller {bamfile} {reference} {vcf} -minth {minth}' if cmd_args['config']: diff --git a/requirements-pip.txt b/requirements-pip.txt index 6da78f48..a27d99fe 100644 --- a/requirements-pip.txt +++ b/requirements-pip.txt @@ -10,3 +10,6 @@ Counter==1.0.0 schema==0.4.0 PyVCF==0.6.6 python-dateutil==2.1 +docopt +plumbum +toolz diff --git a/setup.py b/setup.py index 9c727764..d1ca0d37 100644 --- a/setup.py +++ b/setup.py @@ -40,6 +40,7 @@ 'trim_reads = ngs_mapper.trim_reads:main', 'vcf_consensus = ngs_mapper.vcf_consensus:main', 'vcf_diff = ngs_mapper.vcf_diff:main', + 'lf_consensus = ngs_mapper.lofreq_consensus:main' ] }, package_data = {