# Path: kleborate/blastn.py
# def run_blastn(db, query, min_cov, min_ident):
#     build_blast_database_if_needed(db)
# 
#     cmd = 'blastn -task blastn -db {} -query {}'.format(db, query)
#     cmd += " -outfmt '6 sacc pident slen length bitscore qseq sstrand sstart send" \
#            " qacc qstart qend qframe'"
#     cmd += ' -dust no -evalue 1E-20 -word_size 32 -max_target_seqs 10000'
#     cmd += ' -perc_identity {}'.format(min_ident)
# 
#     # TODO: switch this over to subprocess
#     blast_hits = []
#     f = os.popen(cmd)
#     for line in f:
#         blast_hits.append(BlastHit(line))
#     f.close()
# 
#     # Toss out low identity and low coverage hits.
#     if min_ident is not None:
#         blast_hits = [h for h in blast_hits if h.pcid * 100 >= min_ident]
#     if min_cov is not None:
#         blast_hits = [h for h in blast_hits if h.ref_cov * 100 >= min_cov]
# 
#     # Clean up redundant hits so we only have one hit for each part of the genome.
#     blast_hits = cull_redundant_hits(blast_hits)
# 
#     # Get rid of any hits that start with 'delete_'. This is used in cases where there can be
#     # confounding homologous hits. E.g. searching for rmpA2 alleles but not wanting to be misled by
#     # homologous rmpA alleles.
#     blast_hits = [h for h in blast_hits if not h.gene_id.startswith('delete_')]
# 
#     return blast_hits
# 
# Path: kleborate/truncation.py
# def truncation_check(hit, cov_threshold=90.0):
#     """
#     Checks to see if the gene is truncated at the amino acid level.
#     """
#     nucl_seq, ref_start, _ = hit.get_seq_start_end_pos_strand()
# 
#     # The hit must start at the first base of the gene. If not, the gene is considered 0%.
#     if ref_start != 1:
#         return '-0%', 0.0, ''
# 
#     # If there are any ambiguous bases in the sequence, then they will break translation, probably
#     # resulting in truncation call.
#     ambiguous_bases = set(b for b in nucl_seq) - {'A', 'C', 'G', 'T'}
#     for b in ambiguous_bases:
#         nucl_seq = nucl_seq.split(b)[0]
# 
#     # BioPython doesn't like it if the sequence isn't a multiple of 3.
#     nucl_seq = nucl_seq[:len(nucl_seq) // 3 * 3]
# 
#     # The assumption is that the reference allele is a full CDS with a stop codon at the end. This
#     # isn't always true (the reference sequence is sometimes broken) but will serve to make our
#     # denominator for coverage.
#     ref_aa_length = (hit.ref_length - 3) // 3
# 
#     coding_dna = Seq(nucl_seq)
#     translation = str(coding_dna.translate(table='Bacterial', to_stop=True))
# 
#     coverage = 100.0 * len(translation) / ref_aa_length
#     if coverage >= cov_threshold:
#         return '', coverage, translation
#     else:
#         return '-{:.0f}%'.format(coverage), coverage, translation


from .blastn import run_blastn
from .truncation import truncation_check
If an annotation column is provided (such as clonal complex) in the final column of the profiles
file, this annotation will be reported in column 2 of the output table

NOTE there is a bug with the culling_limit parameter in older versions of BLAST+. This code has
been tested with BLAST+2.2.30. It does not work with BLAST2.2.25. Not sure about other versions.

Copyright 2020 Kat Holt
Copyright 2020 Ryan Wick (rrwick@gmail.com)
https://github.com/katholt/Kleborate/

This file is part of Kleborate. Kleborate is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by the Free Software Foundation,
either version 3 of the License, or (at your option) any later version. Kleborate is distributed in
the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
details. You should have received a copy of the GNU General Public License along with Kleborate. If
not, see <http://www.gnu.org/licenses/>.
"""



def rmpa2_blast(seqs, assemblies, min_cov, min_ident):
    for contigs in assemblies:
        rmpa2_calls = []
        hits = run_blastn(seqs, contigs, min_cov, min_ident)
        for hit in hits:
            if hit.alignment_length > hit.ref_length / 2:
                gene_id = hit.gene_id
                if hit.pcid < 100.00 or hit.alignment_length < hit.ref_length:
                    gene_id += '*'

# Path: kleborate/blastn.py
# def run_blastn(db, query, min_cov, min_ident):
#     build_blast_database_if_needed(db)
# 
#     cmd = 'blastn -task blastn -db {} -query {}'.format(db, query)
#     cmd += " -outfmt '6 sacc pident slen length bitscore qseq sstrand sstart send" \
#            " qacc qstart qend qframe'"
#     cmd += ' -dust no -evalue 1E-20 -word_size 32 -max_target_seqs 10000'
#     cmd += ' -perc_identity {}'.format(min_ident)
# 
#     # TODO: switch this over to subprocess
#     blast_hits = []
#     f = os.popen(cmd)
#     for line in f:
#         blast_hits.append(BlastHit(line))
#     f.close()
# 
#     # Toss out low identity and low coverage hits.
#     if min_ident is not None:
#         blast_hits = [h for h in blast_hits if h.pcid * 100 >= min_ident]
#     if min_cov is not None:
#         blast_hits = [h for h in blast_hits if h.ref_cov * 100 >= min_cov]
# 
#     # Clean up redundant hits so we only have one hit for each part of the genome.
#     blast_hits = cull_redundant_hits(blast_hits)
# 
#     # Get rid of any hits that start with 'delete_'. This is used in cases where there can be
#     # confounding homologous hits. E.g. searching for rmpA2 alleles but not wanting to be misled by
#     # homologous rmpA alleles.
#     blast_hits = [h for h in blast_hits if not h.gene_id.startswith('delete_')]
# 
#     return blast_hits
# 
# Path: kleborate/truncation.py
# def truncation_check(hit, cov_threshold=90.0):
#     """
#     Checks to see if the gene is truncated at the amino acid level.
#     """
#     nucl_seq, ref_start, _ = hit.get_seq_start_end_pos_strand()
# 
#     # The hit must start at the first base of the gene. If not, the gene is considered 0%.
#     if ref_start != 1:
#         return '-0%', 0.0, ''
# 
#     # If there are any ambiguous bases in the sequence, then they will break translation, probably
#     # resulting in truncation call.
#     ambiguous_bases = set(b for b in nucl_seq) - {'A', 'C', 'G', 'T'}
#     for b in ambiguous_bases:
#         nucl_seq = nucl_seq.split(b)[0]
# 
#     # BioPython doesn't like it if the sequence isn't a multiple of 3.
#     nucl_seq = nucl_seq[:len(nucl_seq) // 3 * 3]
# 
#     # The assumption is that the reference allele is a full CDS with a stop codon at the end. This
#     # isn't always true (the reference sequence is sometimes broken) but will serve to make our
#     # denominator for coverage.
#     ref_aa_length = (hit.ref_length - 3) // 3
# 
#     coding_dna = Seq(nucl_seq)
#     translation = str(coding_dna.translate(table='Bacterial', to_stop=True))
# 
#     coverage = 100.0 * len(translation) / ref_aa_length
#     if coverage >= cov_threshold:
#         return '', coverage, translation
#     else:
#         return '-{:.0f}%'.format(coverage), coverage, translation


# Path: kleborate/rmpA.py
from .blastn import run_blastn
from .truncation import truncation_check

If an annotation column is provided (such as clonal complex) in the final column of the profiles
file, this annotation will be reported in column 2 of the output table

NOTE there is a bug with the culling_limit parameter in older versions of BLAST+. This code has
been tested with BLAST+2.2.30. It does not work with BLAST2.2.25. Not sure about other versions.

Copyright 2020 Kat Holt
Copyright 2020 Ryan Wick (rrwick@gmail.com)
https://github.com/katholt/Kleborate/

This file is part of Kleborate. Kleborate is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by the Free Software Foundation,
either version 3 of the License, or (at your option) any later version. Kleborate is distributed in
the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
details. You should have received a copy of the GNU General Public License along with Kleborate. If
not, see <http://www.gnu.org/licenses/>.
"""



def rmpa2_blast(seqs, assemblies, min_cov, min_ident):
    for contigs in assemblies:
        rmpa2_calls = []
        hits = run_blastn(seqs, contigs, min_cov, min_ident)
        for hit in hits:
            if hit.alignment_length > hit.ref_length / 2:
                gene_id = hit.gene_id
                if hit.pcid < 100.00 or hit.alignment_length < hit.ref_length:
                    gene_id += '*'

                gene_id += truncation_check(hit)[0]