its-fun-2-map

A pipeline for processing fungal genome skims from museum specimens. Includes quality control, reference mapping (as a read-baiting step), assembly, BLAST taxonomic assignment, primer-based ITS region extraction, and validation.

Contents

1. Installation and quick start
2. Workflow
3. Configuration
4. Pipeline Overview
5. Output Directory Structure
6. Dependencies
7. Authors

Installation and quick start

Everything the pipeline needs (tools, R, Python and their packages) are installed by the conda environment in its-fun-2-map.yaml. The only things you need to fetch yourself are the reference databases.

1. Install conda if you don't already have it.

2. Clone the repository and create the environment:

git clone https://github.com/museomics/its-fun-2-map.git
cd its-fun-2-map
conda env create -f its-fun-2-map.yaml
conda activate its-fun-2-map

3. Install the pipeline package into the activated environment:

pip install -e .

(or pip install its-fun-2-map to install from PyPI). The conda environment is still required either way, as the pipeline shells out to fastp, seqkit, bwa, samtools, BBMap, SPAdes, BLAST+ and R, none of which are pip-installable.

4. Download the reference databases from UNITE:

• UNITE Full "UNITE+INSD" - reference sequence retrieval (Step 2)
• UNITE General Release (sh_general_release_dynamic) - BLAST Round 1 (Step 5)
• UNITE UCHIME ITS2 reference dataset - BLAST Round 2a (Step 6)
• UNITE UCHIME ITS1 reference dataset - BLAST Round 2b (Step 7)

5. Check the install:

itsfun-qc --help

6. Prepare a tracking sheet (CSV or XLSX) with sample IDs and taxids - see Tracking Sheet Format.

7. Run the pipeline. All steps can be run sequentially with run_pipeline.sh: set your paths and parameters in its configurable parameters block, then run it. See the tutorial for a worked example.

Full version details for every dependency are in Dependencies.

Workflow

Configuration

Default Parameters

Tracking Sheet Format

The pipeline accepts CSV or XLSX tracking sheets with the following columns:

Required columns:

• Sample ID column (configurable via --column_name or --id_header)
• Taxid column (for UNITEd step, configurable via --taxid_header)

Optional but recommended columns:

• forward / fwd - path to forward reads
• reverse / rev - path to reverse reads

If forward/reverse columns are present, file paths are read directly from the tracking sheet. Otherwise, files are searched for in the input directory using sample ID-based string matching.

For taxonomy validation:

• Kingdom, Phylum, Class, Order, Family, Genus, Species - taxonomic hierarchy columns

If your spreadsheet doesn't contain taxonomic hierarchy columns, you can use the helper command itsfun-lineage. This will take any spreadsheet (both CSV and XLSX allowed) and queries all taxonomic IDs within in using NCBI Entrez. It then appends the taxonomic hierarchy columns to the spreadsheet. Outputs a log file and a CSV to your working directory for downstream use.

If using the pipeline shell script, this can be opted in by changing "need_lineage" to "YES" in Advanced Usage: Configurable paramaters.

## Configurable parameters
# itsfun-refs parameters
need_lineage="NO"     #change to "YES" to use itsfun-lineage

Pipeline Overview

Raw Reads → QC → Reference Retrieval → Read Mapping → Assembly → BLAST Taxonomy Validation → ITS Extraction

All pipeline steps can be run sequentially using the run_pipeline.sh script.

Installing the package puts one command per pipeline step on your PATH, so steps can also be run individually from any directory:

Every command takes --help. Each is also runnable as a module, e.g. python -m its_fun_2_map.fastp_module --help.

Step 1: Quality Control (itsfun-qc)

Quality control processing of raw paired-end reads using fastp with a two-stage approach.

• Input: Raw paired-end sequencing reads
• Output:
  • {sample}_trimmed_1.fq, {sample}_trimmed_2.fq - trimmed paired reads
  • {sample}_merged.fq - merged overlapping reads
  • {sample}_unmerged_1.fq, {sample}_unmerged_2.fq - reads that couldn't be merged
  • {sample}_trim.json, {sample}_merge.json - QC metrics
  • {sample}_overlaps.html - overlap distribution visualisation
  • fastp_summary.csv - aggregated summary statistics
• Features:
  • Adapter trimming and filtering (auto-detection for PE reads)
  • Quality filtering (Q30 threshold by default, configurable via --qualified_quality_phred)
  • Read merging (Minimum merged length: 30 bp)
  • Read statistics generation

Parameters:

Note on threading. --threads and --fastp_threads multiply. Total worker threads is approximately --threads x --fastp_threads, so the default --threads 4 uses around 12 threads, not 4. The effective figure is written to the log at startup. On a scheduler, set --threads $SLURM_CPUS_PER_TASK --fastp_threads 1 to make the CPU request and actual usage match.

Note on --fastp_extra_args. This argument takes one quoted string:

--fastp_extra_args "--trim_front1 10 --trim_front2 10"

It cannot be used to set -q/--qualified_quality_phred or -w/--thread, which the module sets itself - use --qualified_quality_phred and --fastp_threads instead. Supplying either via --fastp_extra_args causes the run to exit with an error before processing begins.

Step 2: Pseudo-Reference Sequence Retrieval (itsfun-refs)

Retrieves taxonomically relevant reference sequences from the UNITE database using NCBI taxonomy matching.

Process:

1. Queries NCBI Entrez API for taxonomic lineage of each specimen
2. Matches NCBI taxonomy against UNITE database taxonomy
3. Extracts sequences based on user-specified parameters
4. Supports fallback to higher taxonomic ranks if matches not found

Key Features:

• Diversity mode: Distributes sequences across child taxa, prioritising target species matches (up to 25% of requested sequences)

• Traverse mode: Searches up the taxonomic hierarchy to collect target number of sequences

• Flexible matching: Case-insensitive matching with normalisation for species name variations

• Input: Tracking sheet with specimen taxonomy

• Output:

  • {sample_id}_seed.fasta - reference sequences for each specimen
  • Summary CSV with retrieval statistics

• Database: UNITE public (Full "UNITE+INSD") database

Parameters:

Step 3: Read Mapping & Baiting (itsfun-map)

Maps quality-filtered reads to retrieved reference sequences using BWA to enrich for target ITS regions. Supports both modern DNA (BWA-MEM) and ancient/historical DNA (BWA-ALN) alignment strategies.

Alignment Algorithms:

Process:

1. Automatically indexes reference FASTA files if not already indexed
2. Maps merged reads (single-end mode)
3. Maps unmerged read pairs (paired-end mode)
4. Extracts mapped reads to FASTQ format
5. Repairs paired-end read files to ensure synchronisation
6. Generates mapping statistics

Parameters:

Note on threading. As in Step 1, --threads and --bwa_threads multiply. Total worker threads is approximately --threads x --bwa_threads, so the defaults use around 32 threads, not 4. On a scheduler, set --threads $SLURM_CPUS_PER_TASK --bwa_threads 1 to make the CPU request and actual usage match.

Input:

• Quality-filtered FASTQ files from Step 1 (*_merged.fq, *_unmerged_1.fq, *_unmerged_2.fq)
• Per-sample reference FASTA files from Step 2 (*_seed.fasta)

Output:

• {sample}_mapped.fastq - mapped merged reads
• {sample}_mapped_unmerged_1.fastq, {sample}_mapped_unmerged_2.fastq - mapped paired reads (repaired)
• {sample}_mapped_flagstats.txt - alignment statistics for merged reads
• {sample}_unmerged_1_flagstats.txt, {sample}_unmerged_2_flagstats.txt - alignment statistics for paired reads
• mapping_summary.csv - aggregated mapping metrics for all samples

Step 4: Contig Assembly (itsfun-assemble)

Assembles mapped reads into contigs using SPAdes with a multi-stage fallback strategy to maximise assembly success.

Assembly Strategy:

Process:

1. Attempts assembly with merged reads using default k-mer sizes
2. If assembly fails or produces no valid scaffolds, retries with k=21 only
3. If still unsuccessful, includes unmerged paired reads in the assembly
4. Calculates assembly metrics including N50, scaffold counts, and lengths

Input:

• Mapped merged reads from Step 3 (*_mapped.fastq)
• Original unmerged reads from Step 1 (*_unmerged_1.fq, *_unmerged_2.fq)

Output:

• {sample}.spades.out/ - SPAdes output directory containing:
  • scaffolds.fasta - assembled scaffolds
  • Assembly graphs and logs
• assembly_summary.csv - metrics for all samples including:
  • Assembly method used (INITIAL/FALLBACK1/FALLBACK2)
  • Number of scaffolds
  • Scaffold lengths
  • Mean scaffold length
  • N50 value

Step 5: BLASTn - Round 1 (General Database) (itsfun-blast1)

Searches assembled scaffolds against the general UNITE database for initial taxonomic identification.

Process:

1. Optionally creates BLAST database from reference FASTA
2. Locates scaffolds.fasta files for each sample in assembly directories
3. Runs BLASTn searches in parallel
4. Outputs results in tabular format (outfmt 6)

Input:

• Assembled scaffolds from Step 4 (*.spades.out/scaffolds.fasta)
• UNITE general release database
• Tracking sheet with sample IDs

Output:

• {sample}_blast.tsv - BLAST results in tabular format with columns:
  • qseqid, sseqid, pident, length, mismatch, gapopen, qstart, qend, sstart, send, evalue, bitscore

Database: UNITE general release (sh_general_release_dynamic)

Step 5.5: BLASTn Round 1 Parsing & Taxonomic Validation (itsfun-blast1-parse)

Parses BLAST Round 1 output and validates taxonomic assignments against expected taxonomy with cascading taxonomic rank matching.

Taxonomic Validation Strategy: The parser attempts to match BLAST hits against expected taxonomy using a cascading approach:

1. Species-level match - checks if hits match expected species
2. Genus-level match - if no species match, checks genus
3. Family-level match - if no genus match, checks family
4. If no matches at any level, the sample is marked as FAIL

Filtering Criteria:

Validation Outcomes:

• PASS: Exactly one contig matches expected taxonomy after filtering
• FAIL: Zero contigs match, or multiple contigs remain after filtering

Input:

• BLAST TSV files from Step 5
• Taxonomy CSV with expected taxonomic assignments (must contain ID column and taxonomic hierarchy: Kingdom, Phylum, Class, Order, Family, Genus, Species)
• Assembly directory for contig extraction

Output:

• blast_validation_summary.csv - comprehensive validation results including:
  • ID - sample identifier
  • blast_round1_found_taxon - taxonomic level where match was found
  • blast_round1_expected_taxonomy - full expected taxonomy string
  • blast_round1_n_contigs_in - total contigs in assembly
  • blast_round1_n_contigs_hits - contigs matching expected taxonomy
  • blast_round1_contigs - contig IDs that matched
  • blast_round1_hit_taxonomy - taxonomy of best BLAST hit
  • blast_round1_hit_coords - alignment coordinates
  • blast_round1_correct_taxonomy - PASS/FAIL status
• filtered_{sample}_blast.tsv - filtered BLAST results for passing samples
• {sample}_parsed_contig.fasta - extracted contig sequences for passing samples

Step 6: BLASTn - Round 2a (ITS2-Specific) (itsfun-blast2)

Searches validated contigs against the UCHIME ITS2 reference database for region-specific validation. This step uses seqkit grep to extract only the contigs that passed Round 1 validation before running BLAST.

Process:

1. Reads BLAST Round 1 TSV results to identify validated contig IDs
2. Extracts unique contig headers from each sample's Round 1 results
3. Uses seqkit grep to extract matching scaffolds from assembly FASTA files
4. Runs BLASTn against ITS2-specific UCHIME database
5. Processes samples in parallel

Input:

• Assembly scaffolds from Step 4 (*.spades.out/scaffolds.fasta)
• BLAST Round 1 TSV results from Steps 5/5.5
• UCHIME ITS2 reference database
• Tracking sheet with sample IDs

Output:

• {sample}_blast.tsv - ITS2-specific BLAST results
• Temporary matched scaffold FASTAs (cleaned up after processing)

Database: UCHIME ITS2 reference dataset

Step 6.5: BLASTn Parsing & Validation (ITS2) (itsfun-blast-parse)

Parses ITS2-specific BLAST results and applies multi-tier taxonomic validation with cascading fallbacks from family to genus to species level.

Multi-Tier Taxonomy Matching Logic: When expected family is present:

1. 0 family matches → FAIL
2. 1 family match → PASS
3. >1 family matches → Check genus level:
   • No expected genus → FAIL
   • 1 genus match → PASS
   • 0 or >1 genus matches → Check species level:
     • No expected species → FAIL
     • 1 species match → PASS
     • 0 or >1 species matches → FAIL

When expected family is empty (fallback to genus):

1. 0 genus matches → FAIL
2. 1 genus match → PASS
3. >1 genus matches → Check species level (same logic as above)

Input:

• ITS2 BLAST TSV files from Step 6
• Taxonomy CSV with expected taxonomic assignments
• Assembly directory for contig extraction

Output:

• taxonomy_validation_summary.csv - validation results including:
  • ID - sample identifier
  • expected_family - expected taxonomic family (or genus/species if family unavailable)
  • expected_taxonomy - full expected taxonomy string
  • n_contigs_in - total contigs in assembly
  • n_contigs_hits - contigs matching expected taxonomy
  • contigs - matching contig IDs
  • hit_taxonomy - taxonomy from best BLAST hit
  • its_coords - ITS region coordinates on contig
  • correct_taxonomy - PASS/FAIL status
  • contig_path - path to extracted contig FASTA
• {sample}_parsed_contig.fasta - extracted contig sequences for passing samples

Step 7: BLASTn - Round 2b (ITS1-Specific) (itsfun-blast2)

Parallel analysis searching against UCHIME ITS1 reference database. Uses the same workflow as Step 6 but with the ITS1-specific database.

Input:

• Assembly scaffolds from Step 4
• BLAST Round 1 TSV results from Steps 5/5.5
• UCHIME ITS1 reference database
• Tracking sheet with sample IDs

Output:

• {sample}_blast.tsv - ITS1-specific BLAST results

Database: UCHIME ITS1 reference dataset

Step 7.5: BLASTn Parsing & Validation (ITS1) (itsfun-blast-parse)

Parses ITS1-specific BLAST results using the same multi-tier taxonomic validation as Step 6.5.

Input:

• ITS1 BLAST TSV files from Step 7
• Taxonomy CSV with expected taxonomic assignments
• Assembly directory for contig extraction

Output:

• taxonomy_validation_summary.csv - ITS1 validation results (same format as Step 6.5)
• {sample}_parsed_contig.fasta - extracted contig sequences for passing samples

Step 8: ITS Primer Binding & Extraction (itsfun-extract)

Identifies primer binding sites and extracts ITS sequences from validated contigs using seqkit amplicon. Uses standard ITS primers from White et al. (1990) by default, with support for custom primer sets.

Default Primers (White et al. 1990):

Target Regions:

Process:

1. Reads sample IDs from tracking sheet
2. Locates FASTA files for each sample
3. Runs seqkit amplicon with each primer pair (max 2 mismatches, strict mode)
4. Generates extraction summary statistics

Input:

• Validated contig FASTAs from Steps 6.5/7.5 (*_parsed_contig.fasta)
• Tracking sheet with sample IDs
• Optional: Custom primers TSV and regions TSV Output:
• ITS1/{sample}_ITS1.fa - extracted ITS1 sequences
• ITS2/{sample}_ITS2.fa - extracted ITS2 sequences
• ITS_complete/{sample}_ITS_complete.fa - extracted complete ITS sequences
• extraction_summary.csv - extraction success per sample and region
• summary_report.txt - detailed extraction statistics

Custom Primers TSV Format: A header row is required, with the columns name and sequence (tab-separated):

name	sequence
FWD1	ATCGATCGATCG
REV1	GCTAGCTAGCTA

Custom Regions TSV Format: A header row is required, with the columns region, forward and reverse (tab-separated). The forward and reverse values must match name entries in the primers TSV:

region	forward	reverse
MyRegion	FWD1	REV1

--primers_tsv and --regions_tsv must be supplied together, or neither (to use the built-in ITS defaults).

Note on custom regions. summary_report.txt reports per-region counts and the "no regions extracted" list using whatever region names your regions TSV defines. The ITS-specific breakdown - complete ITS, ITS1 only, ITS2 only is emitted only when the region set includes ITS1, ITS2 and ITS_complete.

Step 9: Summary Metrics Aggregation (itsfun-summary)

Compiles comprehensive summary statistics across all pipeline steps, performs contig analysis to resolve ITS1/ITS2 results, and prepares final output FASTAs for submission.

Process:

1. CSV Discovery: Finds and merges all pipeline summary CSVs (fastp, UNITEd, mapping, assembly, BLAST, extraction)
2. Contig Analysis: Compares ITS1 and ITS2 results to identify:
   • Contigs with both ITS1 and ITS2 on the same scaffold
   • Overlapping coordinate detection (potential assembly issues)
   • Resolution when ITS1 and ITS2 are on different contigs
3. Decision Making: Runs R-based decision logic (its_decision_making.R) to determine final outcomes
4. FASTA Preparation: Renames and organises final FASTAs for submission. This includes new sequence headers (default header= <ID>|<ITS1/ITS2/ITS_complete/putative_ITS_region> for any passes)

Contig Analysis Logic:

• Both ITS1 and ITS2 PASS on same contig: Check for coordinate overlap
  • Non-overlapping → PASS (contiguous)
  • Overlapping → WARNING (potential issue)
• Both PASS on different contigs: Default to ITS2 contig
• Only ITS1 PASS: Use ITS1 contig
• Only ITS2 PASS: Use ITS2 contig
• Both FAIL: Mark as failed

Input:

• All pipeline summary CSV files from previous steps
• Optional: Naming TSV for custom FASTA header formatting

Output:

• final_results_dir/merged_summary_{date}.csv - comprehensive merged summary with columns:
  • UNITEd metrics (reference retrieval)
  • Trimming/merging metrics (fastp)
  • Assembly metrics (SPAdes)
  • BLAST Round 1 results
  • ITS1 and ITS2 BLAST results
  • Extraction results
  • Mapping statistics
  • Contig analysis results (overlapping_coords, same_contig)
  • Final decision (Final_outcome, Final_contig_desc, Final_contig)
• final_results_dir/pass_fastas/ - renamed FASTAs for samples that passed all validation
• final_results_dir/manual_verification/ - FASTAs requiring manual review

Naming TSV Format (optional):

SAMPLE001	Genus_species_voucher123
SAMPLE002	Genus_species_voucher456

Output Directory Structure

output_base/
├── 00_logs/                          # Detailed logs for each step
├── 01_fastp_processed/               # Step 1: Quality-filtered reads
│   ├── {sample}_trimmed_1.fq
│   ├── {sample}_trimmed_2.fq
│   ├── {sample}_merged.fq
│   ├── {sample}_unmerged_1.fq
│   ├── {sample}_unmerged_2.fq
│   ├── {sample}_trim.json
│   ├── {sample}_merge.json
│   ├── {sample}_overlaps.html
│   ├── fastp_summary.out             (long-formatted version of fastp_summary.csv)
│   ├── fastp_summary.csv
│   ├── trimmed_json_summary.csv
│   └── merged_json_summary.csv
├── 02_UNITEd/                        # Step 2: Retrieved reference sequences
│   ├── {sample}_seed.fasta
│   └── UNITEd_summary.csv
├── 03_mapped_reads/                  # Step 3: Reads mapped to references
│   ├── {sample}_mapped.fastq
│   ├── {sample}_mapped_unmerged_1.fastq
│   ├── {sample}_mapped_unmerged_2.fastq
│   ├── {sample}_mapped_flagstats.txt
│   ├── {sample}_unmerged_1_flagstats.txt
│   ├── {sample}_unmerged_2_flagstats.txt
│   └── mapping_summary.csv
├── 04_assemblies/                    # Step 4: Assembled contigs
│   ├── {sample}.spades.out/
│   │   ├── scaffolds.fasta
│   │   ├── contigs.fasta
│   │   └── ...
│   └── assembly_summary.csv
├── 05_blast_round1/                  # Step 5: Initial BLAST results
│   └── {sample}_blast_results.tsv
├── 05a_blast_parsed1/                # Step 5.5: Filtered and validated results
│   ├── blast_validation_summary.csv
│   ├── filtered_{sample}_blast_results.tsv
│   └── {sample}_parsed_contig.fasta
├── 06_blast_round2a-ITS2/            # Step 6: ITS2-specific BLAST results
│   └── {sample}_blast_results.tsv
├── 06a_blast_parsed2a-ITS2/          # Step 6.5: Validated ITS2 results
│   ├── blast_parser_its2_summary.csv
│   ├── {sample}_blast_results-all_hits.tsv
│   └── {sample}_parsed_contig.fasta
├── 07_blast_round2b-ITS1/            # Step 7: ITS1-specific BLAST results
│   └── {sample}_blast_results.tsv
├── 07a_blast_parsed2b-ITS1/          # Step 7.5: Validated ITS1 results
│   ├── blast_parser_its1_summary.csv
│   ├── {sample}_blast_results-all_hits.tsv
│   └── {sample}_parsed_contig.fasta
├── 08_its_primer_extraction/         # Step 8: Extracted ITS sequences
│   ├── ITS1/
│   │   └── {sample}_ITS1.fa
│   ├── ITS2/
│   │   └── {sample}_ITS2.fa
│   ├── ITS_complete/
│   │   └── {sample}_ITS_complete.fa
│   ├── extraction_summary.csv
│   └── summary_report.txt
└── final_results_dir/                # Step 9: Final aggregated outputs
    ├── merged_summary_{date}.csv
    ├── pass_fastas/
    │   └── {sample}_{description}.fasta
    └── manual_verification/
        └── {sample}_manual_verification_needed.fasta

Dependencies

All of the dependencies below are installed by the conda environment (its-fun-2-map.yaml) and the its-fun-2-map package - see Installation and quick start. They are listed here for reference.

Languages

• r-base 4.5.2 - keep on 4.5.x; rpy2's conda build encodes the R series, so bumping R alone leaves rpy2 with no matching build
• python 3.14 - minimum supported is 3.12, which seqpy-tools requires

Conda packages

BBMap supplies repair.sh, which itsfun-map uses to resynchronise mapped read pairs in Step 3; openjdk is its Java runtime. htslib is no longer listed - it is pulled in automatically as a library dependency of samtools.

Python packages (not built-in)

R packages

• jsonlite 2.0.0 - used by parse_fastp_json.R; its_decision_making.R needs only base R

Custom tool packages

• its_fun_2_map.its_fun_tools - shared helpers, installed as part of the its-fun-2-map package. No PYTHONPATH fiddling and no need to run from the repository directory.
• seqpy-tools 0.1.0 (available via pip) - imported by every module in the pipeline

Authors

M. Kamouyiaros & D. Parsons @ the Natural History Museum, London (NHMUK)