|
| 1 | +#!/usr/bin/env python |
| 2 | +import logging |
| 3 | +import argparse |
| 4 | +import re |
| 5 | +import statistics |
| 6 | +from typing import Set |
| 7 | + |
| 8 | +# Create a logger |
| 9 | +logging.basicConfig(format="%(name)s - %(asctime)s %(levelname)s: %(message)s") |
| 10 | +logger = logging.getLogger("fasta_gtf_filter") |
| 11 | +logger.setLevel(logging.INFO) |
| 12 | + |
| 13 | + |
| 14 | +def extract_fasta_seq_names(fasta_name: str) -> Set[str]: |
| 15 | + """Extracts the sequence names from a FASTA file.""" |
| 16 | + with open(fasta_name) as fasta: |
| 17 | + return {line[1:].split(None, 1)[0] for line in fasta if line.startswith(">")} |
| 18 | + |
| 19 | + |
| 20 | +def tab_delimited(file: str) -> float: |
| 21 | + """Check if file is tab-delimited and return median number of tabs.""" |
| 22 | + with open(file, "r") as f: |
| 23 | + data = f.read(1024) |
| 24 | + return statistics.median(line.count("\t") for line in data.split("\n")) |
| 25 | + |
| 26 | + |
| 27 | +def filter_gtf(fasta: str, gtf_in: str, filtered_gtf_out: str, skip_transcript_id_check: bool) -> None: |
| 28 | + """Filter GTF file based on FASTA sequence names.""" |
| 29 | + if tab_delimited(gtf_in) != 8: |
| 30 | + raise ValueError("Invalid GTF file: Expected 8 tab-separated columns.") |
| 31 | + |
| 32 | + seq_names_in_genome = extract_fasta_seq_names(fasta) |
| 33 | + logger.info(f"Extracted chromosome sequence names from {fasta}") |
| 34 | + logger.debug("All sequence IDs from FASTA: " + ", ".join(sorted(seq_names_in_genome))) |
| 35 | + |
| 36 | + seq_names_in_gtf = set() |
| 37 | + try: |
| 38 | + with open(gtf_in) as gtf, open(filtered_gtf_out, "w") as out: |
| 39 | + line_count = 0 |
| 40 | + for line in gtf: |
| 41 | + seq_name = line.split("\t")[0] |
| 42 | + seq_names_in_gtf.add(seq_name) # Add sequence name to the set |
| 43 | + |
| 44 | + if seq_name in seq_names_in_genome: |
| 45 | + if skip_transcript_id_check or re.search(r'transcript_id "([^"]+)"', line): |
| 46 | + out.write(line) |
| 47 | + line_count += 1 |
| 48 | + |
| 49 | + if line_count == 0: |
| 50 | + raise ValueError("All GTF lines removed by filters") |
| 51 | + |
| 52 | + except IOError as e: |
| 53 | + logger.error(f"File operation failed: {e}") |
| 54 | + return |
| 55 | + |
| 56 | + logger.debug("All sequence IDs from GTF: " + ", ".join(sorted(seq_names_in_gtf))) |
| 57 | + logger.info(f"Extracted {line_count} matching sequences from {gtf_in} into {filtered_gtf_out}") |
| 58 | + |
| 59 | + |
| 60 | +if __name__ == "__main__": |
| 61 | + parser = argparse.ArgumentParser(description="Filters a GTF file based on sequence names in a FASTA file.") |
| 62 | + parser.add_argument("--gtf", type=str, required=True, help="GTF file") |
| 63 | + parser.add_argument("--fasta", type=str, required=True, help="Genome fasta file") |
| 64 | + parser.add_argument("--prefix", dest="prefix", default="genes", type=str, help="Prefix for output GTF files") |
| 65 | + parser.add_argument( |
| 66 | + "--skip_transcript_id_check", action="store_true", help="Skip checking for transcript IDs in the GTF file" |
| 67 | + ) |
| 68 | + |
| 69 | + args = parser.parse_args() |
| 70 | + filter_gtf(args.fasta, args.gtf, args.prefix + ".filtered.gtf", args.skip_transcript_id_check) |
0 commit comments