Skip to content

Commit 35a6007

Browse files
committed
.
1 parent e117f8d commit 35a6007

2 files changed

Lines changed: 89294 additions & 53437 deletions

File tree

recmpox/recmpox.py

Lines changed: 16 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1380,14 +1380,21 @@ def _extract_tract_sequences(
13801380

13811381
# ref1 (Ia) file: keep only ia tract bases; Ib + other → N.
13821382
# ref2 (Ib) file: keep only ib tract bases; Ia + other → N.
1383+
safe_id = _safe_fasta_id(sample_id)
13831384
seq1 = _extract_tracts_as_n_full_length(aligned_seq, merged_tracts, keep_clade="ia")
1384-
fh1.write(f">{sample_id} {ref1_label}_tracts_only\n")
1385-
for i in range(0, len(seq1), line_len):
1385+
len1 = len(seq1)
1386+
non_n1 = sum(1 for b in seq1.upper() if b in "ACGT")
1387+
cov1 = (100.0 * non_n1 / len1) if len1 else 0.0
1388+
fh1.write(f">{safe_id}_{ref1_label}_tract_HC_{cov1:.2f}%\n")
1389+
for i in range(0, len1, line_len):
13861390
fh1.write(seq1[i : i + line_len] + "\n")
13871391

13881392
seq2 = _extract_tracts_as_n_full_length(aligned_seq, merged_tracts, keep_clade="ib")
1389-
fh2.write(f">{sample_id} {ref2_label}_tracts_only\n")
1390-
for i in range(0, len(seq2), line_len):
1393+
len2 = len(seq2)
1394+
non_n2 = sum(1 for b in seq2.upper() if b in "ACGT")
1395+
cov2 = (100.0 * non_n2 / len2) if len2 else 0.0
1396+
fh2.write(f">{safe_id}_{ref2_label}_tract_HC_{cov2:.2f}%\n")
1397+
for i in range(0, len2, line_len):
13911398
fh2.write(seq2[i : i + line_len] + "\n")
13921399

13931400
n_written += 1
@@ -1432,7 +1439,7 @@ def _run_phylogeny_pipeline(
14321439
logger.warning("--phylogeny: extracted tract FASTAs not found (%s, %s); skipping phylogeny.", ref1_fa, ref2_fa)
14331440
return
14341441

1435-
# Use only the bundled reference set (no user override)
1442+
# Use bundled reference set (no clade-based restriction)
14361443
refs_path = PHYLOGENY_REFS_FASTA
14371444
if not refs_path.exists():
14381445
logger.error("--phylogeny: bundled references not found at %s", refs_path)
@@ -1499,13 +1506,15 @@ def _run_phylogeny_pipeline(
14991506
cmd_iqtree = [
15001507
"iqtree",
15011508
"-s", str(aln_in_phylogeny),
1502-
"-m", "GTR",
15031509
"-bb", "10000",
15041510
"-pre", str(iqtree_prefix),
15051511
"-czb",
15061512
]
1507-
if getattr(args, "threads", 1) and int(args.threads) > 1:
1513+
# Use all available CPUs by default; if user requests >1 threads, honor it.
1514+
if getattr(args, "threads", None) is not None and int(args.threads) > 1:
15081515
cmd_iqtree.extend(["-nt", str(args.threads)])
1516+
else:
1517+
cmd_iqtree.extend(["-nt", "AUTO"])
15091518
logger.info("Running IQ-TREE: %s", " ".join(cmd_iqtree))
15101519
result = subprocess.run(cmd_iqtree, timeout=7200)
15111520
if result.returncode != 0:

0 commit comments

Comments
 (0)