PlasmidScreen detects engineered DNA in sequencing reads using a Kraken2 minimizer/kmer scan (taxid 32630 in sliding windows). Optional codon adaptation index (CAI) scoring uses DIAMOND blastx for ORF coordinates and host taxonomy, compared against host codon usage tables from the Codon Statistics Database.
Design specification (inputs/outputs): docs/DESIGN.md
| Step | Tool | Purpose |
|---|---|---|
| 1 | CSDB build | codon_tables.json + optional taxonomy_parents.json |
| 2 | Kraken2 | Engineered k-mer labels (Natural / Synthetic) |
| 3 | DIAMOND blastx | ORF intervals (qstart/qend) and host taxid (staxids) |
| 4 | CAI | Codon adaptation vs host reference (Natural reads only in full screen) |
| 5 | Combined call | ScreenResult.per_read[].engineered_overall / overall_label |
Each read in screen_result.per_read gets an overall_label (Natural or Synthetic) from the thresholds you
pass to run_screen which can be a result from the two components K-mer scan, and Codon CAI which each are optionally
ran. If either are found to contain engineering, ScreenResult object will report engineering.
| Signal | Parameter | Counts toward overall engineered when |
|---|---|---|
| K-mer scan | --threshold, --window_size |
Engineered minimizers (taxid 32630) in a window ≥ threshold → k-mer Synthetic |
| Codon CAI | --codon-cai-threshold (optional) |
CAI < threshold on a scored Natural read (requires codon step) |
- Without
--codon-cai-threshold, overall matches the k-mer scan only (engineered_scan.synthetic_count==overall_synthetic_count). - With
--codon-cai-threshold, a read can be k-mer Natural but overall Synthetic if CAI is low. - Use
screen_result.overall_synthetic_count,engineered_read_ids, andnatural_read_ids_overallfor run-level summaries (not onlyengineered_scan.synthetic_count).
from plasmidScreen import (
build_codon_database,
run_screen,
)
# Step 1 — build reference JSON from Codon Statistics Database,
# please specify path otherwise download will occur in ~/.local/share/PlasmidScreen/
build_codon_database(output_dir="/path/to/codon_usage_db")
# Full pipeline: Kraken engineered scan + codon optimization detection on Natural reads
screen_result = run_screen(
"reads.fa",
kraken_db="/path/to/kraken/db", # Needed for Kraken2
diamond_db="/path/to/protein.dmnd",
codon_usage_dir="/path/to/codon_usage_db",
engineered_report_path="engineered_report.txt", # omit for in-memory results only
codon_cai_engineered_threshold=0.6,
)
# Here is the found engineered reads
print(screen_result.engineered_read_ids)
# Here is the engineered reads based on kmer scanning
print(screen_result.engineered_scan.engineered_read_ids)
# Here is the engineered reads from codon optimization detection
print(screen_result.codon_adaptation.engineered_read_ids)
# Here is the natural read ids found:
print(screen_result.natural_read_ids_overall)
# If you want to get more granular information of each read's call you can look at the list of reads
# This contains the ReadFlagDetail object that holds additonal info
for r in screen_result.per_read:
print(
r.read_id,
r.overall_label, # Returns literal of natural or synthetic
r.engineered_overall, # boolean True if engineered
r.engineered_methods, # returns list of strings of type of engineering detection performed on the read
r.engineered_by_kmer_scan, # boolean if engineered by kmer scan
r.cai_vs_host, # boolean if engineered by codon optimization threshold
)# Default: import all taxids in the CSDB archive (~5.2 GB download on first run)
python plasmidScreen.py build
# Use an existing CSDB archive (no download)
python plasmidScreen.py build \
--csdb-archive /data/codonstatsdb_March2022.tar.gz \
--no-download-csdb
The archive is cached under ~/.local/share/PlasmidScreen/ unless --csdb-archive is set.
Output: codon_usage/codon_tables.json and taxonomy_parents.json.
pip install -r requirements.txt
We recommend custom databases for both kraken2 and diamond. The suggested kraken2 database contains a curated set of sequences that is tailored for finding the most amount of signal for detecting engineering. In order to best indentify codon optimization, identifying the host taxonomy of the found ORF is crucial. Therefore, a comprehensive diamond database is recommended. The suggested diamond database is composed of all UniRef100 sequences with TAXIDs assigned.
- K-mer scan: Custom kraken2 database (~20GB)
- Codon Optimization: UniRef100 Diamond Database (~235GB)
python plasmidScreen.py screen reads.fa report.txt /path/to/kraken/db \
--diamond-db /path/to/protein.dmnd \
--codon-usage-dir ~/.local/share/PlasmidScreen/codon_usage \
--codon-usage-output codon_usage.tsv
# Reuse Kraken classifications (no Kraken subprocess)
python plasmidScreen.py screen reads.fa report.txt /path/to/kraken/db \
--no-run-kraken --kraken-output-path kraken.out \
--diamond-db /path/to/protein.dmnd
# Save DIAMOND TSV for debugging / codon-only reruns
python plasmidScreen.py screen reads.fa report.txt /path/to/kraken/db \
--diamond-db /path/to/protein.dmnd \
--debug-write-diamond-out --diamond-output-path diamond.tsv
# Re-run codon step from saved DIAMOND output
python plasmidScreen.py screen reads.fa report.txt /path/to/kraken/db \
--no-run-diamond --diamond-output-path diamond.tsv \
--no-run-kraken --kraken-output-path kraken.out
# Low CAI can flag overall engineered (on Natural reads with CAI scores)
python plasmidScreen.py screen reads.fa report.txt /path/to/kraken/db \
--diamond-db /path/to/protein.dmnd \
--codon-cai-threshold 0.6 \
--threshold 10
# Engineered k-mer scan only (no codon / DIAMOND; overall == k-mer labels)
python plasmidScreen.py screen reads.fa report.txt /path/to/kraken/db --skip-codon-usageIf a host taxid has no resolvable codon table, screening raises MissingCodonReferenceError — build or extend the reference first (no online fetch at runtime).
The output contains two columns representing the confidence that a given read is of engineered origin. The first column reports the percentage of k-mers classified as engineered out of all k-mers in the read. The second column provides a p-value estimating the probability of observing this fraction of engineered k-mers by chance.
The p-value is calculated using a binomial distribution based on the total number of classified k-mers, the number of classified engineered k-mers, and a background rate of engineered k-mers in non-engineered sequences. This background rate represents the probability that a random, naturally occurring k-mer is misclassified as engineered. We recommend a conservative default background rate of 0.1%; lower values result in more aggressive filtering. This rate was established on a curated dataset of non-engineered sequences, comprising bacterial chromosomes (>22,000) and naturally occurring incomplete plasmids from IMG-PR (>561,000), each chopped into 1000 bp fragments.
Reported p-values are adjusted for multiple testing across all classified reads using the Benjamini-Hochberg False Discovery Rate (FDR) correction.