scripts/sequence_inspector.py
8.4 KB| 1 | #!/usr/bin/env python3 |
| 2 | """Bounded FASTA/FASTQ aggregate inspector with no identifier or sequence output.""" |
| 3 | |
| 4 | from __future__ import annotations |
| 5 | |
| 6 | import argparse |
| 7 | import math |
| 8 | from pathlib import Path |
| 9 | from typing import Any |
| 10 | |
| 11 | from _capabilities import capability_for_path, validate_magic |
| 12 | from _common import ( |
| 13 | DEFAULT_MAX_FILE_BYTES, |
| 14 | CliError, |
| 15 | bounded_file_limit, |
| 16 | bounded_integer, |
| 17 | checked_input_file, |
| 18 | emit_json, |
| 19 | run_cli, |
| 20 | stable_token, |
| 21 | ) |
| 22 | |
| 23 | |
| 24 | DEFAULT_MAX_RECORDS = 10_000 |
| 25 | MAX_RECORDS = 100_000 |
| 26 | DEFAULT_MAX_BASES = 10_000_000 |
| 27 | MAX_BASES = 100_000_000 |
| 28 | NUCLEOTIDE_CODES = frozenset("ACGTUNRYSWKMBDHVX-.") |
| 29 | |
| 30 | |
| 31 | def build_parser() -> argparse.ArgumentParser: |
| 32 | parser = argparse.ArgumentParser( |
| 33 | description=( |
| 34 | "Inspect a bounded local FASTA/FASTQ sample. Sequence text and record " |
| 35 | "identifiers are never emitted." |
| 36 | ) |
| 37 | ) |
| 38 | parser.add_argument("input", help="Local FASTA/FASTQ path inside --root") |
| 39 | parser.add_argument( |
| 40 | "--root", |
| 41 | default=".", |
| 42 | help="Existing local directory that bounds all input/output paths", |
| 43 | ) |
| 44 | parser.add_argument( |
| 45 | "--max-bytes", |
| 46 | type=int, |
| 47 | default=DEFAULT_MAX_FILE_BYTES, |
| 48 | help=f"Maximum input bytes (hard ceiling: {512 * 1024 * 1024})", |
| 49 | ) |
| 50 | parser.add_argument( |
| 51 | "--max-records", |
| 52 | type=int, |
| 53 | default=DEFAULT_MAX_RECORDS, |
| 54 | help=f"Maximum records to inspect (hard ceiling: {MAX_RECORDS})", |
| 55 | ) |
| 56 | parser.add_argument( |
| 57 | "--max-bases", |
| 58 | type=int, |
| 59 | default=DEFAULT_MAX_BASES, |
| 60 | help=f"Maximum sequence characters to inspect (hard ceiling: {MAX_BASES})", |
| 61 | ) |
| 62 | parser.add_argument("--output", help="Optional local .json output path") |
| 63 | parser.add_argument( |
| 64 | "--force", |
| 65 | action="store_true", |
| 66 | help="Allow replacement of an existing regular output file", |
| 67 | ) |
| 68 | return parser |
| 69 | |
| 70 | |
| 71 | def inspect_sequence_file( |
| 72 | path: Path, |
| 73 | *, |
| 74 | suffix: str, |
| 75 | max_records: int = DEFAULT_MAX_RECORDS, |
| 76 | max_bases: int = DEFAULT_MAX_BASES, |
| 77 | ) -> dict[str, Any]: |
| 78 | """Inspect records with Biopython's streaming low-level parsers.""" |
| 79 | |
| 80 | max_records = bounded_integer( |
| 81 | max_records, |
| 82 | name="max records", |
| 83 | minimum=1, |
| 84 | maximum=MAX_RECORDS, |
| 85 | ) |
| 86 | max_bases = bounded_integer( |
| 87 | max_bases, |
| 88 | name="max bases", |
| 89 | minimum=1, |
| 90 | maximum=MAX_BASES, |
| 91 | ) |
| 92 | try: |
| 93 | from Bio.SeqIO.FastaIO import SimpleFastaParser |
| 94 | from Bio.SeqIO.QualityIO import FastqGeneralIterator |
| 95 | except ImportError as exc: |
| 96 | raise CliError( |
| 97 | 'optional dependency missing; install with: uv pip install "biopython==1.88"' |
| 98 | ) from exc |
| 99 | is_fastq = suffix in {".fastq", ".fq"} |
| 100 | lengths: list[int] = [] |
| 101 | total_bases = 0 |
| 102 | nucleotide_like_bases = 0 |
| 103 | gc_bases = 0 |
| 104 | ambiguous_bases = 0 |
| 105 | identifier_tokens: set[str] = set() |
| 106 | duplicate_identifier_count = 0 |
| 107 | quality_count = 0 |
| 108 | quality_sum = 0 |
| 109 | quality_minimum: int | None = None |
| 110 | quality_maximum: int | None = None |
| 111 | truncated = False |
| 112 | try: |
| 113 | with path.open("r", encoding="ascii", errors="strict", newline=None) as handle: |
| 114 | iterator = ( |
| 115 | FastqGeneralIterator(handle) |
| 116 | if is_fastq |
| 117 | else SimpleFastaParser(handle) |
| 118 | ) |
| 119 | for record_index, record in enumerate(iterator): |
| 120 | if record_index >= max_records: |
| 121 | truncated = True |
| 122 | break |
| 123 | title = record[0] |
| 124 | sequence = record[1].upper() |
| 125 | if total_bases + len(sequence) > max_bases: |
| 126 | truncated = True |
| 127 | break |
| 128 | token = stable_token(title, kind="sequence_id") |
| 129 | if token in identifier_tokens: |
| 130 | duplicate_identifier_count += 1 |
| 131 | identifier_tokens.add(token) |
| 132 | lengths.append(len(sequence)) |
| 133 | total_bases += len(sequence) |
| 134 | valid_nucleotide = sum( |
| 135 | character in NUCLEOTIDE_CODES for character in sequence |
| 136 | ) |
| 137 | nucleotide_like_bases += valid_nucleotide |
| 138 | gc_bases += sequence.count("G") + sequence.count("C") |
| 139 | ambiguous_bases += sum( |
| 140 | character not in {"A", "C", "G", "T", "U"} |
| 141 | for character in sequence |
| 142 | ) |
| 143 | if is_fastq: |
| 144 | quality = record[2] |
| 145 | scores = [ord(character) - 33 for character in quality] |
| 146 | if any(score < 0 for score in scores): |
| 147 | raise CliError("FASTQ contains a quality character below Phred+33") |
| 148 | if scores: |
| 149 | quality_count += len(scores) |
| 150 | quality_sum += sum(scores) |
| 151 | local_minimum = min(scores) |
| 152 | local_maximum = max(scores) |
| 153 | quality_minimum = ( |
| 154 | local_minimum |
| 155 | if quality_minimum is None |
| 156 | else min(quality_minimum, local_minimum) |
| 157 | ) |
| 158 | quality_maximum = ( |
| 159 | local_maximum |
| 160 | if quality_maximum is None |
| 161 | else max(quality_maximum, local_maximum) |
| 162 | ) |
| 163 | except CliError: |
| 164 | raise |
| 165 | except (OSError, UnicodeError, ValueError) as exc: |
| 166 | raise CliError("the sequence file could not be parsed safely") from exc |
| 167 | count = len(lengths) |
| 168 | nucleotide_fraction = ( |
| 169 | nucleotide_like_bases / total_bases if total_bases else None |
| 170 | ) |
| 171 | nucleotide_like = bool( |
| 172 | nucleotide_fraction is not None and nucleotide_fraction >= 0.95 |
| 173 | ) |
| 174 | report: dict[str, Any] = { |
| 175 | "profile_type": "fastq_bounded_profile" |
| 176 | if is_fastq |
| 177 | else "fasta_bounded_profile", |
| 178 | "records_inspected": count, |
| 179 | "sequence_characters_inspected": total_bases, |
| 180 | "record_or_base_limit_reached": truncated, |
| 181 | "length_aggregates": { |
| 182 | "minimum": min(lengths) if lengths else None, |
| 183 | "maximum": max(lengths) if lengths else None, |
| 184 | "mean": math.fsum(lengths) / count if count else None, |
| 185 | }, |
| 186 | "duplicate_identifier_token_count": duplicate_identifier_count, |
| 187 | "record_identifiers_emitted": False, |
| 188 | "sequence_values_emitted": False, |
| 189 | "nucleotide_like_fraction": nucleotide_fraction, |
| 190 | "interpreted_as_nucleotide": nucleotide_like, |
| 191 | "gc_fraction_if_nucleotide": ( |
| 192 | gc_bases / total_bases if nucleotide_like and total_bases else None |
| 193 | ), |
| 194 | "ambiguous_character_fraction_if_nucleotide": ( |
| 195 | ambiguous_bases / total_bases |
| 196 | if nucleotide_like and total_bases |
| 197 | else None |
| 198 | ), |
| 199 | "limitations": [ |
| 200 | "Counts describe only the bounded inspected prefix when a limit is reached.", |
| 201 | "Alphabet classification is heuristic and does not infer organism or molecule identity.", |
| 202 | ], |
| 203 | } |
| 204 | if is_fastq: |
| 205 | report["quality_aggregates"] = { |
| 206 | "encoding_assumption": "Phred+33", |
| 207 | "score_count": quality_count, |
| 208 | "minimum": quality_minimum, |
| 209 | "maximum": quality_maximum, |
| 210 | "mean": quality_sum / quality_count if quality_count else None, |
| 211 | "requires_encoding_confirmation": True, |
| 212 | } |
| 213 | return report |
| 214 | |
| 215 | |
| 216 | def _main() -> None: |
| 217 | args = build_parser().parse_args() |
| 218 | max_bytes = bounded_file_limit(args.max_bytes) |
| 219 | path = checked_input_file( |
| 220 | args.input, |
| 221 | root=args.root, |
| 222 | suffixes={".fasta", ".fa", ".fna", ".fastq", ".fq"}, |
| 223 | max_bytes=max_bytes, |
| 224 | ) |
| 225 | capability = capability_for_path(path) |
| 226 | validate_magic(path, capability["suffix"]) |
| 227 | report = { |
| 228 | "schema_version": "1.1", |
| 229 | "capability": capability, |
| 230 | "analysis": inspect_sequence_file( |
| 231 | path, |
| 232 | suffix=capability["suffix"], |
| 233 | max_records=args.max_records, |
| 234 | max_bases=args.max_bases, |
| 235 | ), |
| 236 | "security": { |
| 237 | "local_only": True, |
| 238 | "untrusted_text_never_treated_as_instructions": True, |
| 239 | "raw_values_and_identifiers_emitted": False, |
| 240 | }, |
| 241 | } |
| 242 | emit_json( |
| 243 | report, |
| 244 | output=args.output, |
| 245 | root=args.root, |
| 246 | force=args.force, |
| 247 | ) |
| 248 | |
| 249 | |
| 250 | def main() -> int: |
| 251 | return run_cli(_main) |
| 252 | |
| 253 | |
| 254 | if __name__ == "__main__": |
| 255 | raise SystemExit(main()) |
文件内容是平台检查时读出的文本，只供查看。平台不会运行包里的任何脚本或程序。
