# How to Extract Sequences from BLAST Hits in Python By the end of this guide you will be able to run a BLAST search from the command line, parse the results in Python, filter the hits by identity, coverage, and e-value, and write the matching sequences from your source FASTA files into a clean output FASTA. You will also be able to pull out just the aligned region of a hit, handle reverse-strand matches correctly, and verify that the extraction actually worked. Everything here runs on a standard laptop with no cluster access required. ## What You Need Before Starting ![Five step flow from running BLAST in tabular format through parsing, indexing, slicing and writing a verified FASTA file.](/article-images/illustrations/extract-sequences-from-blast-hits-python-blast-pipeline.webp) Install and have on hand the following: - **Python 3.8 or newer.** Check with `python3, version`. - **Biopython.** Install with `pip install biopython`. Biopython is the standard library for python for bioinformatics work and handles FASTA reading, writing, and reverse complementation. - **NCBI BLAST+ command line tools.** These provide `blastn`, `blastp`, `blastx`, `tblastn`, and `makeblastdb`. On macOS use `brew install blast`, on Ubuntu use `sudo apt install ncbi-blast+`, or download the binaries from the NCBI BLAST+ download page. - **A query FASTA file** with the sequences you searched with, for example `queries.fasta`. - **A subject FASTA file** with the sequences you searched against, for example `subjects.fasta`. If you searched an NCBI online database instead of a local one, download the matching database sequences or use the NCBI web interface to fetch them by accession. - **A BLAST results file.** If you do not have one yet, this guide shows you how to produce it in the format that is easiest to parse. A quick note on why a scripted approach beats manual work. As the number of sequences grows, manual inspection of BLAST results becomes tedious, and researchers typically need tools that parse and filter results programmatically [1]. Large BLAST searches are also slow to analyze by hand, which is exactly why dedicated parsing tools keep getting written for this problem [2]. A short Python script is often the most flexible option because you control every filter and every output field. ## Step 1: Run BLAST in a Parse-Friendly Format BLAST can emit many output formats. The single most useful choice for scripting is tabular output, requested with `-outfmt 6`. Tabular output places one hit per line with fields separated by tabs, which means one line of Python splits it into columns. First, if you are searching your own FASTA file rather than a prebuilt database, turn the subject file into a BLAST database: ```bash makeblastdb -in subjects.fasta -dbtype nucl -out mydb ``` Use `-dbtype nucl` for nucleotide sequences and `-dbtype prot` for proteins. Then run the search: ```bash blastn -query queries.fasta -db mydb \ -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qlen slen" \ -evalue 1e-5 -max_target_seqs 5 \ -out hits.tsv ``` Here is what each flag does: - `-query` names the FASTA file with your sequences of interest. - `-db` points at the database you built with `makeblastdb`. Use the same base name you gave `-out`. - `-outfmt "6 ..."` requests tabular output and lets you list the exact columns, in the exact order you want them. Listing columns explicitly protects you from format changes and makes the parser self-documenting. - `-evalue 1e-5` keeps only hits with an expect value of 0.00001 or lower. This is a conventional first-pass cutoff, not a universal rule. - `-max_target_seqs 5` limits the number of aligned sequences kept per query. This flag behaves in surprising ways, covered later in the mistakes section, because it can affect the search process itself rather than acting purely as a post-search display filter [2]. - `-out` writes results to a file instead of standard output. If you prefer the default tabular columns, plain `-outfmt 6` gives you these twelve fields in this order: | Column number | Field name | Meaning | Python index | |--|--|--|--| | 1 | qseqid | Query sequence identifier | 0 | | 2 | sseqid | Subject sequence identifier | 1 | | 3 | pident | Percent identical positions | 2 | | 4 | length | Alignment length | 3 | | 5 | mismatch | Number of mismatches | 4 | | 6 | gapopen | Number of gap openings | 5 | | 7 | qstart | Query alignment start (1-based) | 6 | | 8 | qend | Query alignment end (1-based) | 7 | | 9 | sstart | Subject alignment start (1-based) | 8 | | 10 | send | Subject alignment end (1-based) | 9 | | 11 | evalue | Expect value | 10 | | 12 | bitscore | Bit score | 11 | The Python index column is the position of that field after splitting a line, and it is one less than the human-readable column number because Python counts from zero. This one-offset difference is the single most common source of bugs when you extract sequences from BLAST hits in Python, so internalize it now. Peek at your results to confirm the file looks right: ```bash head -n 5 hits.tsv ``` Expected output looks like this (your values will differ): ``` contig_101 scaffold_22 98.403 312 5 1 1 312 45021 45332 2.1e-152 538 contig_101 scaffold_87 92.949 312 21 2 1 310 120455 120766 4.7e-118 424 contig_214 scaffold_3 99.678 154 0 0 12 165 8092 8245 3.3e-76 277 ``` Every line is one query-subject pairwise hit. A single query can produce many lines. ## Step 2: Parse the Tabular File in Plain Python You do not strictly need Biopython to read the hits file. Plain Python handles it cleanly: ```python def parse_hits(path): """Read a tabular BLAST file and return a list of hit dictionaries.""" fields = ["qseqid", "sseqid", "pident", "length", "mismatch", "gapopen", "qstart", "qend", "sstart", "send", "evalue", "bitscore", "qlen", "slen"] hits = [] with open(path) as handle: for line in handle: line = line.rstrip("\n") if not line: continue values = line.split("\t") hit = dict(zip(fields, values)) # Convert the fields you will compute with for key in ("pident", "evalue", "bitscore"): hit[key] = float(hit[key]) for key in ("length", "mismatch", "gapopen", "qstart", "qend", "sstart", "send", "qlen", "slen"): hit[key] = int(hit[key]) hits.append(hit) return hits hits = parse_hits("hits.tsv") print(len(hits), "hits parsed") print(hits) ``` Expected output: ``` 3 hits parsed {'qseqid': 'contig_101', 'sseqid': 'scaffold_22', 'pident': 98.403, 'length': 312, 'mismatch': 5, 'gapopen': 1, 'qstart': 1, 'qend': 312, 'sstart': 45021, 'send': 45332, 'evalue': 2.1e-152, 'bitscore': 538.0, 'qlen': 312, 'slen': 61234} ``` The dictionary approach beats positional indexing because if you later add or reorder columns in `-outfmt`, your code keeps working as long as you update the `fields` list to match. This is one of the highest-value tips for anyone learning how to extract sequences from BLAST hits in Python: never rely on bare numeric indices like `line.split("\t")` scattered through a script. One column insertion silently shifts every index and produces wrong coordinates with no error message. ## Step 3: Load the Subject Sequences with Biopython To extract sequences from BLAST hits you need the actual sequence text, and that lives in your source FASTA files, not in the BLAST output. Load the subject FASTA into a dictionary-like object: ```python from Bio import SeqIO # index() reads headers eagerly but loads sequence text lazily, # which matters when the subject file is large subject_db = SeqIO.index("subjects.fasta", "fasta") print(len(subject_db), "subjects indexed") record = subject_db["scaffold_22"] print(record.id, "length", len(record.seq)) ``` Expected output: ``` 412 subjects indexed scaffold_22 length 61234 ``` Two functions are available and the choice matters: - `SeqIO.index(filename, "fasta")` builds a lightweight on-disk or in-memory index and only reads a record's sequence when you request it. Use this for genome-sized or otherwise large FASTA files. - `SeqIO.to_dict(SeqIO.parse(filename, "fasta"))` loads every sequence into memory at once. This is fine for a few thousand short sequences and simpler to reason about, but it will exhaust RAM on a full genome. The same applies to the query file if you want to extract query sequences too: ```python query_db = SeqIO.index("queries.fasta", "fasta") ``` One caution: the identifiers in the BLAST output are the first whitespace-delimited token of each FASTA header. If your headers look like `>scaffold_22 whole genome shotgun`, the BLAST hit id is `scaffold_22` and `subject_db["scaffold_22 whole genome shotgun"]` will raise a `KeyError`. `SeqIO.index` has an optional `key_function` argument that lets you transform headers during indexing if your files use full descriptions as identifiers. ## Step 4: Filter the Hits Extracting everything is rarely the goal. Filter first, extract second. A typical filter keeps hits that are long, similar, and statistically significant: ```python def passes_filter(hit, min_pident=90.0, min_aln_frac=0.5, max_evalue=1e-20): """A hit must be similar, cover a reasonable fraction of the query, and be statistically significant.""" if hit["pident"] < min_pident: return False if hit["evalue"] > max_evalue: return False # Alignment length as a fraction of the query length if hit["length"] < min_aln_frac * hit["qlen"]: return False return True kept = [h for h in hits if passes_filter(h)] print(f"{len(kept)} of {len(hits)} hits passed the filter") ``` Expected output: ``` 2 of 3 hits passed the filter ``` Adjust the thresholds to your question. For primers or species identification you might demand 98 to 100 percent identity over nearly the full query. For remote homology discovery you would relax all three thresholds and lean more on bit score and domain evidence. Percent identity is the fraction of aligned positions that match exactly, the e-value estimates how many hits with a score at least this good you would expect to see by chance in a database of this size, and alignment length counts aligned columns including gaps. The filtering-and-retrieval pattern shown here is the same pattern behind server-based BLAST result management systems, which let users interrogate stored results with logical operations instead of rerunning searches [1]. Doing it locally in a script gives you the same flexibility with no server dependency. ## Step 5: Extract Full Sequences for the Hits The most common request is simple: give me the complete sequence of every subject that matched. Keep the order stable, deduplicate subjects that were hit by multiple queries, and write a FASTA: ```python from Bio import SeqIO from Bio.SeqRecord import SeqRecord subject_db = SeqIO.index("subjects.fasta", "fasta") hits = parse_hits("hits.tsv") kept = [h for h in hits if passes_filter(h)] seen = set() records = [] for hit in kept: sid = hit["sseqid"] if sid in seen: continue seen.add(sid) record = subject_db[sid] # Annotate the new header so you keep the evidence trail records.append(SeqRecord( record.seq, id=record.id, description=(f"query={hit['qseqid']} pident={hit['pident']:.1f} " f"evalue={hit['evalue']:.1e}") )) SeqIO.write(records, "hit_subjects.fasta", "fasta") print(f"Wrote {len(records)} sequences to hit_subjects.fasta") ``` Expected output: ``` Wrote 2 sequences to hit_subjects.fasta ``` The output file looks like this: ``` >scaffold_22 query=contig_101 pident=98.4 evalue=2.1e-152 GATCGATCGGTAC... ``` Carrying the query id and the key statistics into the FASTA description is one of the most practical tips in this biopython tutorial. Months later, when you find `scaffold_22` in a downstream alignment, the header tells you which query pulled it out and how strong the evidence was. Without that annotation, extracted sequence files become anonymous piles that no one trusts. ## Step 6: Extract the Aligned Region Instead of the Whole Sequence ![Decision path checking coordinate order and search task to pick a forward slice, a reverse complement slice, or translation.](/article-images/illustrations/extract-sequences-from-blast-hits-python-strand-slicing.webp) Often you do not want the entire 60 kilobase scaffold. You want the exact stretch that aligned. This is where the coordinate columns earn their keep, and where the 1-based versus 0-based conversion matters: ```python def extract_aligned_region(subject_db, hit): """Return the subject sequence segment covered by the alignment, reverse complemented when the match is on the minus strand.""" record = subject_db[hit["sseqid"]] start, end = hit["sstart"], hit["send"] # BLAST coordinates are 1-based and inclusive. # Python slices are 0-based and half-open. if start <= end: subseq = record.seq[start - 1:end] else: # Minus strand hit: sstart is greater than send subseq = record.seq[end - 1:start].reverse_complement() return subseq sub_db = SeqIO.index("subjects.fasta", "fasta") for hit in kept: subseq = extract_aligned_region(sub_db, hit) print(hit["sseqid"], hit["sstart"], hit["send"], len(subseq), subseq[:30]) ``` Expected output: ``` scaffold_22 45021 45332 312 GATCGATCGGTACATTGCACTGGATCGATC scaffold_3 8092 8245 154 TTGACCATGGCTAGCTTACGGATCCATGA ``` Walk through the arithmetic once so it never bites you again. BLAST reports positions counting from 1, and both the start and end positions are included in the alignment. Python slices count from 0, and the end position is excluded. For an alignment from position 45021 to 45332, the correct slice is `seq[45020:45332]`, which has length 312, matching the alignment length column exactly. The general rule is `seq[start - 1:end]`. The strand logic works because BLAST reports coordinates on the subject's plus strand orientation. When the subject match is on the minus strand, `sstart` is the higher coordinate and `send` the lower one, so the plus-strand interval still runs from `send` to `sstart`. Slicing that interval and applying `reverse_complement()` recovers the sequence exactly as the query saw it. If you skip the reverse complementation, your extracted fragments will fail to align to the query in downstream steps, which is a classic symptom worth remembering. For protein searches such as `blastp`, strand is irrelevant and the simple forward slice is all you need. For translated searches such as `blastx` and `tblastn`, the query or subject coordinates are in amino acid space while the FASTA is in nucleotide space, so you must translate the extracted nucleotide region using the correct genetic code table rather than slicing alone. If you are just starting out, do your extraction exercises with `blastn` against a nucleotide database, where coordinates and sequences live in the same units. ## Step 7: Extract Matching Query Sequences Too The identical pattern retrieves the query-side sequence, useful when you want both partners of each hit for later alignment: ```python query_db = SeqIO.index("queries.fasta", "fasta") for hit in kept: qstart, qend = hit["qstart"], hit["qend"] record = query_db[hit["qseqid"]] if qstart <= qend: qsubseq = record.seq[qstart - 1:qend] else: qsubseq = record.seq[qend - 1:qstart].reverse_complement() print(hit["qseqid"], len(qsubseq)) ``` For nucleotide-to-nucleotide searches the query strand logic mirrors the subject logic. In `blastn` output a query on the minus strand still shows `qstart > qend`. You can also request the `sstrand` column explicitly in `-outfmt`, where `blastn` writes `plus` or `minus`. Relying on coordinate ordering is more portable because some task combinations do not populate `sstrand` the way you might expect. ## Step 8: Write Everything to One Clean FASTA Combine full-subject extraction and aligned-region extraction into one script you can reuse: ```python """extract_blast_hits.py Extract sequences for BLAST hits from a tabular results file. Usage: python extract_blast_hits.py hits.tsv subjects.fasta full python extract_blast_hits.py hits.tsv subjects.fasta aligned """ import sys from Bio import SeqIO from Bio.SeqRecord import SeqRecord def parse_hits(path): fields = ["qseqid", "sseqid", "pident", "length", "mismatch", "gapopen", "qstart", "qend", "sstart", "send", "evalue", "bitscore", "qlen", "slen"] hits = [] with open(path) as handle: for line in handle: if not line.strip(): continue values = line.rstrip("\n").split("\t") hit = dict(zip(fields, values)) hit["pident"] = float(hit["pident"]) hit["evalue"] = float(hit["evalue"]) for key in ("length", "qstart", "qend", "sstart", "send", "qlen", "slen"): hit[key] = int(hit[key]) hits.append(hit) return hits def passes_filter(hit, min_pident=90.0, max_evalue=1e-20): return (hit["pident"] >= min_pident and hit["evalue"] <= max_evalue and hit["length"] >= 0.5 * hit["qlen"]) def main(): hits_path, fasta_path, mode = sys.argv[2], sys.argv[1], sys.argv db = SeqIO.index(fasta_path, "fasta") hits = [h for h in parse_hits(hits_path) if passes_filter(h)] seen = set() records = [] for hit in hits: sid = hit["sseqid"] if sid in seen: continue seen.add(sid) source = db[sid] if mode == "full": seq = source.seq else: s, e = hit["sstart"], hit["send"] if s <= e: seq = source.seq[s - 1:e] else: seq = source.seq[e - 1:s].reverse_complement() records.append(SeqRecord(seq, id=sid, description=f"query={hit['qseqid']} pident={hit['pident']:.1f}")) out_name = f"hit_{mode}.fasta" SeqIO.write(records, out_name, "fasta") print(f"Wrote {len(records)} sequences to {out_name}") if __name__ == "__main__": main() ``` Run it: ```bash python extract_blast_hits.py hits.tsv subjects.fasta full python extract_blast_hits.py hits.tsv subjects.fasta aligned ``` Expected output: ``` Wrote 2 sequences to hit_full.fasta Wrote 2 sequences to hit_aligned.fasta ``` This script is a solid scaffold. In real projects you will extend it with additional filters, a two-input mode for query and subject extraction, or per-query hit limits. The core loop, coordinate math, and deduplication rarely change. ## An Alternative: Parsing BLAST XML Output Tabular output is compact, but it loses some information such as the full alignment strings and per-HSP detail. XML output, requested with `-outfmt 5`, preserves everything, and Biopython's `Bio.Blast.NCBIXML` parser reads it: ```python from Bio.Blast import NCBIXML with open("hits.xml") as handle: for record in NCBIXML.parse(handle): for alignment in record.alignments: for hsp in alignment.hsps: print(record.query, alignment.hit_id, hsp.expect, hsp.sbjct_start, hsp.sbjct_end) ``` Expected output: ``` contig_101 scaffold_22 2.1e-152 45021 45332 ``` One documented pitfall with XML output: the coordinate fields in NCBIXML are 0-based, while tabular output coordinates are 1-based. If you parse both formats in one project, write the slicing rule down in one helper function and use it everywhere rather than scattering slice arithmetic across your code. Biopython also offers `Bio.SearchIO`, a unified parsing interface that reads tabular, XML, and many other alignment formats with one API: ```python from Bio import SearchIO # Works with default outfmt 6 (the standard 11 or 12 columns) for qresult in SearchIO.parse("hits_default.tsv", "blast-tab"): for hit in qresult: for hsp in hit.hsps: print(qresult.id, hit.id, hsp.hit_start, hsp.hit_end, hsp.evalue) ``` `SearchIO` normalizes coordinates to 0-based half-open intervals internally, so `hit_start` and `hit_end` can be used directly as Python slice bounds without adjustment. If your tabular file uses a custom column order, `SearchIO.parse` accepts a `fields` keyword listing your column names in order. The tradeoff is that `SearchIO` has a learning curve, and for the plain extraction task in this guide the simple dictionary parser above is easier to audit. Performance is also a real consideration at scale. BLAST-QC was developed specifically because parsing large BLAST result sets was tedious with existing options, and its authors reported better runtimes than parsers built on Biopython and BioPerl for their workloads [2]. If you routinely process result files with tens of millions of lines, benchmark your parser and consider streaming approaches that never hold the whole file in memory, such as iterating over lines and writing output as you go. ## Verifying That the Extraction Worked Never trust a script that writes sequence files without at least three checks. **Count and deduplication check.** The number of FASTA records written should equal the number of unique subject ids among filtered hits: ```bash grep -c ">" hit_full.fasta ``` ```python print(len({h["sseqid"] for h in kept}), "unique subjects", len(kept), "filtered hits") ``` These two numbers match only when your deduplication is correct. **Length check.** For aligned-region extraction, the extracted length should equal `send - sstart + 1` for plus-strand hits, and should match the `length` column from BLAST unless the alignment contains gaps: ```python for hit in kept: span = abs(hit["send"] - hit["sstart"]) + 1 print(hit["sseqid"], "span", span, "aln_length", hit["length"]) ``` When the alignment contains gaps, the span along the subject can exceed the alignment length column, because gaps consume query positions without consuming subject positions. Seeing a span slightly larger than the alignment length is normal. Seeing a span wildly different from the alignment length means your slicing is wrong. **Round-trip check.** Take one extracted sequence, search it back against the same database, and confirm you recover the original hit coordinates: ```bash blastn -query hit_aligned.fasta -db mydb \ -outfmt "6 qseqid sseqid pident length sstart send evalue" \ -max_target_seqs 1 ``` Expected output: ``` scaffold_22 scaffold_22 100.000 312 45021 45332 2.1e-152 ``` A perfect self-hit at the original coordinates confirms the extraction, the strand handling, and the coordinate math all at once. This round-trip trick catches reverse-complement errors and off-by-one slicing faster than any amount of staring at code. ## Common Mistakes and Limitations These are the errors that actually show up in real analysis pipelines, with their symptoms and fixes. **Off-by-one slicing.** Using `seq[start:end]` on 1-based BLAST coordinates drops the first base of every extraction. The fix is always `seq[start - 1:end]`. The symptom is subtle: extracted regions are one base short and self-hits report an alignment one shorter than expected. **Ignoring strand direction.** Treating a minus-strand hit as a plus-strand slice gives you the reverse complement of what you wanted. Symptom: extracted fragments contain plausible sequence but downstream BLAST or alignment finds no match to the query. Fix: check whether `sstart > send` and reverse complement when it is. **Mixing up query and subject coordinates.** With `tblastn` and `blastx`, one side is in amino acid coordinates and the other in nucleotide coordinates. Slicing a [nucleotide sequence](/knowledge/molecular-biology/nucleotide-sequence) with protein coordinates produces fragments about three times too short. Fix: know which task you ran, and convert or translate accordingly. **Hardcoding column indices.** `values` breaks silently the day you add a column to `-outfmt`. Fix: parse into named fields as shown in Step 2, and keep the `fields` list synchronized with the `-outfmt` string. **Trusting `-max_target_seqs` as a strict top-N guarantee.** This flag can influence the heuristic search itself, not merely the number of displayed results, which means the hits you keep may depend on the flag's value in ways that are hard to predict [2]. If you need a reliable top-N per query, retrieve a generous number of hits and apply your own ranking and truncation in Python based on bit score or e-value. **Trusting hit order in multi-query files.** Lines are grouped by query, and within a query group, ordering by e-value is generally expected for a single-query search, but relying on input line order as a ranking across the whole file is fragile. Sort explicitly in Python: ```python by_query = {} for h in kept: by_query.setdefault(h["qseqid"], []).append(h) for qid, group in by_query.items(): group.sort(key=lambda h: (-h["bitscore"], h["evalue"])) ``` **Assuming e-value alone means biological relevance.** A tiny e-value on a 15-base alignment is weak evidence. Short alignments can be statistically significant yet biologically meaningless, especially in repeat-rich sequences. Always require a minimum alignment length or coverage fraction alongside the e-value. **Masking low-complexity regions without knowing it.** BLAST applies default filtering of low-complexity sequence. Alignments may exclude positions your query contains, so an alignment that does not span the full query is not automatically a truncated match. Check whether filtered or low-complexity sequence explains the gap before assuming a truncated gene. **Loading a full genome into memory with `SeqIO.to_dict`.** Symptom: your script grinds to a halt or the process is killed. Fix: use `SeqIO.index`, which loads sequence text on demand. **Extracting from an online search without local data.** If you ran BLAST on the NCBI web interface, the results contain identifiers and coordinates but not sequence text for subjects. You still need the subject FASTA, either downloaded by accession or via a local database, before any Python extraction can work. **Forgetting that a query can hit the same subject many times.** Repetitive regions produce multiple HSPs per query-subject pair, each on its own line. A naive full-sequence deduplication handles this gracefully, but naive aligned-region extraction will write overlapping fragments. Decide whether you want the single best HSP per pair, or all HSPs merged. Finally, a limitation that no script fixes: coordinate-based extraction reproduces exactly what BLAST aligned, and BLAST alignments are heuristics. A locally alignable region outside the reported HSPs will be missed. When in doubt, extract with generous flanks, for example 50 bases on each side, and re-align. ## Sensible Next Steps Once extraction works, the extracted FASTA files become inputs to the next stage of a standard workflow: - **[Multiple sequence alignment](/blog/guides/multiple-sequence-alignment)** of extracted homologs with MAFFT or MUSCLE, followed by phylogenetic analysis. - **Reciprocal best hit analysis** for orthology. Extract the best hit of A against B, and the best hit of B against A, and keep pairs that agree. Filtering results by logical operations such as reciprocal matches, shared targets, or multi-region matches is exactly the kind of question dedicated BLAST result systems were built to answer [1]. The same logic is a short Python function once both directions of hits are parsed. - **Primer or probe design** from high-identity extracted regions, after checking the flanks. - **Annotation pipelines** where extracted candidate regions get scanned for open reading frames or compared against protein domains. ## Frequently Asked Questions ### Can I extract sequences from BLAST hits without the original FASTA files? No. The BLAST output contains identifiers, coordinates, and scores, but not the full subject sequence text. You need the subject FASTA, or a way to fetch those accessions from a database, before extraction is possible. ### Why are some of my extracted sequences reversed? Because the hit was on the minus strand. When `sstart` is greater than `send` in tabular output, the match lies on the reverse strand, so you must slice from the lower coordinate to the higher one and apply `reverse_complement()`. ### What BLAST output format is best for Python parsing? Tabular output with `-outfmt 6` and an explicit column list is the best default, because it is small, one hit per line, and trivially parsed. Use `-outfmt 5` XML when you need alignment strings and full HSP detail that tabular output omits. ### How do I extract only the aligned part of a hit rather than the whole sequence? Slice the subject sequence with `seq[sstart - 1:send]` for plus-strand hits, or `seq[send - 1:sstart].reverse_complement()` when `sstart` exceeds `send`. Remember BLAST coordinates are 1-based and inclusive while Python slices are 0-based and exclude the end position. ### Do I need Biopython or can I use plain Python? Plain Python is fully sufficient for parsing the tabular hit file and filtering rows. You need Biopython for convenient FASTA reading and writing, reverse complementation, XML parsing, and the unified SearchIO interface, which makes it worth installing for any real project. ### Why do I get fewer or different hits than my -max_target_seqs setting suggests? The flag interacts with the search heuristic itself rather than acting only as a display cutoff, so the retained hits can vary in unexpected ways. Retrieve a larger number and apply your own top-N ranking by bit score in Python for reproducible results. ### How do I avoid writing the same subject sequence many times? Track subject identifiers in a Python `set` as you loop over hits, and skip any identifier you have already written. Keep the first, best-scoring occurrence, or merge coordinates if you need a consensus region. ### How do I keep only reciprocal best hits between two sequence sets? Run BLAST in both directions, parse both result files with the same parser, select the top hit per query in each direction, and keep pairs where each sequence is the other's best hit. This reciprocal best hit criterion is a standard way to find candidate orthologs and is easy to implement once both hit tables are in dictionaries. ## Sources 1. [The PARIGA server for real time filtering and analysis of reciprocal BLAST results.](https://pubmed.ncbi.nlm.nih.gov/23667459/) 2. [BLAST-QC: automated analysis of BLAST results.](https://pubmed.ncbi.nlm.nih.gov/33902722/)