Virasca

Virasca (Virus Scaffolding) is a CLI tool for scaffolding/genome assembly of viral contigs, using Snakemake for workflow management.

Installation

It is recommended to use micromamba (or conda/mamba) to manage dependencies.

1. Clone virasca repository and submodules:

   git clone https://github.com/WallauBioinfo/virasca.git
   cd virasca

2. Create the environment:

   micromamba env create -f environment.yml

3. Activate the environment:

   micromamba activate virasca
   conda config --set channel_priority strict

Usage

Configure Database

This command downloads and configures the necessary databases (NCBI Virus genomes and taxonomy).

virasca configure-database

Options:

• --cores: Number of cores to use (default: 1).
• --dry-run: Print the workflow without executing.

Using Custom Databases

Instead of using the default RefSeq viral databases downloaded via configure-database, you can use your own custom databases in Virasca. Ensure you have the following two files ready:

1. Database FASTA (--database-seq): A standard FASTA file containing all reference sequences you want Virasca to align against and use for reference-guided assembly/scaffolding.
2. Metadata TSV (--database-metadata): A tab-separated file mapping each FASTA accession to its taxonomic lineage. This file is strictly required for Virasca's taxonomic grouping logic. It must contain at least the following columns:
   • accession: The exact sequence ID present in the FASTA headers.
   • Taxonomic level columns (e.g., species_name, genus_name, family_name, order_name, class_name, phylum_name). If a lineage level is unknown, you can leave it empty or write NA.

When running your analysis, simply replace the default paths with your custom files:

virasca run \
  --input contigs.fasta \
  --reads-r1 R1.fastq.gz \
  --reads-r2 R2.fastq.gz \
  --output custom_analysis_out \
  --database-seq path/to/custom_sequences.fasta \
  --database-metadata path/to/custom_metadata.tsv \
  --tax-level species

Run Analysis

This command runs the viral assembly and identification pipeline.

virasca run --input <scaffolds.fasta> --reads-r1 <R1.fastq.gz> --reads-r2 <R2.fastq.gz>

Arguments:

• --input: Input FASTA file with contigs/scaffolds.
• --reads-r1: Path to R1 reads (FASTQ).
• --reads-r2: Path to R2 reads (FASTQ).

Options:

• --database-seq: Path to blast database FASTA (required).
• --database-metadata: Path to blast database metadata file (required).
• --output: Output directory for analysis results (required).
• --use-virseqimprover: Enable Virseqimprover step for iterative assembly improvement.
• --tax-level: Taxonomic level for reference selection (default: species). Choices: phylum, class, order, family, genus, species.
• --fastp-threads: Number of threads for fastp (default: 4).
• --bwa-threads: Number of threads for BWA (default: 4).
• --min-len: Minimum read length for fastp (default: 50).
• --trim-len: Trim length for fastp (default: 0).
• --min-base-quality: Minimum base quality for fastp (default: 20).
• --mapping-quality: Minimum mapping quality for iVar consensus (default: 20).
• --min-depth: Minimum depth to call consensus in iVar (default: 1).
• --cores: Number of cores for Snakemake (default: 4).
• --dry-run: Print the workflow without executing.

Final Outputs

The virasca run command will generate a detailed output folder containing:

• output.tsv: A tabular summary of each input contig/scaffold, its taxonomic classification, sequence identity, coverage, and the specific refinement action taken.
• output.info.md: A detailed Markdown report tracing the exact fate of each sequence throughout the pipeline. It aggregates sequences by fate and includes a Consensus Mapping table that correlates the generated assembly consensus sequences directly back to the original input contigs that triggered them.
• read_assembly_*/: Directories containing the reference-guided assembly results (cleaned reads, BAM mappings, and the consensus.fa files generated by iVar) for sequences classified as intermediario or singleton unrepresented incompleto sequences.
• ragtag_*/: Directories containing the scaffolding output generated by RagTag for sequences classified as incompleto (groups that have $\ge 2$ contigs assigned to their taxonomic group).

Testing

To verify the installation and workflow structure (dry-run):

# Ensure you are in the virasca environment
micromamba activate virasca

# Run the test script
python3 tests/test_cli.py

Workflow Logic

graph TD
    Start["Input: Contigs  & Reads"] --> Decision{"Use Virseqimprover?"}
    Decision -- Yes --> Virseq["Virseqimprover (Iterative Extension)"]
    Decision -- No --> Blast["BLASTn vs Database"]
    Virseq --> Blast
    Blast --> Classify["Classification"]
    
    Classify -->|Coverage >= 98%| Comp["Completo"]
    Classify -->|Identity >= 90%| Inter["Intermediario"]
    Classify -->|Identity >= 60%| Inc["Incompleto"]
    
    Comp --> Output["Final Output"]
    
    Inter --> RefSel1["Select Best Reference"]
    RefSel1 --> Assembly["Reference-Guided Assembly (fastp -> bwa -> ivar)"]
    Assembly --> Output
    
    Inc --> RefSel2["Select Best Reference"]
    RefSel2 -->|>= 2 Contigs| RagTag["Scaffolding (RagTag)"]
    RefSel2 -->|< 2 Contigs| Assembly
    RagTag --> Output
    
    style Start fill:#f9f,stroke:#333,stroke-width:2px
    style Output fill:#9f9,stroke:#333,stroke-width:2px
    style Decision fill:#ff9,stroke:#333,stroke-width:2px

Loading

Configure Database

1. Download: Downloads RefSeq viral genomes using NCBI Datasets CLI.
2. Taxonomy: Downloads and configures NCBI taxonomy dump using TaxonKit.
3. Metadata: Processes metadata to link genomic accessions with full taxonomic lineage (Phylum to Species).
4. Database: Creates a local BLAST nucleotide database.

Run Analysis

1. Preprocessing (Optional): If --use-virseqimprover is set, runs Virseqimprover to iteratively extend input contigs using the provided reads.
2. Identification:
   • Aligns contigs against the local viral database using BLASTn.
   • Classifies each hit based on Identity (pident) and Coverage (scov):
     • Completo: Coverage >= 98%
     • Intermediario: Identity >= 90% (and not complete)
     • Incompleto: Identity >= 60% (and not intermediate)
3. Refinement Strategy:
   • For "Intermediario" Genomes: Selects the best reference (highest bitscore) per taxonomic group. Performs reference-guided assembly:
     • Read cleaning (fastp)
     • Mapping (BWA MEM)
     • Consensus generation (iVar)
   • For "Incompleto" Genomes: Selects the best reference per taxonomic group. Performs scaffolding:
     • Scaffolds input contigs against the reference using RagTag.
4. Output: Generates a final summary report linking original contigs, classification status, and refinement results.