enable gapped blast hits for porins, rev guided by @alawsin

Author: chrisgulvik

c-SSTAR

@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 
-__version__ = '1.1'
+__version__ = '1.2'
 
 import argparse
 import logging
@@ -30,12 +30,16 @@ def syscall(syscmd):
 			sys.exit('ERROR: failed syscall ' + syscmd)
 
 def translateSeq(nuclSeq):
+	'''use biopython to take in a nucleotide sequence (string) and return a protein sequence'''
 	proteinSeq = Seq(nuclSeq, generic_dna).translate(cds=False, to_stop=False, stop_symbol='*')
 	return proteinSeq
 
 def internalSTOPcodon(candidate):
 	'''counts number of internal stop codons; requires sequence from database to be in frame'''
-	nucSeq = candidate[10]
+	if '-' in candidate[10]:
+		nucSeq = candidate[10].replace('-', '')
+	else:
+		nucSeq = candidate[10]
 	if int(int(candidate[8])%3) == 1:
 		frameStart = int(0)
 	elif int(int(candidate[8])%3) == 0:
@@ -44,13 +48,12 @@ def internalSTOPcodon(candidate):
 		frameStart = int(2)
 	if frameStart > 0:
 		nucSeq = nucSeq[frameStart:]
-	frameEnd = int(int(candidate[9])%3)
-	if frameEnd > 0:
-		nucSeq = nucSeq[:-frameEnd]
 	if len(nucSeq)%3 == 0:
 		protein = translateSeq(nucSeq)
+	elif len(nucSeq)%3 > 0:
+		protein = translateSeq(nucSeq[:-(len(nucSeq)%3)])
 	else:
-		sys.exit('ERROR: incorrect nucleotide length (%s) after trimming' % len(nucSeq))
+		sys.exit('ERROR: incorrect nucleotide length ({}) after trimming'.format(len(nucSeq)))
 	numInternalSTOPcodons = len(re.findall(r'\*[ABCDEFGHIKLMNPQRSTVWYZ]', str(protein)))
 	return (protein, str(numInternalSTOPcodons))
 
@@ -74,22 +77,22 @@ def main(args=None):
 		baseGenome = os.path.splitext(os.path.basename(genome))[0]
 	if not os.path.exists(outdir): 
 		os.mkdir(outdir)
-	logging.basicConfig(filename='%s/c-SSTAR_%s.log' % (outdir, baseGenome), format='%(asctime)s: %(levelname)s: %(message)s', datefmt='%d-%m-%Y %I:%M:%S %p', level=logging.INFO)
-	logging.info('user: %s' % pwd.getpwuid(os.getuid()).pw_name)
-	logging.info('release: %s' % os.uname()[3])
-	logging.info('shell env: %s' % pwd.getpwuid(os.getuid()).pw_shell)
-	logging.info('cwd: %s' % pwd.getpwuid(os.getuid()).pw_dir)
-	logging.info('python version: %s' % sys.version)
+	logging.basicConfig(filename='{}/c-SSTAR_{}.log'.format(outdir, baseGenome), format='%(asctime)s: %(levelname)s: %(message)s', datefmt='%d-%m-%Y %I:%M:%S %p', level=logging.INFO)
+	logging.info('user: {}'.format(pwd.getpwuid(os.getuid()).pw_name))
+	logging.info('release: {}'.format(os.uname()[3]))
+	logging.info('shell env: {}'.format(pwd.getpwuid(os.getuid()).pw_shell))
+	logging.info('cwd: {}'.format(pwd.getpwuid(os.getuid()).pw_dir))
+	logging.info('python version: {}'.format(sys.version))
 	logging.info(subprocess.check_output('command -v blastn', shell=True).rstrip())
 	logging.info(subprocess.check_output('blastn -version | tail -n1', shell=True).rstrip())
 	database = args.database
 	similarity = args.similarity
-	syscall('makeblastdb -in ' + genome + ' -out ' + os.path.join(outdir, baseGenome) + ' -dbtype nucl')
-	syscall('blastn -task blastn -query {0} -db {out} -out {out}.blastn.tsv -evalue 1e-5 -max_target_seqs 1 -ungapped -perc_identity {1} -culling_limit 1 -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qlen sseq"'.format(database, similarity, out=os.path.join(outdir, baseGenome)))
-	syscall('rm ' + os.path.join(outdir, baseGenome) + '.nin')
-	syscall('rm ' + os.path.join(outdir, baseGenome) + '.nsq')
-	syscall('rm ' + os.path.join(outdir, baseGenome) + '.nhr')
-	with open(os.path.join(outdir, baseGenome) + '.blastn.tsv') as infile:
+	syscall('makeblastdb -in {} -out {} -dbtype nucl'.format(genome, os.path.join(outdir, baseGenome)))
+	syscall('blastn -task blastn -query {0} -db {out} -out {out}.blastn.tsv -evalue 1e-5 -max_target_seqs 1 -perc_identity {1} -culling_limit 1 -outfmt "6 qseqid sseqid pident length mismatch gaps qstart qend sstart send evalue bitscore qlen sseq"'.format(database, similarity, out=os.path.join(outdir, baseGenome)))
+	os.remove(os.path.join(outdir, baseGenome + '.nin'))
+	os.remove(os.path.join(outdir, baseGenome + '.nsq'))
+	os.remove(os.path.join(outdir, baseGenome + '.nhr'))
+	with open(os.path.join(outdir, baseGenome + '.blastn.tsv')) as infile:
 		best = ['-1','a','a','a',0,'-1','a',0,'a','a','a']
 		currentClusterNr = '-1'
 		topHits = []