Plant_from_bug

Pipeline for extracting plant barcodes (for 2 at a time) from insect Illumina metagenomes, collapsing unique hits, and identifying candidate host plants via GenBank nt. Though initially intended for plant barcodes and insect genomes, the pipeline works for barcodes and genomes from any group of organisms (see Customization Tips section).

Repository Layout

• plant_barcodes/: provide reference barcode FASTA files.
  • files must match MARKER1*.fasta and MARKER2*.fasta
  • e.g., for matK and rbcL files must match matK*.fasta and rbcL*.fasta
• plant_genes_Nov25/data/: place insect metagenomes (*_R1_R2.fastq.gz)
• run_pipeline.sh: orchestrates the workflow from genomes to hits
• NCBI_search_run.sh: allows for batch remote BLASTN submissions against the GenBank nt database for hit identification
• scripts/prepare_fastas.py: converts compressed FASTQ libraries to FASTA
• scripts/collect_unique_hits.py: filters BLAST output to per-subject unique hits (longest alignment wins) and emits updated TSV/FASTA files
• scripts/blast_nt_hits.sh: helper to BLAST unique hits against GenBank nt
• scripts/slurm_blast_task.sh: Slurm array worker for per-sample BLAST jobs
• scripts/slurm_fastq_to_fasta.sh: Slurm array worker for FASTQ to FASTA conversion

Dependencies

• Python 3.8+
• NCBI BLAST+ (blastn, makeblastdb)
• GNU coreutils (sed, cat, etc.; available on most Linux distros)
• Slurm client utilities (sbatch, squeue; sacct is optional but used when present)

Running the Pipeline

# optional: customize threads/e-value
module load ncbi-blast/2.17.0+       # if not already in your shell startup
export THREADS=32
export EVALUE=1e-3
# optional: Slurm/resource tuning
export FASTQ_BATCH_SIZE=10          # max concurrent FASTQ->FASTA jobs
export GENOME_BATCH_SIZE=10        # max concurrent insect genomes
export BLAST_THREADS=3             # threads per BLAST task
export SLURM_MEM_PER_CPU=10G              # whatever per-core memory you need
# (you can mix in other options, e.g., SLURM_PARTITION, SLURM_TIME, etc.)
# export SLURM_PARTITION=general   # uncomment to target a partition

# Analyses stay inside each dataset directory (e.g., plant_genes_Nov25/results)
bash run_pipeline.sh plant_barcodes plant_genes_Nov25 matK rbcL
# or (matK and rbcL are the default)
bash run_pipeline.sh plant_barcodes plant_genes_Nov25

# For a different batch (e.g., December 2025 genomes)
bash run_pipeline.sh plant_barcodes plant_genes_Dec25 matK rbcL
# or (matK and rbcL are the default)
bash run_pipeline.sh plant_barcodes plant_genes_Dec25

# For a different set of barcodes (e.g., ITS and psbA-trnH)
bash run_pipeline.sh plant_barcodes plant_genes_Dec25 ITS psbA-trnH

Notes:

• For reproducability, save a temp.txt file with the exported values to keep track of your settings.
• Pass the dataset root (plant_genes_Nov25). The script looks for data/ inside that directory; if you instead provide the data/ path directly, it is detected automatically.
• Outputs land in <dataset>/results by default. Provide a third argument (e.g. analysis_run2) to place them in <dataset>/analysis_run2. Absolute paths are honored as-is.
• run_pipeline.sh submits Slurm job arrays for both the FASTQ conversion and BLAST stages; run it from a login/submit node with access to your shared filesystem.
• NCBI_search_run.sh submits remote BLASTN searches against the GenBank nt database for hit identification. Results land in <dataset>/results/nt
• BLASTN jobs can fail due to BLAST engine error: Database memory map file error, or due to connection issues with GenBank. Failed runs are automatically rerun 2 times by default, with a lag of 5 minutes before the next attempt. The lag time and number of attempts can be adjusted.

What happens:

1. *_R1_R2.fastq.gz libraries are streamed into results/fastas/*.fasta through a Slurm array (10 concurrent conversions by default)
2. All MARKER1 and MARKER2 barcode references are concatenated per gene (default matK and rbcL if no barcode names provided)
3. Each FASTA becomes its own makeblastdb target
4. Slurm job arrays process the FASTA manifest ~10 samples at a time (tunable) and run blastn (default e-value 1e-3) for matK/rbcL queries
5. BLAST results stay per-sample (no combined tables) and are deduplicated per gene using the subject ID (sseqid). When multiple hits share a subject, the longest alignment (then higher bitscore/lower e-value) is retained, yielding per-sample unique hit tables/FASTAs.

Key outputs (relative to each dataset directory):

• logs (run_pipeline_YMD_HM.log, NCBI_search_YMD_HM.log) are automatically dated and stored in the main repository folder
• results/blast/: raw BLAST tables per sample (one file per insect genome)
• results/unique/by_sample/<gene>/*_{gene}_unique_hits.(fasta|tsv): per-insect unique hits with duplicates (same sseqid) removed by keeping the longest alignment, preserving which genome each barcode hit originated from. These are the primary inputs for downstream analyses (no repository-level combined unique files are produced).
• Intermediates (results/fastas, results/blastdbs, results/manifests) are deleted automatically at the end of each run to save space—rerunning the pipeline regenerates them.
• Slurm stdout/stderr for each BLAST task lands in results/blast/<gene>/slurm-<job>_<task>.out.
• Slurm stdout/stderr for each FASTQ conversion task lands in results/fastas/slurm_fastq/slurm-<job>_<task>.out
• results/nt/: results of BLASTN against GenBank nt

Slurm scheduling details

• FASTQ_BATCH_SIZE (default 10) controls how many FASTQ conversion jobs run simultaneously. Set to 0 or a negative value to allow unlimited concurrency.
• GENOME_BATCH_SIZE (default 10) caps how many insect genomes run simultaneously. Set to 0 or a negative value to allow unlimited concurrency.
• BLAST_THREADS (or THREADS) controls the per-task blastn -num_threads. Ensure GENOME_BATCH_SIZE * BLAST_THREADS fits within your allocation.
• MAX_RETRIES (deafult 2) sets a number of attempts for failed runs
• RETRY_DELAY (deafult 300s) allocates a time lag between running the next attempt
• Optional environment variables forwarded to sbatch: SLURM_PARTITION, SLURM_TIME, SLURM_MEM_PER_CPU, SLURM_CPUS_PER_TASK, SLURM_QOS, plus any extra flags via SLURM_SBATCH_OPTS (space-separated string).
• SLURM_WAIT_POLL (seconds, default 30) sets how often the pipeline polls squeue to wait for arrays to finish.
• Set TASK_INDEX to a zero-based value if you need to run scripts/slurm_blast_task.sh manually for debugging outside of Slurm.

Identifying Plant Taxa via GenBank nt

Remote BLAST requires an email (and benefits from an NCBI API key).

export NCBI_EMAIL="you@example.com"
# optional
export NCBI_API_KEY="XXXX"
# optional for search in other taxa (default Viridiplantae)
export TAX_FILTER="TAXON_OF_CHOICE[ORGN]"

# The script allows for subsampling the specimens in the analysis folder based on the naming pattern
# bash NCBI_search_run.sh PATTERN-FOR-SOME-PARTS GENE ANALYSIS_FOLDER
# example:
bash NCBI_search_run.sh 00- matK plant_genes_Nov25

Outputs are tab-delimited tables enriched with taxonomy columns (staxids, sscinames, sskingdoms, stitle). Adjust MAX_TARGET_SEQS or TAX_FILTER (defaults to Viridiplantae[ORGN]) through environment variables if needed. Iterate over the per-sample FASTA outputs if you want to BLAST every insect genome.

Customization Tips

• Change FASTQ_SUFFIX if your libraries follow a different naming scheme.
• Override THREADS to match the number of CPU cores available.
• Use scripts/prepare_fastas.py --force ... to re-generate FASTA files.
• Tweak EVALUE from the environment when invoking run_pipeline.sh or scripts/blast_nt_hits.sh.
• GENOME_BATCH_SIZE, BLAST_THREADS, and the SLURM_* environment variables let you shape how aggressively the BLAST stage submits work to your cluster.
• Though initially intended for plant barcodes and insect genomes, the pipeline works for barcodes and genomes from any group of organisms. The only adjustment needed is while using NCBI_search_run.sh, make sure to change Viridiplantae[ORGN] to a database that matches your barcode organisms.