|
| 1 | +#!/usr/bin/env python3 |
| 2 | +import os |
| 3 | +import re |
| 4 | +import sys |
| 5 | +from pathlib import Path |
| 6 | + |
| 7 | + |
| 8 | +def read_fasta(path, sample_id): |
| 9 | + entries = [] |
| 10 | + header = None |
| 11 | + seq_lines = [] |
| 12 | + with open(path, "r") as handle: |
| 13 | + for raw in handle: |
| 14 | + line = raw.strip() |
| 15 | + if not line: |
| 16 | + continue |
| 17 | + if line.startswith(">"): |
| 18 | + if header is not None: |
| 19 | + entries.append((header, "".join(seq_lines).upper())) |
| 20 | + header = line[1:].strip() or f"{sample_id}_chain_{len(entries) + 1}" |
| 21 | + seq_lines = [] |
| 22 | + else: |
| 23 | + seq_lines.append(line.replace(" ", "").upper()) |
| 24 | + if header is not None: |
| 25 | + entries.append((header, "".join(seq_lines).upper())) |
| 26 | + return entries |
| 27 | + |
| 28 | + |
| 29 | +def infer_type(header, sequence): |
| 30 | + type_aliases = { |
| 31 | + "protein": "P", |
| 32 | + "prot": "P", |
| 33 | + "aa": "P", |
| 34 | + "pep": "P", |
| 35 | + "peptide": "P", |
| 36 | + "p": "P", |
| 37 | + "rna": "R", |
| 38 | + "r": "R", |
| 39 | + "double": "D", |
| 40 | + "ds": "D", |
| 41 | + "dsdna": "D", |
| 42 | + "double_dna": "D", |
| 43 | + "single": "S", |
| 44 | + "ss": "S", |
| 45 | + "ssdna": "S", |
| 46 | + "single_dna": "S", |
| 47 | + "single-strand": "S", |
| 48 | + "singlestrand": "S", |
| 49 | + } |
| 50 | + header_lower = header.lower() |
| 51 | + match = re.search(r"(?:type|entity|molecule|mol)[:=]\s*([A-Za-z0-9_-]+)", header_lower) |
| 52 | + if match: |
| 53 | + candidate = match.group(1).lower() |
| 54 | + if candidate in type_aliases: |
| 55 | + return type_aliases[candidate] |
| 56 | + for alias, code in type_aliases.items(): |
| 57 | + if re.search(r"\b" + re.escape(alias) + r"\b", header_lower): |
| 58 | + return code |
| 59 | + |
| 60 | + seq_set = set(sequence) |
| 61 | + if not sequence: |
| 62 | + return None |
| 63 | + if seq_set <= set("ACUGN"): |
| 64 | + return "R" |
| 65 | + # Default DNA to double-stranded unless explicitly marked single-strand. |
| 66 | + if seq_set <= set("ACTGN"): |
| 67 | + return "D" |
| 68 | + protein_letters = set("ACDEFGHIKLMNPQRSTVWYBXZOU") |
| 69 | + if seq_set <= protein_letters and not (seq_set <= set("ACUGTN")): |
| 70 | + return "P" |
| 71 | + return "P" |
| 72 | + |
| 73 | + |
| 74 | +def main(): |
| 75 | + if len(sys.argv) != 3: |
| 76 | + sys.stderr.write("Usage: fasta_to_rosettafold.py <sample_id> <fasta_path>\n") |
| 77 | + return 1 |
| 78 | + |
| 79 | + sample_id, fasta_path = sys.argv[1], sys.argv[2] |
| 80 | + allowed_ext = (".fa", ".fasta", ".fas", ".faa", ".fna") |
| 81 | + if not fasta_path.lower().endswith(allowed_ext): |
| 82 | + sys.stderr.write( |
| 83 | + f"[ROSETTAFOLD2NA_FASTA] Input file '{fasta_path}' must be a FASTA file.\n" |
| 84 | + ) |
| 85 | + return 1 |
| 86 | + |
| 87 | + if not os.path.exists(fasta_path): |
| 88 | + sys.stderr.write( |
| 89 | + f"[ROSETTAFOLD2NA_FASTA] Input FASTA '{fasta_path}' does not exist.\n" |
| 90 | + ) |
| 91 | + return 1 |
| 92 | + |
| 93 | + entries = read_fasta(fasta_path, sample_id) |
| 94 | + if not entries: |
| 95 | + sys.stderr.write( |
| 96 | + f"[ROSETTAFOLD2NA_FASTA] No sequences found in '{fasta_path}'.\n" |
| 97 | + ) |
| 98 | + return 1 |
| 99 | + |
| 100 | + output_dir = Path("rf2na_input") |
| 101 | + output_dir.mkdir(parents=True, exist_ok=True) |
| 102 | + |
| 103 | + chain_records = [] |
| 104 | + observed_files = set() |
| 105 | + for idx, (header, sequence) in enumerate(entries, start=1): |
| 106 | + chain_type = infer_type(header, sequence) |
| 107 | + if chain_type is None: |
| 108 | + sys.stderr.write( |
| 109 | + f"[ROSETTAFOLD2NA_FASTA] Unable to determine entity type for entry '{header}'. " |
| 110 | + "Please include a token such as 'type=protein', 'type=double_dna', or 'type=single_dna'.\n" |
| 111 | + ) |
| 112 | + return 1 |
| 113 | + if chain_type not in {"P", "R", "D", "S"}: |
| 114 | + sys.stderr.write( |
| 115 | + f"[ROSETTAFOLD2NA_FASTA] Unable to determine entity type for entry '{header}'. " |
| 116 | + "Allowed types: protein (P), rna (R), double_dna (D), single_dna (S).\n" |
| 117 | + ) |
| 118 | + return 1 |
| 119 | + safe_name = re.sub(r"[^A-Za-z0-9_.-]+", "_", header) or f"chain_{idx}" |
| 120 | + filename = f"chain_{idx:03d}_{safe_name[:40]}.fa" |
| 121 | + if filename in observed_files: |
| 122 | + filename = f"chain_{idx:03d}_{idx}.fa" |
| 123 | + observed_files.add(filename) |
| 124 | + with open(output_dir / filename, "w") as fh: |
| 125 | + fh.write(f">{header}\n") |
| 126 | + for start in range(0, len(sequence), 80): |
| 127 | + fh.write(sequence[start : start + 80] + "\n") |
| 128 | + chain_records.append((chain_type, filename, header)) |
| 129 | + |
| 130 | + with open(output_dir / "chain_map.tsv", "w") as mapping: |
| 131 | + mapping.write("type\tfilename\theader\n") |
| 132 | + for chain_type, filename, header in chain_records: |
| 133 | + mapping.write(f"{chain_type}\t{filename}\t{header}\n") |
| 134 | + |
| 135 | + return 0 |
| 136 | + |
| 137 | + |
| 138 | +if __name__ == "__main__": |
| 139 | + sys.exit(main()) |
0 commit comments