Skip to content

Commit

Permalink
new --ungapped alignment opt; fixed tmp path use to avoid retaining m…
Browse files Browse the repository at this point in the history
…akeblastdb files
  • Loading branch information
chrisgulvik committed Oct 2, 2020
1 parent 52b638d commit 05f5ad2
Showing 1 changed file with 18 additions and 10 deletions.
28 changes: 18 additions & 10 deletions c-SSTAR
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
#!/usr/bin/env python

__version__ = '2.0.0'
__version__ = '2.1.0'

import argparse
import gzip
Expand Down Expand Up @@ -52,6 +52,9 @@ def parse_args():
'--min-identity', type=require_float_1_to_100,
metavar='FLOAT', default=95,
help='minimum alignment nucleotide identity [%(default)s]')
aln.add_argument(
'--ungapped', action='store_true', default=False,
help='store ungapped-only alignments for BLASTn [off]')
opt = parser.add_argument_group('I/O Options')
opt.add_argument(
'--basename', type=str, metavar='STR',
Expand All @@ -72,6 +75,9 @@ def parse_args():
'--report', default=None, metavar='FILE',
help='summary tab-delimited report [stdout]')
msc = parser.add_argument_group('Misc Options')
msc.add_argument(
'--debug', default=False, action='store_true',
help='turn on debugging')
msc.add_argument(
'-h', '--help', action='help',
help='show this help message and exit')
Expand Down Expand Up @@ -291,7 +297,7 @@ def main():
'command -v blastn', shell=True).rstrip())
logging.info(subprocess.check_output(
'blastn -version | tail -n 1', shell=True).rstrip())
tmp_query_idx = os.path.join(outdir, base_genome)
tmp_query_idx = os.path.join(tmp, base_genome)
sys_call('makeblastdb -in {} -out {} -dbtype nucl'.format(
infile_query, tmp_query_idx))
for x in ['.nin', '.nsq', '.nhr']:
Expand All @@ -302,14 +308,16 @@ def main():
' search\n{} file empty\n'.
format(os.path.join(outdir, base_genome + x)))
sys.exit(1)
sys_call(
'blastn -task blastn -query {} -db {} -out {} -evalue 1e-5'
' -max_target_seqs 1 -perc_identity {} -culling_limit 1'
' -outfmt "6 qseqid sseqid pident length mismatch gaps qstart qend'
' sstart send evalue bitscore qlen slen sseq" -num_threads {}'.
format(infile_db, tmp_query_idx,
os.path.join(outdir, base_genome + '.blastn.tsv'),
min_identity, args.cpus))
cmd = ('blastn -task blastn -query {} -db {} -out {} -evalue 1e-5'
' -max_target_seqs 1 -perc_identity {} -culling_limit 1'
' -outfmt "6 qseqid sseqid pident length mismatch gaps qstart qend'
' sstart send evalue bitscore qlen slen sseq" -num_threads {}'.
format(infile_db, tmp_query_idx,
os.path.join(outdir, base_genome + '.blastn.tsv'),
min_identity, args.cpus))
if args.ungapped:
cmd += ' -ungapped'
sys_call(cmd)
shutil.rmtree(tmp)
empty = require_file_exists_and_test_empty(
os.path.join(outdir, base_genome + '.blastn.tsv'))
Expand Down

0 comments on commit 05f5ad2

Please sign in to comment.