|
| 1 | +#!/usr/bin/env python |
| 2 | + |
| 3 | +"""Provide a command line tool to filter blast results.""" |
| 4 | + |
| 5 | +import argparse |
| 6 | +import logging |
| 7 | +import sys |
| 8 | +from pathlib import Path |
| 9 | + |
| 10 | +import pandas as pd |
| 11 | +from Bio import SeqIO |
| 12 | + |
| 13 | +logger = logging.getLogger() |
| 14 | + |
| 15 | + |
| 16 | +def parse_args(argv=None): |
| 17 | + """Define and immediately parse command line arguments.""" |
| 18 | + parser = argparse.ArgumentParser( |
| 19 | + description="Provide a command line tool to filter blast results.", |
| 20 | + epilog="Example: python blast_filter.py in.clstr prefix", |
| 21 | + ) |
| 22 | + |
| 23 | + parser.add_argument( |
| 24 | + "-i", |
| 25 | + "--mash", |
| 26 | + metavar="MASH FILE", |
| 27 | + type=Path, |
| 28 | + help="Mash screen result file.", |
| 29 | + ) |
| 30 | + |
| 31 | + parser.add_argument( |
| 32 | + "-r", |
| 33 | + "--references", |
| 34 | + metavar="REFERENCE FILE", |
| 35 | + type=Path, |
| 36 | + help="Contig sequence file was screened against", |
| 37 | + ) |
| 38 | + |
| 39 | + parser.add_argument( |
| 40 | + "-p", |
| 41 | + "--prefix", |
| 42 | + metavar="PREFIX", |
| 43 | + type=str, |
| 44 | + help="Output file prefix", |
| 45 | + ) |
| 46 | + |
| 47 | + parser.add_argument( |
| 48 | + "-l", |
| 49 | + "--log-level", |
| 50 | + help="The desired log level (default WARNING).", |
| 51 | + choices=("CRITICAL", "ERROR", "WARNING", "INFO", "DEBUG"), |
| 52 | + default="INFO", |
| 53 | + ) |
| 54 | + return parser.parse_args(argv) |
| 55 | + |
| 56 | +def to_dict_remove_dups(sequences): |
| 57 | + return {record.id: record for record in sequences} |
| 58 | + |
| 59 | + |
| 60 | +def extract_hits(df, references, prefix): |
| 61 | + """ |
| 62 | + Extracts contigs hits from a DataFrame and writes them to a FASTA file. |
| 63 | +
|
| 64 | + Args: |
| 65 | + df (pandas.DataFrame): DataFrame containing the hits information. |
| 66 | + contigs (str): Path to the contigs file. |
| 67 | + references (str): Path to the references file in FASTA format. |
| 68 | + prefix (str): Prefix for the output file. |
| 69 | +
|
| 70 | + Returns: |
| 71 | + None |
| 72 | + """ |
| 73 | + try: |
| 74 | + ref_records = SeqIO.to_dict(SeqIO.parse(references, "fasta")) |
| 75 | + except ValueError as e: |
| 76 | + logger.warning( |
| 77 | + "Indexing the reference pool file causes an error: %s \n Make sure all fasta headers are unique and it is in fasta format! \n AUTOFIX: Taking last occurence of duplicates to continue analysis", |
| 78 | + e, |
| 79 | + ) |
| 80 | + ref_records = to_dict_remove_dups(SeqIO.parse(references, "fasta")) |
| 81 | + with open(f"{prefix}_reference.fa", "w") as f: |
| 82 | + init_position = f.tell() |
| 83 | + for hit in df["query-ID"].unique(): |
| 84 | + hit_name = hit.split(" ")[0] |
| 85 | + if hit_name in ref_records: |
| 86 | + # Sometimes reads can have illegal characters in the header |
| 87 | + ref_records[hit_name].id = ref_records[hit_name].id.replace("\\","_") |
| 88 | + ref_records[hit_name].description = ref_records[hit_name].description.replace("\\","_") |
| 89 | + SeqIO.write(ref_records[hit_name], f, "fasta") |
| 90 | + if f.tell() == init_position: |
| 91 | + logger.error("No reference sequences found in the hits. Exiting...") |
| 92 | + |
| 93 | + |
| 94 | +def read_mash_screen(file): |
| 95 | + """ |
| 96 | + Read in the file and return a pandas DataFrame |
| 97 | + File format: |
| 98 | + [identity, shared-hashes, median-multiplicity, p-value, query-ID, query-comment] |
| 99 | + 0.996341 3786/4000 121 0 USANYPRL230901_81A112023 |
| 100 | + 0.997144 3832/4000 124 0 EnglandQEUH3267E6482022 |
| 101 | + 0.997039 3826/4000 121 0 OP958840 |
| 102 | + 0.997022 3825/4000 122 0 OP971202 |
| 103 | + """ |
| 104 | + |
| 105 | + logger.info("Reading in the mash screen file...") |
| 106 | + df = pd.read_csv(file, sep="\t", header=None) |
| 107 | + df.columns = ["identity", "shared-hashes", "median-multiplicity", "p-value", "query-ID", "query-comment"] |
| 108 | + |
| 109 | + logger.info("Removing duplicates and sorting by identity and shared-hashes...") |
| 110 | + df['shared-hashes_num'] = df['shared-hashes'].str.split('/').str[0].astype(float) |
| 111 | + df = df.sort_values(by=["identity", "shared-hashes_num"], ascending=False) |
| 112 | + df = df.drop(columns=['shared-hashes_num']) |
| 113 | + |
| 114 | + return df.iloc[[0]] |
| 115 | + |
| 116 | + |
| 117 | +def main(argv=None): |
| 118 | + """Coordinate argument parsing and program execution.""" |
| 119 | + args = parse_args(argv) |
| 120 | + logging.basicConfig(level=args.log_level, format="[%(levelname)s] %(message)s") |
| 121 | + if not args.mash.is_file(): |
| 122 | + logger.error(f"The given input file {args.mash} was not found!") |
| 123 | + sys.exit(2) |
| 124 | + if not args.references.is_file(): |
| 125 | + logger.error(f"The given input file {args.references} was not found!") |
| 126 | + sys.exit(2) |
| 127 | + |
| 128 | + df = read_mash_screen(args.mash) |
| 129 | + |
| 130 | + extract_hits(df, args.references, args.prefix) |
| 131 | + |
| 132 | + df.to_json(f"{args.prefix}.json",orient="records", lines=True) |
| 133 | + |
| 134 | + return 0 |
| 135 | + |
| 136 | +if __name__ == "__main__": |
| 137 | + sys.exit(main()) |
0 commit comments