Skip to content

cbalbin-bio/Submatch-BLAST-Parser

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 

History

5 Commits
 
 
 
 
 
 
 
 

Repository files navigation

If using outside of this module import the module.

from blastparser import BLASTParser

Instantiate BLASTParser

bp = BLASTParser("path/to/query/seq", "path/to/blast/db")

Alternatively, to filter blast results using taxids, pass a taxid or a list of taxids to exclude from blast search. This requires that you have get_species_taxids.sh in your path and you used a taxid map when creating your db.

get_species_taxids.sh is part of blast+

For example the call below will exlude all coronaviruses

bp = BLASTParser("path/to/query/seq", "path/to/blast/db", taxids=11118)

Retrieve the parsed hits.

This returns an instance of HitsDataContainer which is essentially just a list with some convience methods I wrote for filtering. The hits list contains instances of HitData.

hits = bp.gethits()

Filters out any hits not containing 5 consecutive matches somewhere in the hit.
hits= hits.filterbymatchlen(5)

You can chain filters together.

hits = hits.filterbymatchlen(5).filterbyident(.60).filterbycover(60)

Only printing first 3 hits for this example.

for hit in hits[:3]:
    # you can access the properties of the HitData for each
    print(f"Subject accession {hit.subject_accession}")
    print(" Query Match range(s):", "".join([f"{match_range[0]}-{match_range[1]}, " for match_range in hit.match_ranges]), "  Match len(s):", ", ".join([f"{match_len}" for match_len in hit.match_lengths]))
    print()
    print(f"Query  {hit.query_start:<5d} {hit.aln_query_seq} {hit.query_end:<5d}")
    print(f"            ", "".join(hit.cigar))
    print(f"Subj   {hit.subject_start:<5d} {hit.aln_subject_seq} {hit.subject_end:<5d}")
    print()


    # to handle subbhits (mutliple consecutive hits of >= 5 in a single blast hit)
    # remember hit.matches ranges is a lists of lists as structured as below
    # [[start_pos_int, end_pos_int], [start_pos_int, end_pos_int], ...]
    for index, (start, end) in enumerate(hit.match_ranges):

        # subtract 1 from start to go back to 0 based index (blast uses 1 based index) 
        if end - (start-1) >=5:
            #only print subhits greater than 5
            print(f"{start} {hit.submatch_seqs[index]} {end}")



    print()
    print()
    print()
    print()