Pathogen Detection Using Machine Learning

Reference implementation for the book chapter "Pathogen Detection Using Machine Learning" (Senivarapu S, Coxe T, Murari A, Chandi A, Azad RK).

The protocol builds a labeled dataset of bacterial genomes from BacSPaD, downloads the corresponding assemblies from BV-BRC, and compares four classification pipelines on the same train/test split:

Pipeline	Features	Model
A	6-mer frequencies	LinearSVC(penalty="l1", dual=False, C=1.0)
B	6-mer frequencies	RandomForestClassifier(n_estimators=500, random_state=0)
C	250 bp DNA fragments	Reverse-complement CNN, 10 epochs
D	6-mer frequencies	GradientBoostingClassifier(n_estimators=200, learning_rate=0.05, max_depth=3, random_state=0)

The example dataset is a balanced subset of 100 HP (pathogenic to humans) and 100 NHP (non-pathogenic to humans) genomes.

Requirements

• Python 3.10, 16 GB RAM, ~50 GB free disk space
• A stable internet connection for downloading genome assemblies

Setup

git clone https://github.com/sudeeps1/pathogenprediction.git
cd pathogenprediction
conda create -n pathogen python=3.10 -y
conda activate pathogen
pip install -r requirements.txt

Repository layout

pathogenprediction/
├── labels/
│   ├── Genomes_labeled.csv           # BacSPaD table export (input to build_labels.py)
│   ├── accessions_pathogenic.txt     # 100 HP genome identifiers
│   ├── accessions_nonpathogenic.txt  # 100 NHP genome identifiers
│   ├── genome_list.txt               # all 200 identifiers
│   └── labels.csv                    # id,label (1 = HP, 0 = NHP)
├── data/genomes/                     # genome FASTA files, downloaded in step 2
├── results/                          # split, predictions and metrics
├── figures/                          # chapter figures
└── scripts/
    ├── build_labels.py               # Section 3.2
    ├── make_split.py                 # Section 3.4
    ├── kmer_table.py                 # Section 4.1
    ├── train_svm.py                  # Pipeline A
    ├── train_rf.py                   # Pipeline B
    ├── make_fragments.py             # Section 4.3
    ├── train_rc_cnn.py               # Pipeline C
    ├── train_gb.py                   # Pipeline D
    └── score.py                      # Section 5

Running the protocol

1. Build the labeled dataset (Section 3.2)

labels/Genomes_labeled.csv is the BacSPaD table export and is included in this repository, so this step runs without any manual download. To rebuild it from source, export the table from https://bacspad.altrabio.com/#data and save it to that path.

python scripts/build_labels.py

This samples 100 HP and 100 NHP genome identifiers with a fixed random seed and writes labels/accessions_pathogenic.txt, labels/accessions_nonpathogenic.txt, labels/genome_list.txt and labels/labels.csv. The generated files are already committed, so the download step below can be run directly after cloning.

2. Download the genome assemblies (Section 3.3)

BV-BRC requires an SSL-enabled FTP client, so pass --ssl-reqd to curl and use the ftp.bv-brc.org host name (the TLS certificate does not cover the ftp.bvbrc.org alias):

mkdir -p data/genomes
for i in $(cat labels/genome_list.txt); do
  curl -s --ssl-reqd -f -o "data/genomes/$i.fna" "ftp://ftp.bv-brc.org/genomes/$i/$i.fna"
done

Verify that 200 files were downloaded (see Note 1 and Note 2):

ls data/genomes | wc -l

3. Create the train/test split (Section 3.4)

python scripts/make_split.py
awk -F, 'NR>1 {n[$3"_"$2]++} END {for (k in n) print k, n[k]}' results/split.csv

Expected: TRAIN_1 80, TRAIN_0 80, TEST_1 20, TEST_0 20.

4. Run the four pipelines (Section 4)

python scripts/kmer_table.py       # results/k6_features.csv
python scripts/train_svm.py        # results/svm_predictions.csv
python scripts/train_rf.py         # results/rf_predictions.csv
python scripts/train_gb.py         # results/gb_predictions.csv
python scripts/make_fragments.py   # results/fragments.tsv  (~850 MB, see Note 8)
python scripts/train_rc_cnn.py     # results/cnn_predictions.csv

train_rc_cnn.py uses CUDA or Apple MPS when available and falls back to the CPU otherwise (see Note 6). To shrink the fragment file and shorten CNN training, set MAX_FRAGMENTS_PER_GENOME in scripts/make_fragments.py to a positive value.

5. Evaluate (Section 5)

python scripts/score.py

This writes results/metrics.txt. A copy of the expected output is committed to the repository so a run can be checked against it.

Reproducibility

results/split.csv, the four prediction files and results/metrics.txt are committed. The large intermediates results/k6_features.csv and results/fragments.tsv are regenerated by the scripts and are not tracked.

License

MIT — see LICENSE.