Skip to content

Commit 8dea52f

Browse files
committed
.
1 parent 48abd57 commit 8dea52f

1 file changed

Lines changed: 51 additions & 6 deletions

File tree

virasign/virasign.py

Lines changed: 51 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -3706,8 +3706,10 @@ def extract_accessions_from_fasta(database_fasta: Path) -> set:
37063706

37073707
def build_header_mapping(database_fasta: Path) -> dict:
37083708
"""
3709-
Build a mapping from accession to full header description.
3710-
This helps match SAM headers (which may be truncated) to full FASTA headers.
3709+
Build a mapping from accession to FASTA title (text after '>', stripped).
3710+
3711+
For selected-reference FASTAs this title is usually the accession only, so the map
3712+
doubles as accession -> display name for remapped SAM stats.
37113713
"""
37123714
header_map = {}
37133715

@@ -4090,6 +4092,29 @@ def _canonical_curated_accession(to_canon: Dict[str, str], raw: str) -> str:
40904092
return to_canon.get(base, raw)
40914093

40924094

4095+
def _fasta_seq_name_for_export(header: str, fallback: str = "") -> str:
4096+
"""
4097+
Single-word FASTA name for remapping and exported reference FASTAs (consensus-friendly).
4098+
4099+
Uses NCBI-style accession when parsable; avoids pipe/colon-heavy DB headers in SAM/BAM SN.
4100+
"""
4101+
h = (header or "").strip()
4102+
acc = extract_accession_from_header(h)
4103+
if acc:
4104+
return acc
4105+
if h:
4106+
first = h.split()[0].strip()
4107+
if first:
4108+
return first
4109+
fb = (fallback or "").strip()
4110+
if fb:
4111+
acc2 = extract_accession_from_header(fb)
4112+
if acc2:
4113+
return acc2
4114+
return fb.split()[0].strip() or "reference"
4115+
return "reference"
4116+
4117+
40934118
def extract_selected_references(database_fasta: Path, selected_headers: list, out_fasta: Path) -> int:
40944119
"""
40954120
Extract multiple reference sequences from database FASTA matching the given headers.
@@ -4134,7 +4159,8 @@ def extract_selected_references(database_fasta: Path, selected_headers: list, ou
41344159
header_acc_base = header_acc.split('.')[0] if '.' in header_acc else header_acc
41354160
write = (header_acc in selected_accessions_set or header_acc_base in selected_accessions_set)
41364161
if write:
4137-
out.write(line)
4162+
seq_name = _fasta_seq_name_for_export(header)
4163+
out.write(f">{seq_name}\n")
41384164
found_count += 1
41394165
else:
41404166
if write:
@@ -4162,7 +4188,8 @@ def _updated_stats_from_remapped_sam(
41624188
min_identity: float,
41634189
) -> Dict[str, dict]:
41644190
"""
4165-
Parse remapped SAM (same alignment set as sorted BAM) into stats keyed by full FASTA description.
4191+
Parse remapped SAM (same alignment set as sorted BAM) into stats keyed by FASTA title
4192+
(typically the accession when references were exported with consensus-friendly headers).
41664193
"""
41674194
selected_header_map = build_header_mapping(selected_refs_fasta)
41684195
for acc_key, hdr in list(selected_header_map.items()):
@@ -6630,7 +6657,8 @@ def extract_fasta_record(database_fasta: Path, target_header: str, out_fasta: Pa
66306657
# Match on full header
66316658
write = (header == target_header)
66326659
if write:
6633-
out.write(line)
6660+
seq_name = _fasta_seq_name_for_export(header, fallback=target_header)
6661+
out.write(f">{seq_name}\n")
66346662
found = True
66356663
else:
66366664
if write:
@@ -6665,7 +6693,8 @@ def extract_fasta_record_by_accession(database_fasta: Path, accession: str, out_
66656693
header_acc = extract_accession_from_header(header)
66666694
write = (header_acc == accession)
66676695
if write:
6668-
out.write(line)
6696+
# Use curated accession so folder name, FASTA, and BAM RNAME stay identical.
6697+
out.write(f">{accession}\n")
66696698
found = True
66706699
else:
66716700
if write:
@@ -11438,6 +11467,22 @@ def main(args=None):
1143811467
for database_fasta_path in database_fasta_paths:
1143911468
prepare_database_context(Path(database_fasta_path))
1144011469

11470+
# Taxonomy SQLite lives under each database directory. Without this, parallel sample
11471+
# workers each call ensure_taxonomy_resources() and can race: duplicate NCBI dump
11472+
# parsing ("Processed ... lines... (kept ... virus accessions)" interleaved in logs).
11473+
# RVDB/RefSeq FASTA download still happens only in resolve_database_path() above.
11474+
logger.info("Ensuring taxonomy resources (parent process, once per database)...")
11475+
for database_fasta_path in database_fasta_paths:
11476+
try:
11477+
ensure_taxonomy_resources(
11478+
Path(database_fasta_path),
11479+
force_rebuild=getattr(args, "rebuild", False),
11480+
)
11481+
except Exception as e:
11482+
logger.warning(
11483+
f"Taxonomy preparation failed for {database_fasta_path} (workers may retry): {e}"
11484+
)
11485+
1144111486
# Prepare DB only mode: download/unpack/index completed, no sample processing.
1144211487
if getattr(args, "prepare_db", False):
1144311488
print("Database preparation complete.")

0 commit comments

Comments
 (0)