Metadata-Version: 2.4
Name: plasrisk
Version: 1.6.0
Summary: PlasRisk: data-driven weighted risk assessment for bacterial plasmids from FASTA sequences (raw WGS/long-read/metagenomic assembly prediction front end, isolate and metagenomic modes; lite 5-dim FASTA-only core by default, full 10-dim model, and original PIPdb ordinal index)
Author: PlasRisk Team
Maintainer: PlasRisk Team
License: MIT
Project-URL: Homepage, https://github.com/LLQ95/PlasRisk
Project-URL: Repository, https://github.com/LLQ95/PlasRisk
Project-URL: Issues, https://github.com/LLQ95/PlasRisk/issues
Project-URL: Documentation, https://github.com/LLQ95/PlasRisk#readme
Project-URL: Changelog, https://github.com/LLQ95/PlasRisk/releases
Keywords: bioinformatics,plasmid,antimicrobial resistance,AMR,risk assessment,microbial genomics,One Health,biocide resistance,co-selection,horizontal gene transfer
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.8
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Operating System :: OS Independent
Classifier: Environment :: Console
Requires-Python: >=3.8
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: pandas>=1.3
Requires-Dist: numpy>=1.20
Provides-Extra: dev
Requires-Dist: pytest>=7.0; extra == "dev"
Requires-Dist: pytest-cov; extra == "dev"
Requires-Dist: build; extra == "dev"
Requires-Dist: twine; extra == "dev"
Provides-Extra: full
Requires-Dist: biopython>=1.79; extra == "full"
Dynamic: license-file

# PlasRisk

**Data-driven weighted risk assessment for bacterial plasmids from FASTA sequences**
**(three modes: 5-dimension FASTA-only lite core by default, 10-dimension full model, and the original PIPdb ordinal index)**

[![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](LICENSE)
[![Python](https://img.shields.io/badge/python-3.9+-blue.svg)](https://python.org)

![PlasRisk graphical abstract](docs/graphical_abstract.svg)

PlasRisk computes a composite risk score for bacterial plasmids based on
biologically motivated dimensions: antimicrobial resistance gene (ARG) burden,
virulence factors (VFs), mobility/conjugation potential, host range, replicon
type, plasmid size, biocide/metal resistance (BMRG), geographic spread, habitat
breadth, and temporal growth rate. It accepts plasmid FASTA sequences and
automatically annotates them using [abricate](https://github.com/tseemann/abricate).

The default output is the FASTA-only lite core:

```
S_lite = 0.253*S_ARG + 0.237*S_BMG + 0.201*S_MOB + 0.175*S_SIZE + 0.135*S_VF
```

The full 10-dimension model (`--mode full`) is:

```
S_full = 0.237*S_ARG + 0.222*S_BMG + 0.189*S_MOB + 0.164*S_SIZE
       + 0.126*S_VF + 0.030*S_HOST + 0.019*S_HAB + 0.005*S_GEO
       + 0.004*S_REP + 0.002*S_GROW

Weights derived by data-driven consensus (Random Forest MDG, LASSO, and
grid-search optimization) on 792,964 PIPdb PSCs; sum ≈ 1.0.
```

The original PIPdb ordinal index (`--mode ordinal`) reproduces the published
8-item formula (Zhu et al., 2025, Table 1):

```
CRI = clamp(Round((pathogenic_phylum + pathogenic_species + habitats
      + ARGs + VFGs + 2*WHO_ARGs + ISs + growth_rate)/8 + 0.6), 1, 5)
```

Each item is binned to an ordinal value of 1–5; the index runs from 1
(Minimal) to 5 (Very High), with a normalized value of (CRI−1)/4.

Risk grades: **A** (Very High, S >= 0.60), **B** (High, >= 0.45),
**C** (Moderate, >= 0.30), **D** (Low, >= 0.15), **E** (Minimal, < 0.15).

### AWaRe category table for S_ARG

S_ARG weights genes by the WHO AWaRe categories (Access, Watch, Reserve).
The frozen component matrix used the legacy keyword mapping, which remains
the default so published scores are exactly reproducible. The WHO AWaRe
2025 table can be selected instead:

```bash
plasrisk --aware-version 2025 plasmid.fasta
```

```python
scorer = PlasRiskScorer(mode="lite", aware_version="2025")
```

Under the 2025 mapping, biocide/metal genes contribute 0 to S_ARG because
they are scored under S_BMG, and the separate 1.5x last-resort multiplier
is not applied (Reserve already carries the highest weight). An audit of
all 792,964 PSCs found that the 2025 mapping changes S_ARG by at least
0.05 for 4.9% of PSCs and moves the high-risk-ARG AUC on the locked test
from 0.9807 to 0.9789 (delta -0.0018); the frozen values are retained
(see Method S14 and Table S56 in the Supplementary Materials).

---

## Installation

### Option 1: conda (recommended)

```bash
# Create a dedicated environment
conda create -n plasrisk -c bioconda -c conda-forge plasrisk
conda activate plasrisk
```

This installs PlasRisk together with `abricate` and `blast` for full
annotation capability.

### Option 2: pip + manual abricate

```bash
pip install plasrisk

# Install abricate separately for annotation
conda install -c bioconda abricate
# or on Debian/Ubuntu: apt-get install abricate
```

### Option 3: from source

```bash
git clone https://github.com/LLQ95/PlasRisk.git
cd plasrisk
pip install .

# Install annotation dependencies
conda install -c bioconda abricate blast
```

### Set up abricate databases

After installing abricate, download the databases you need:

```bash
# Download/update all default databases
abricate-get_db --db card --force
abricate-get_db --db vfdb --force
abricate-get_db --db plasmidfinder --force
abricate-get_db --db resfinder --force
abricate-get_db --db ncbi --force

# Verify
abricate --list
```

For BacMet (biocide/metal resistance) database, see:
https://github.com/tseemann/abricate#making-your-own-database

---

## Quick start

```bash
# Score a single plasmid
plasrisk plasmid.fasta

# Score multiple plasmids
plasrisk *.fasta

# Score all FASTA files in a directory
plasrisk /path/to/plasmids/

# Specify output directory
plasrisk -o results *.fasta

# Use specific abricate databases
plasrisk --db card,vfdb,plasmidfinder,bacmet plasmid.fasta

# Sequence-only mode (no abricate needed; scores based on length + replicon lookup)
plasrisk --no-abricate contigs.fasta

# JSON output
plasrisk --json -o results plasmid.fasta

# Full 10-dimension model with epidemiological context (default is lite)
plasrisk --mode full *.fasta

# Original PIPdb 8-item ordinal index (grades 1-5)
plasrisk --mode ordinal *.fasta

# Verify that the package S_SIZE transform and the discovery-pipeline
# transform rank plasmids identically (writes rank_concordance.tsv)
plasrisk --rank-concordance *.fasta
```

### Lite (default), Full and Ordinal modes

The default output is the 5-dimension FASTA-only lite core; no PIPdb-style
metadata are required. The full 10-dimension model is available with
`--mode full`, and the original PIPdb 8-item ordinal index with
`--mode ordinal`.

| | Lite (default) | Full | Ordinal |
|---|---|---|---|
| Flag | `--mode lite` | `--mode full` | `--mode ordinal` |
| Output scale | continuous 0–1, grades A–E | continuous 0–1, grades A–E | ordinal index 1–5 |
| Components | S_ARG, S_BMG, S_MOB, S_SIZE, S_VF | lite 5 plus S_HOST, S_HAB, S_GEO, S_REP, S_GROW | pathogenic phylum/species, habitats, ARGs, VFGs, WHO ARGs (2×), ISs, growth rate |
| Weights / bins | 0.253, 0.237, 0.201, 0.175, 0.135 | 0.237, 0.222, 0.189, 0.164, 0.126, 0.030, 0.019, 0.005, 0.004, 0.002 | published PIPdb ordinal bins, 1–5 per item |
| Mean AUC (4 outcomes) | 0.919 | 0.919 | CRI AUC 0.81–0.83 (internal consistency) |
| Grade agreement | 92.7% exact / 100% within one grade vs. full | reference | not applicable |
| Required annotations | ARG + VF + mobility + length + BacMet | ARG + VF + mobility + replicon + BacMet + metadata | ARG + VF + WHO ARG + IS counts; unknown metadata default to bin 1 |
| Use case | Rapid FASTA-only screening, resource-limited settings | Comprehensive risk assessment with epidemiological context | Reproducing the published PIPdb index for direct comparison |

### S_SIZE normalization and rank-concordance check

The released package computes the length component with a monotonic sigmoid
on log10 length, centered at 30 kb:

```
S_SIZE = 1 / (1 + exp(-4 * (log10(L) - log10(30,000))))
```

The discovery pipeline used the log-ratio `clip(log10(L)/log10(1,000,000), 0, 1)`
for weight fitting. Both transforms are strictly increasing in L, so plasmid
rankings are identical (Spearman rho = 1.000, Kendall tau = 1.000 over the
792,964 PIPdb PSCs); only the absolute scaling differs. The equivalence can be
reproduced on any input set with `plasrisk --rank-concordance` (writes
`rank_concordance.tsv`) or via the Python API
(`from plasrisk import size_rank_concordance`). See Supplementary Text S1.

### Example output

```
  Grade A:     3 (  3.2%) ######
  Grade B:    12 ( 12.9%) ##########################
  Grade C:    28 ( 30.1%) ############################################################
  Grade D:    35 ( 37.6%) ############################################################################
  Grade E:    15 ( 16.1%) ################################

  Top 10 highest-risk plasmids:
    pNDM-1_260kb                            S=0.712  grade A  8 ARG  IncX3 [blaNDM]
    pMCR-1_33kb                             S=0.581  grade B  4 ARG  IncX4 [mcr]
    pKPC-2_110kb                            S=0.534  grade B  6 ARG  IncFII(K) [blaKPC]
```

### Output files

| File | Description |
|------|-------------|
| `plasrisk_results.tsv` | Per-sequence scores: active components, S_total, S_norm, grade, gene lists |
| `rank_concordance.tsv` | S_SIZE rank-concordance check (with `--rank-concordance`) |
| `plasrisk_summary.tsv` | Per-file summary: counts, grade distribution, mean/max scores |
| `plasrisk_results.json` | JSON format (with `--json`) |

---

## Metagenomic mode (`plasrisk meta`)

PlasRisk extends to fragmented metagenomic assemblies with three aggregation
layers. The frozen consensus weights are never refit; every output carries
an explicit completeness/confidence label so that fragments are not scored
as if they were closed plasmids.

- **Layer 1, contig/element level.** S_ARG, S_BMG and S_VF are always
  evaluated from observed genes. S_MOB is evaluated only from observed
  transfer apparatus (MPF/T4CP, relaxase, oriT) or insertion-sequence
  evidence; isolate-derived replicon/global priors are never injected when
  the apparatus was simply not covered. S_SIZE is evaluated only for
  circularised headers (`circular=true`) or contigs at/above
  `--complete-length` (default 10,000 bp); fragment lengths are never
  passed through S_SIZE, which would penalise megaplasmids that short-read
  assembly split. Renormalised scores are reported as S_cargo3 (cargo
  dimensions), S_frag4 (cargo + mobility) and S_lite5 (all five lite
  dimensions); A-E grades are emitted only when all five dimensions are
  assessable.
- **Layer 2, plasmid-bin level.** Contigs assigned to the same predicted
  plasmid (geNomad components, SCAPP plasmids, MOB-suite bins or any
  generic contig-to-bin table) are aggregated: gene sets are unioned with
  deduplication, lengths are summed, mobility markers are OR-combined, and
  per-contig coverage is length-weighted. Bins are scored with the standard
  lite model and labelled complete / probable / fragmented.
- **Layer 3, sample/community level.** Descriptive abundance-aware
  statistics: coverage-length weighted mean score (CPRI_mean), P90/P99,
  grade prevalence, and coverage-weighted grade A/B and critical-ARG
  loads (with optional per-Mb normalisation). CPRI is a descriptive
  extension that was not trained or validated on metagenomic data; this is
  stated in every sample summary.

```bash
# (1) Predict plasmids with geNomad (v1.10+, aggregate mode)
genomad aggregate --clean-up output_dir/ assembly.fasta genomad_db/

# (2) Score all three layers; coverage is optional but recommended
plasrisk meta assembly.fasta \
  --genomad-bin output_dir/ \
  --coverage jgi_depth.tsv \
  --min-plasmid-score 0.7 \
  -o plasrisk_meta/

# Generic bin table instead of geNomad (works with SCAPP/MOB-suite output):
plasrisk meta assembly.fasta --bins contig_to_plasmid.tsv \
  --coverage depth.tsv -o plasrisk_meta/

# Contig-only scoring without binning (each contig is a singleton element):
plasrisk meta assembly.fasta -o plasrisk_meta/
```

The generic bins table is a TSV/CSV with `seq_id,bin_id` columns (optional
`circular`/`topology` and `plasmid_score`). For SCAPP, one bin is one
plasmid FASTA in `plasmid_fastas/`; build the table from the record headers
and pass it with `--bins`. Coverage tables accept the
jgi_summarize_bam_contig_depths format; use `--coverage-column` for
multi-sample depth files.

Outputs: `meta_contig_scores.tsv` (Layer 1), `meta_bin_scores.tsv`
(Layer 2), and `meta_sample_summary.tsv` (Layer 3, plus JSON with
`--json`).

### Fragmentation recovery simulation

Before applying PlasRisk to a new assembly type, quantify how fragmentation
affects scores with a controlled grid on closed reference plasmids:

```bash
plasrisk-fragment-sim closed_plasmids.fasta \
  --frag-lengths 500,1000,2000,5000,10000,20000 \
  --replicates 10 -o fragment_sim/
```

For each fragment length it reports Spearman rank correlation of
oracle-binned scores against intact-plasmid scores, grade exact/within-one
agreement, median score bias, and critical-ARG recovery. This is the
upper-bound (oracle binning) layer; pair it with ART/wgsim reads,
metaSPAdes/MEGAHIT assembly and geNomad binning for the realistic layer.

---

## Python API

```python
from plasrisk import PlasRiskScorer, PlasmidFeatures, annotate_fasta, load_replicon_lookup

# Option A: annotate a FASTA file directly (default = 5-dim FASTA-only lite)
lookup = load_replicon_lookup()
result = annotate_fasta("plasmid.fasta", lookup=lookup)
scorer = PlasRiskScorer(replicon_lookup=lookup)  # lite by default
df = scorer.score_dataframe(result.features)
print(df[["seq_id", "S_norm", "grade", "high_risk_genes"]])

# Full 10-dimension model with epidemiological context
scorer_full = PlasRiskScorer(replicon_lookup=lookup, mode="full")
df_full = scorer_full.score_dataframe(result.features)

# Original PIPdb 8-item ordinal index (three ways to obtain it)
from plasrisk import PIPdbScorer, get_scorer
scorer_ord = PIPdbScorer()                       # direct class
scorer_ord = get_scorer("pipdb")                 # factory, model name
scorer_ord = get_scorer("plasrisk", mode="ordinal")  # factory, mode
df_ord = scorer_ord.score_dataframe(result.features)

# S_SIZE rank-concordance check (package sigmoid vs pipeline log-ratio)
from plasrisk import size_rank_concordance
print(size_rank_concordance([f.length_bp for f in result.features]))

# Option B: construct features manually
feat = PlasmidFeatures(
    seq_id="pExample",
    length_bp=85000,
    arg_names=["NDM-1", "CTX-M-15", "TEM-1"],
    vf_names=["aerobactin"],
    vf_categories=["Nutritional/Metabolic factor"],
    bm_gene_names=["merA", "qacEdelta1"],
    replicon="IncX3",
    has_t4cp=True,
    has_relaxase=True,
    has_oriT=True,
    has_auxiliary=True,
)
scores = scorer.score(feat)
print(f"S_norm = {scores['S_norm']:.3f}, grade = {scores['grade']}")
```

---

## The 10 risk dimensions

| Component | Weight | What it measures | Scoring basis |
|-----------|--------|------------------|---------------|
| **S_ARG** | 0.237 | ARG count, WHO-priority genes, high-risk genes (mcr, NDM, KPC, CTX-M, tetX, etc.) | Base + per-gene + high-risk bonuses |
| **S_BMG** | 0.222 | Biocide/metal resistance (mer, qac, ars/cop/sil), co-selection potential | Base + per-gene + family bonuses |
| **S_MOB** | 0.189 | T4CP, relaxase, oriT, auxiliary transfer proteins | Element-based additive score |
| **S_SIZE** | 0.164 | Plasmid length (cargo capacity) | Sigmoid: midpoint 30 kb |
| **S_VF** | 0.126 | VF count, exotoxins, secretion systems (T3SS/T4SS) | Base + per-gene + category bonuses |
| **S_HOST** | 0.030 | Number of host genera / replicon prior | Empirical host range or lookup |
| **S_HAB** | 0.019 | Habitat breadth (human/animal/environment) | PIPdb-derived lookup |
| **S_GEO** | 0.005 | Number of countries observed | PIPdb-derived lookup |
| **S_REP** | 0.004 | Replicon backbone risk (IncX3, IncN, ColKP3 high; ColpVC low) | PIPdb-derived lookup table |
| **S_GROW** | 0.002 | Annual growth rate of the replicon | PIPdb-derived lookup |

---

## Command-line options

```
plasrisk [options] <fasta1> [fasta2 ...]

positional arguments:
  FASTA                 FASTA file(s) or directory

options:
  -o, --output DIR      Output directory (default: ./plasrisk_output)
  -t, --threads N       Number of abricate threads (default: 4)
  --min-id FLOAT        Minimum abricate identity % (default: 75)
  --min-cov FLOAT       Minimum abricate coverage % (default: 50)
  --no-abricate         Skip abricate; sequence-only scoring
  --mode {lite,full,ordinal}
                        Scoring mode: lite (5-dim FASTA-only core, default),
                        full (10-dim weighted), or ordinal (PIPdb 8-item index)
  --model {plasrisk,pipdb}
                        Deprecated alias; pipdb is equivalent to --mode ordinal
  --rank-concordance    Verify S_SIZE rank concordance (sigmoid vs log-ratio)
  --db LIST             Comma-separated abricate databases (default: auto)
  --json                Also write JSON output
  -q, --quiet           Suppress progress messages
  -v, --version         Show version
  -h, --help            Show help
```

---

## Model validation

The PlasRisk model was developed and validated using 792,964 plasmid sequence
clusters from PIPdb (Zhu et al., *Nucleic Acids Res.*, 2025). Validation
included:

- **Quartile stratification**: Q1 (highest risk) plasmids had 65.6% ARG prevalence,
  12.2% high-risk ARG rate, 13.3% conjugative rate (vs. 0% in Q4).
- **Data-driven weights**: RF-MDG, LASSO, and grid-search optimization across four
  outcomes (high-risk ARG, MDR-VF fusion, conjugative capacity, BMRG carriage)
  converged on S_ARG (0.237), S_BMG (0.222), S_MOB (0.189), and S_SIZE (0.164)
  as dominant predictors.
- **AUC validation**: Final weights achieved AUC 0.955 (high-risk ARG), 0.967
  (MDR-VF fusion), 0.846 (conjugation), 0.910 (BMRG); mean 0.919.
- **Dimensionality analysis**: All-subsets evaluation (1,023 subsets, 5-fold CV)
  showed that 5 core dimensions (ARG+BM+MOB+SIZE+VF) achieve equivalent mean AUC
  to the full 10-dim model (0.919 vs. 0.919), provided as `--mode lite`.
  Overfitting diagnostics (train-test gap, bootstrap optimism, learning curves)
  confirmed no excess optimism in either model.
- **Leave-one-replicon-out CV**: mean AUC = 0.962 across 40 replicons.
- **Weight perturbation sensitivity** (100 iterations, +/-30%): mean Spearman
  rho = 0.994, mean top-10 overlap = 9.2/10.
- **External case-control challenge**: 367 independent plasmids (167 carrying
  critical ARGs and 200 without ARG/VF), deduplicated at 95% ANI over 80%
  coverage, were scored with labels derived only from abricate annotation;
  AUC 0.998 (95% CI 0.994-1.000).

---

## Uploading to conda (bioconda)

To make PlasRisk installable via `conda install -c bioconda plasrisk`:

### Step 1: Upload to PyPI

```bash
# Install build tools
pip install build twine

# Build distributions
python -m build

# Upload to PyPI
twine upload dist/*
```

### Step 2: Fork and clone bioconda-recipes

```bash
git clone https://github.com/bioconda/bioconda-recipes.git
cd bioconda-recipes
```

### Step 3: Create the recipe

```bash
# Create recipe directory
mkdir -p recipes/plasrisk
```

Create `recipes/plasrisk/meta.yaml`:

```yaml
{% set version = "1.0.0" %}

package:
  name: plasrisk
  version: {{ version }}

source:
  url: https://pypi.io/packages/source/p/plasrisk/plasrisk-{{ version }}.tar.gz
  sha256: <SHA256 from PyPI>

build:
  number: 0
  noarch: python
  entry_points:
    - plasrisk = plasrisk.cli:main
  script: "{{ PYTHON }} -m pip install . --no-deps --ignore-installed -vv"

requirements:
  host:
    - python >=3.9
    - pip
    - setuptools >=61.0
    - wheel
  run:
    - python >=3.9
    - pandas >=1.3
    - numpy >=1.20
    # abricate/blast optional; CLI falls back to --no-abricate mode

test:
  imports:
    - plasrisk
  commands:
    - plasrisk --help
    - plasrisk --version

about:
  home: https://github.com/LLQ95/PlasRisk
  license: MIT
  license_file: LICENSE
  summary: "Ten-dimension data-driven weighted risk assessment for bacterial plasmids"
```

> **Note:** The complete, ready-to-submit recipe is in the `bioconda/` directory.
> See `UPLOAD_GUIDE.md` for the full step-by-step release process.

### Step 4: Test locally

```bash
# Install bioconda-utils
conda install -c bioconda bioconda-utils

# Test the recipe
bioconda-utils build recipes/plasrisk --docker
```

### Step 5: Submit a pull request

```bash
git checkout -b plasrisk
git add recipes/plasrisk/
git commit -m "Add plasrisk recipe"
git push origin plasrisk
# Open PR at https://github.com/bioconda/bioconda-recipes
```

Once the PR is merged and CI passes, PlasRisk will be installable via:

```bash
conda install -c bioconda plasrisk
```

### Local conda build (without bioconda)

```bash
# Build from the conda/ directory in this repo
conda build conda/

# Install locally
conda install --use-local plasrisk
```

---

## Running tests

```bash
cd plasrisk_py
python -m pytest tests/ -v
# or
python tests/test_scoring.py
```

---

## Citation

If you use PlasRisk, please cite:

> [Authors]. PlasRisk: a ten-dimension weighted risk assessment framework for
> bacterial plasmids. *Journal*, 2025. doi: [to be added]

The model is based on data from:
> Zhu Q, Chen Q, Lu X, et al. PIPdb: a comprehensive plasmid sequence resource
> for tracking the horizontal transfer of pathogenic factors and antimicrobial
> resistance genes. *Nucleic Acids Research*, 2025, 53(D1):D169-D178.
> doi:10.1093/nar/gkae952

---

## License

MIT License - see [LICENSE](LICENSE) for details.
