Skip to content

Commit b7df8b7

Browse files
committed
.
1 parent 4a00716 commit b7df8b7

2 files changed

Lines changed: 167 additions & 37 deletions

File tree

README.md

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -94,7 +94,7 @@ virasign --help
9494
- `-o, --output`: Output directory (default: creates `Virasign_output/`).
9595

9696
- **Choose database (auto-downloads on first run)**
97-
- `-d, --database`: Database name `RVDB,RefSeq`, or an accession (e.g. `OZ254622.1`) (default: `RVDB`).
97+
- `-d, --database`: `RVDB` (default), `RefSeq`, `RVDB,RefSeq`, an accession (e.g. `OZ254622.1`), or a species name (e.g. `Orthopoxvirus monkeypox`).
9898
- `--rvdb-version`: Which RVDB release to download (default: `31.0`). See [available versions](https://rvdb.dbi.udel.edu/previous-release).
9999
- `-a, --accession`: Extra NCBI accessions to include in the run (merged with selected database).
100100
- `--db-dir`: Reuse an existing database folder (optional; example: `/path/to/Databases/`).
@@ -144,6 +144,13 @@ virasign -i input_dir -d RVDB,RefSeq --rvdb-version 31.0
144144
# Use a single accession as the database
145145
virasign -i input_dir -d OZ254622.1 -o output_dir
146146

147+
# Use an organism/species-restricted database (downloads a small custom database)
148+
virasign -i input_dir -d "Orthopoxvirus monkeypox" -o output_dir
149+
150+
# Use text file with species names as database
151+
virasign -i input_dir -d species_list.txt -o output_dir
152+
# (species_list.txt contains one species name per line)
153+
147154
# Use text file with accessions as database
148155
virasign -i input_dir -d my_accessions.txt -o output_dir
149156
# (my_accessions.txt contains one accession per line)

virasign/virasign.py

Lines changed: 159 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -1217,6 +1217,106 @@ def download_accession_from_ncbi(accession: str, output_dir: Path = None) -> Pat
12171217
logger.error(f"Failed to download accession {accession}: {e}")
12181218
raise
12191219

1220+
1221+
def download_species_database_from_ncbi(
1222+
species_name: str,
1223+
output_fasta: Path,
1224+
retmax: int = 5000,
1225+
batch_size: int = 500,
1226+
) -> Path:
1227+
"""
1228+
Download a FASTA database for a given viral species/organism name from NCBI nuccore.
1229+
1230+
Uses E-utilities esearch (usehistory) + efetch in batches.
1231+
The search term is restricted to complete genomes/sequences to keep DB size manageable.
1232+
"""
1233+
species_name = (species_name or "").strip().strip('"').strip("'")
1234+
if not species_name:
1235+
raise ValueError("Species name is empty.")
1236+
1237+
output_fasta = Path(output_fasta)
1238+
output_fasta.parent.mkdir(parents=True, exist_ok=True)
1239+
1240+
term = (
1241+
f"\"{species_name}\"[Organism] AND "
1242+
f"(\"complete genome\"[Title] OR \"complete sequence\"[Title] OR complete[Title])"
1243+
)
1244+
logger.info(f"Searching NCBI for species database: {species_name}")
1245+
1246+
search_params = {
1247+
"db": "nuccore",
1248+
"term": term,
1249+
"retmode": "json",
1250+
"usehistory": "y",
1251+
"retmax": str(int(retmax)),
1252+
}
1253+
search_url = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi"
1254+
full_search_url = f"{search_url}?{urllib.parse.urlencode(search_params)}"
1255+
1256+
with urllib.request.urlopen(full_search_url, timeout=60) as resp:
1257+
data = json.loads(resp.read().decode("utf-8"))
1258+
1259+
esr = data.get("esearchresult", {})
1260+
count = int(esr.get("count", 0) or 0)
1261+
webenv = esr.get("webenv")
1262+
query_key = esr.get("querykey")
1263+
idlist = esr.get("idlist") or []
1264+
1265+
if count == 0:
1266+
raise ValueError(f"No NCBI nuccore records found for species '{species_name}'.")
1267+
1268+
if count > int(retmax):
1269+
logger.warning(
1270+
f"NCBI search found {count} records but retmax={retmax}; "
1271+
f"database will be truncated to the first {retmax} records."
1272+
)
1273+
1274+
if not webenv or not query_key:
1275+
# Fallback: fetch by idlist (rare)
1276+
logger.warning("NCBI usehistory not available; falling back to direct id list fetch.")
1277+
ids = idlist
1278+
if not ids:
1279+
raise ValueError("NCBI search returned no IDs to fetch.")
1280+
fetch_params = {"db": "nuccore", "id": ",".join(ids), "rettype": "fasta", "retmode": "text"}
1281+
fetch_url = f"https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?{urllib.parse.urlencode(fetch_params)}"
1282+
tmp = output_fasta.with_suffix(output_fasta.suffix + ".tmp")
1283+
with urllib.request.urlopen(fetch_url, timeout=120) as resp, open(tmp, "wb") as out:
1284+
out.write(resp.read())
1285+
if tmp.stat().st_size == 0:
1286+
tmp.unlink()
1287+
raise ValueError("NCBI returned empty FASTA for species query.")
1288+
tmp.replace(output_fasta)
1289+
logger.info(f"Saved species database FASTA: {output_fasta}")
1290+
return output_fasta
1291+
1292+
fetch_url = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi"
1293+
tmp = output_fasta.with_suffix(output_fasta.suffix + ".tmp")
1294+
max_to_fetch = min(count, int(retmax))
1295+
with open(tmp, "wb") as out:
1296+
for retstart in range(0, max_to_fetch, int(batch_size)):
1297+
fetch_params = {
1298+
"db": "nuccore",
1299+
"query_key": str(query_key),
1300+
"WebEnv": str(webenv),
1301+
"rettype": "fasta",
1302+
"retmode": "text",
1303+
"retstart": str(retstart),
1304+
"retmax": str(int(batch_size)),
1305+
}
1306+
full_fetch_url = f"{fetch_url}?{urllib.parse.urlencode(fetch_params)}"
1307+
with urllib.request.urlopen(full_fetch_url, timeout=120) as resp:
1308+
chunk = resp.read()
1309+
if chunk:
1310+
out.write(chunk)
1311+
1312+
if tmp.stat().st_size == 0:
1313+
tmp.unlink()
1314+
raise ValueError("NCBI returned empty FASTA for species query.")
1315+
1316+
tmp.replace(output_fasta)
1317+
logger.info(f"Saved species database FASTA: {output_fasta}")
1318+
return output_fasta
1319+
12201320
def merge_fasta_files(fasta_files: list, output_fasta: Path) -> Path:
12211321
"""
12221322
Merge multiple FASTA files into a single file.
@@ -1287,6 +1387,25 @@ def resolve_database_path(
12871387
database_arg = database_arg.strip()
12881388
databases_dir = Path(databases_dir) if databases_dir is not None else get_virasign_databases_dir()
12891389
databases_dir.mkdir(parents=True, exist_ok=True)
1390+
1391+
def _resolve_species_database(species_name: str) -> Path:
1392+
species_name = (species_name or "").strip().strip('"').strip("'")
1393+
if not species_name:
1394+
raise ValueError("Species name is empty.")
1395+
species_dir = databases_dir / "Custom" / "Species"
1396+
species_dir.mkdir(parents=True, exist_ok=True)
1397+
out_fasta = species_dir / f"species_{safe_stem(species_name)}.fasta"
1398+
if out_fasta.exists() and not force_named_database_rebuild:
1399+
logger.info(f"Using existing species database: {out_fasta}")
1400+
return out_fasta
1401+
logger.info(f"Building species database for: {species_name}")
1402+
return download_species_database_from_ncbi(species_name, out_fasta)
1403+
1404+
# Species database shortcut: species:<name>
1405+
lowered = database_arg.lower()
1406+
if lowered.startswith("species:") or lowered.startswith("organism:"):
1407+
species_name = database_arg.split(":", 1)[1].strip()
1408+
return _resolve_species_database(species_name)
12901409

12911410
# Check if database_arg is a single accession number
12921411
if is_accession_number(database_arg):
@@ -1326,49 +1445,46 @@ def resolve_database_path(
13261445
# FASTA files typically have extensions: .fasta, .fa, .fna, .fas, .faa, .fq, .fastq
13271446
fasta_extensions = {'.fasta', '.fa', '.fna', '.fas', '.faa', '.fq', '.fastq', '.gz'}
13281447
if db_path.suffix.lower() not in fasta_extensions:
1329-
# It's likely a text file with accessions
1330-
logger.info(f"Reading accessions from file: {db_path}")
1331-
accession_list = []
1332-
with open(db_path, 'r') as f:
1448+
# Text file with either accessions OR species names (one per line, comments allowed)
1449+
logger.info(f"Reading database entries from file: {db_path}")
1450+
items = []
1451+
with open(db_path, "r") as f:
13331452
for line in f:
13341453
line = line.strip()
1335-
if line and not line.startswith('#'): # Skip empty lines and comments
1336-
accession_list.append(line)
1337-
1338-
if not accession_list:
1339-
raise ValueError(f"No accessions found in file: {db_path}")
1340-
1341-
logger.info(f"Found {len(accession_list)} accession(s) in file")
1342-
1343-
# Download all accessions to Custom directory
1454+
if not line or line.startswith("#"):
1455+
continue
1456+
items.append(line)
1457+
1458+
if not items:
1459+
raise ValueError(f"No database entries found in file: {db_path}")
1460+
1461+
accessions_in_file = [x for x in items if is_accession_number(x)]
1462+
species_in_file = [x for x in items if not is_accession_number(x)]
1463+
13441464
custom_dir = databases_dir / "Custom"
13451465
custom_dir.mkdir(parents=True, exist_ok=True)
1346-
1347-
accession_fasta_files = []
1348-
for accession in accession_list:
1349-
try:
1350-
acc_fasta = download_accession_from_ncbi(accession.strip(), custom_dir)
1351-
accession_fasta_files.append(acc_fasta)
1352-
except Exception as e:
1353-
logger.error(f"Failed to download accession {accession}: {e}")
1354-
raise
1355-
1466+
1467+
fasta_files = []
1468+
if accessions_in_file:
1469+
logger.info(f"Found {len(accessions_in_file)} accession(s) in file")
1470+
for accession in accessions_in_file:
1471+
fasta_files.append(download_accession_from_ncbi(accession.strip(), custom_dir))
1472+
1473+
if species_in_file:
1474+
logger.info(f"Found {len(species_in_file)} species/organism name(s) in file")
1475+
for sp in species_in_file:
1476+
fasta_files.append(_resolve_species_database(sp))
1477+
13561478
# If additional accessions provided via -a, add them
13571479
if accessions:
13581480
logger.info(f"Downloading {len(accessions)} additional accession(s) from -a argument...")
13591481
for accession in accessions:
1360-
try:
1361-
acc_fasta = download_accession_from_ncbi(accession.strip(), custom_dir)
1362-
accession_fasta_files.append(acc_fasta)
1363-
except Exception as e:
1364-
logger.error(f"Failed to download accession {accession}: {e}")
1365-
raise
1366-
1367-
# Merge all accessions
1482+
fasta_files.append(download_accession_from_ncbi(accession.strip(), custom_dir))
1483+
13681484
file_stem = db_path.stem
13691485
merged_fasta = custom_dir / f"{file_stem}_database.fasta"
1370-
result = merge_fasta_files(accession_fasta_files, merged_fasta)
1371-
logger.info(f"Created merged database from {len(accession_fasta_files)} accession(s): {merged_fasta}")
1486+
result = merge_fasta_files(fasta_files, merged_fasta)
1487+
logger.info(f"Created merged database from {len(fasta_files)} entry/entries: {merged_fasta}")
13721488
return result
13731489

13741490
# Download accessions if provided (will determine target directory after we know which database)
@@ -1448,7 +1564,14 @@ def resolve_database_path(
14481564
)
14491565
fasta_files.append(fasta_file)
14501566
else:
1451-
raise ValueError(f"Unknown database name: {db_name}. Supported: 'RVDB', 'RefSeq', or an accession number (e.g., 'OZ254622.1')")
1567+
# If user passed a species name directly (typically contains spaces), treat it as a species database.
1568+
# This allows: -d "Orthopoxvirus monkeypox" (quotes required for spaces in shells).
1569+
if " " in database_arg and len(db_names) == 1:
1570+
return _resolve_species_database(database_arg)
1571+
raise ValueError(
1572+
f"Unknown database name: {db_name}. Supported: 'RVDB', 'RefSeq', 'RVDB,RefSeq', "
1573+
f"an accession (e.g., 'OZ254622.1'), a FASTA path, or a species name (e.g., 'Orthopoxvirus monkeypox')."
1574+
)
14521575

14531576
# If accessions provided, merge with each database directly in the database directory
14541577
if accession_fasta_files:
@@ -2735,9 +2858,9 @@ def safe_stem(value: str, max_len: int = 80) -> str:
27352858
"""
27362859
# Remove or replace problematic characters
27372860
cleaned = re.sub(r'[^\w\-_\.]', '_', value)
2861+
suffix = hashlib.md5(cleaned.encode()).hexdigest()[:10]
27382862
# Limit length
27392863
if len(cleaned) > max_len:
2740-
suffix = hashlib.md5(cleaned.encode()).hexdigest()[:10]
27412864
cleaned = cleaned[:max_len-10]
27422865
return f"{cleaned}_{suffix}"
27432866

@@ -7303,7 +7426,7 @@ def main(args=None):
73037426
default="RVDB",
73047427
dest="database",
73057428
metavar="",
7306-
help="RVDB|RefSeq|RVDB,RefSeq|accession|FASTA (default RVDB).",
7429+
help="RVDB|RefSeq|RVDB,RefSeq|accession|FASTA|species name (default RVDB).",
73077430
)
73087431
choose_db.add_argument(
73097432
"--rvdb-version",

0 commit comments

Comments
 (0)