cuteSV

Before generating simulated datasets, you must decompress the provided BED files, set up a working directory, and install the VISOR dependency.

1. Decompress the simulation BED files (e.g., sim_del.bed.gz, sim_ins.bed.gz, etc.).
2. Create a workdir directory and move the decompressed .bed files and LASeR.bed into it.
3. Clone the VISOR repository and install it using python3 setup.py install.

# Decompress files
gzip -d sim_del.bed.gz
gzip -d sim_ins.bed.gz
gzip -d sim_inv.bed.gz
gzip -d sim_dup.bed.gz
gzip -d sim_tra.bed.gz

# Setup work directory
mkdir workdir && cd workdir
mv ../sim_*.bed ../LASeR.bed ./

# Install VISOR dependency
git clone https://github.com/tjiangHIT/VISOR.git
cd VISOR && python3 setup.py install && cd ..

The Force Calling Benchmark is a demo procedure to evaluate SV calling performance by performing regenotyping against a ground truth set (e.g., NIST GIAB SV v0.6).

Workflow Overview:

1. Setup Environment: Install sniffles, cuteSV, svjedi, and truvari via Conda.
2. Prepare Data: Download NIST ground truth (BED/VCF), reference genome (hg19/hs37d5), and PacBio HiFi alignments.
3. Process Alignments: Convert BAM to FASTA, align with pbmm2, and use samtools calmd to prepare the alignment.
4. Call Variants: Run different callers (sniffles, cuteSV, svjedi) on the processed BAM.
5. Format for Evaluation: Sort, compress (bgzip), and index (tabix) the resulting VCFs to make them compatible with truvari.
6. Compare: Use truvari bench to compare the called variants against the NIST ground truth.

To optimize SV detection, use specific parameter sets depending on your sequencing technology:

For PacBio CLR data:

--max_cluster_bias_INS 100 --diff_ratio_merging_INS 0.3 --max_cluster_bias_DEL 200 --diff_ratio_merging_DEL 0.5

For PacBio CCS (HIFI) data:

--max_cluster_bias_INS 1000 --diff_ratio_merging_INS 0.9 --max_cluster_bias_DEL 1000 --diff_ratio_merging_DEL 0.5

For ONT data:

--max_cluster_bias_INS 100 --diff_ratio_merging_INS 0.3 --max_cluster_bias_DEL 100 --diff_ratio_merging_DEL 0.3

To use cuteSV for structural variant (SV) calling and benchmarking, you can set up a dedicated Conda environment. This environment includes cuteSV, sniffles, svjedi, truvari, and standard bioinformatics tools like samtools and tabix.

conda create -n test_fc python=3.10
conda activate test_fc
conda install -c bioconda sniffles==2.0.7 cuteSV==2.1.0 svjedi==1.1.6 truvari==3.5.0 samtools tabix

You can install cuteSV using pip, conda, or by installing from the source repository.

Using pip:

pip install cuteSV

Using conda (Bioconda):

conda install -c bioconda cutesv

From source:

git clone https://github.com/tjiangHIT/cuteSV.git && cd cuteSV/ && python setup.py install

The project uses a log-space normalization for genotype likelihoods to ensure numerical stability.

1. log10sumexp: Calculates the log-sum-exp of a list of log10 probabilities. This is used to find the total probability in log space.
2. normalize_log10_probs: Adjusts the input log10 probabilities by subtracting the log-sum-exp value, ensuring the resulting probabilities sum to 1 (in log space, this means the maximum value is 0.0).

The project follows specific rules for anchoring structural variant (SV) coordinates to ensure compatibility with the VCF 4.2 standard:

• DEL, INS, DUP, INV: These typically anchor on the base before the event. A 0-based cut point is converted to a 1-based position using vcf_pos.
• BND (Breakends):
  • Types A/B: Anchor on the base before the breakpoint. The ALT string carries the anchor base in front of the bracket pair.
  • Types C/D: Anchor on the base after the breakpoint. These are written one base further on compared to A/B types.
• Special Case (Cut point 0): If a cut point is 0, there is no base to its left. In this case, the record anchors on the base after the event (base 1).

cuteSV is a tool for structural variation (SV) detection from long-read sequencing data. The primary command requires a sorted BAM file (from NGMLR or Minimap2), a reference genome in FASTA format, an output VCF file path, and a work directory for distributed jobs.

Basic Usage Pattern: cuteSV [BAM] [REFERENCE] [OUTPUT] [WORK_DIR] [OPTIONS]

Recommended Parameter Presets by Sequencing Technology:

cuteSV input.bam reference.fasta output.vcf ./work_dir --threads 16

After generating SV calls from various tools, use truvari bench to compare them against a ground truth VCF.

Usage Pattern: Iterate through your tool-specific VCF files and run the benchmark command. Ensure the ground truth VCF is compressed and indexed with tabix.

Example Command:

truvari bench -b <ground_truth>.vcf.gz -c <tool_call>.vcf.gz --includebed <truth>.bed -o <output_dir> -p 0 -r 2 -P 1 --sizemax 1000000

truvari bench -b giab/HG002_SVs_Tier1_v0.6.filter.vcf.gz -c tools/$i.sort.vcf.gz --includebed giab/HG002_SVs_Tier1_v0.6.bed -o cmp -p 0 -r 2 -P 1 --sizemax 1000000

Use the VISOR LASeR command to generate simulated sequencing alignments. This command takes a reference genome and a modified donor genome to produce alignment data at various coverage levels.

Arguments:

• -g (or --genome): The reference genome FASTA file.
• -s (or --subject): The modified donor genome FASTA file.
• -bed: The BED file containing the variation information.
• -o (or --output): The output prefix/name for the generated data.
• -c (or --coverage): The target coverage (e.g., 5, 10, 20, 30).
• --threads: Number of CPU threads to use.
• --noaddtag: Flag to prevent adding tags to the output.

# Example: Generate 5x coverage deletion data
VISOR LASeR -g reference_genome.fasta -s donor_genome_del -bed LASeR.bed -o data_del_5x -c 5 --threads 16 --noaddtag