Tutorial

Last updated: 2026-07-25

This tutorial provides a step-by-step guide on how to utilize harmonisv for post-processing the results of various SV calling methods and conducting joint SV calling across multiple samples and methods. The necessary input data and scripts for this tutorial can be found in the test folder on GitHub.

In this tutorial, we use the following tools to showcase the usage of harmonisv:

• aligners: minimap2 and NGMLR
• SV callers: cuteSV, sniffles, and SVIM
• SV merging: jasmine

Table of Contents

• 1. SV discovery using multiple methods
• 2. Harmonize VCFs
• 3. Remove duplicated SVs
• 4. SV merging
• 5. SV re-genotyping
• 6. Combine per-sample SV genotyping results
• 7. Random forest model for SV filtering
• 8. Merge per-sample VCFs into a population VCF
• 9. Run time

1. SV discovery using multiple methods

harmonisv is designed to harmonize and integrate the results from different SV calling methods. Input VCF files must include the following information:

• Type of SV
• Length of SV
• Sequencing depth for both reference and alternative alleles

We have generated the example VCF files of HG002 chr22 in the test/raw folder:

• HG002.minimap2.sniffles.vcf
• HG002.minimap2.svim.vcf
• HG002.minimap2.cuteSV.vcf
• HG002.NGMLR.sniffles.vcf
• HG002.NGMLR.svim.vcf
• HG002.NGMLR.cuteSV.vcf

These VCF files were generated using the following commands:

ref="hs37d5.fa"
input="HG002.fasta"

# align long-read to reference
# minimap2
minimap2 -L -t 36 --MD -a -x map-pb $ref $input | samtools sort -@ 4 -o $bam

# NGMLR
ngmlr --bam-fix -t 40 -x pacbio \
-r $ref \
-q $input \
-o $bam

# SV calling
# sniffles (ver 2.0.6)
sniffles \
-i $bam \
-v $vcf \
--reference $ref \
--tandem-repeats "human_hs37d5.trf.bed"

# SVIM (ver 2.0.0)
svim alignment $outpath $input $ref \
--max_sv_size 1000000

# cuteSV (ver 2.0.3)
cuteSV \
--max_cluster_bias_INS 100 \
--diff_ratio_merging_INS 0.3 \
--max_cluster_bias_DEL 200 \
--diff_ratio_merging_DEL 0.5 \
--genotype \
$input $ref $vcf $workdir

2. Harmonize VCFs

The output VCF files from different SV calling methods usually have different formats in their INFO and FORMAT fields. To integrate the results from different methods, we need to harmonize the VCFs to a standard format. Here, we first use harmonize-header to combine the headers from all input VCFs.

# If one tag is defined in multiple VCFs
# the one from the reference VCF or the first VCF in the list will be used
dir_header="output/harmonize_header/"
[[ ! -d $dir_header ]] && mkdir -p $dir_header

ls raw/*NGMLR* > $dir_header/NGMLR_vcf_list.txt

harmonisv harmonize-header \
-i raw/HG002.minimap2.cuteSV.vcf,raw/HG002.minimap2.sniffles.vcf,raw/HG002.minimap2.svim.vcf \
-f $dir_header/NGMLR_vcf_list.txt \
-o $dir_header/harmonized_header.txt \
-r raw/HG002.minimap2.sniffles.vcf

Then, we can use harmonize to standardize the VCFs:

1. standardize VCF headers and tag names
2. normalize SV types (e.g., DUP:TANDEM to DUP)
3. extract depth and number of supporting reads to INFO/DP and INFO/RE, respectively
4. rename SV ID to make them unique across samples and methods. By using --id-prefix and --rename-id, a unique ID in the format of prefix.chr.svtype.number will be assigned to each SV.

dir_harmonize="output/harmonize/"
[[ ! -d $dir_harmonize ]] && mkdir -p $dir_harmonize

# sniffles
# note: INFO/DP = FORMAT/DR + FORMAT/DV
ls raw/*sniffles* | while read vcf; do
    name=$(basename $vcf)
    name=${name%.vcf}

    harmonisv harmonize \
    -i $vcf \
    -o ${dir_harmonize}/${name}.harmonized.vcf \
    --info SVTYPE,SVLEN,END,STRANDS=STRAND \
    --format-to-info RE=DV \
    --format-to-info-sum DP=DR,DP=DV \
    --header $dir_header/harmonized_header.txt \
    --header-str 'STRANDS,1,String,Strand orientation of supporting reads' \
    --id-prefix $name \
    --rename-id

done

# SVIM
# note: convert SVTYPE "DUP:TANDEM" and "DUP:INT" to "DUP"
ls raw/*svim* | while read vcf; do
    name=$(basename $vcf)
    name=${name%.vcf}

    harmonisv harmonize \
    -i $vcf \
    -o ${dir_harmonize}/${name}.harmonized.vcf \
    --info SVTYPE,SVLEN,END,RE=SUPPORT \
    --format-to-info DP=DP \
    --DUP DUP,DUP:TANDEM,DUP:INT \
    --header $dir_header/harmonized_header.txt \
    --header-str 'STRANDS,1,String,Strand orientation of supporting reads' \
    --id-prefix $name \
    --rename-id
done

# cuteSV
ls raw/*cuteSV* | while read vcf; do
    name=$(basename $vcf)
    name=${name%.vcf}

    harmonisv harmonize \
    -i $vcf \
    -o ${dir_harmonize}/${name}.harmonized.vcf \
    --info SVTYPE,SVLEN,END,RE \
    --format-to-info-sum DP=DR,DP=DV \
    --header $dir_header/harmonized_header.txt \
    --header-str 'STRANDS,1,String,Strand orientation of supporting reads' \
    --id-prefix $name \
    --rename-id
done

3. Remove duplicated SVs

SV calling methods may redundantly call the same SV. We can use jasmine to identify intra-sample duplicated SVs. The results can be found in test/dup_call

dir_dupcall="dup_call/"
[[ ! -d $dir_dupcall ]] && mkdir -p $dir_dupcall

ls output/harmonize/*harmonized.vcf | while read vcf; do
    name=$(basename $vcf)
    name=${name%.vcf}

    jasmine file_list=$vcf \
    out_file=${dir_dupcall}/${name}.dup_call.vcf \
    genome_file=hs37d5.fa \
    --comma_filelist \
    max_dist=200 \
    --allow_intrasample \
    --nonlinear_dist \
    --ignore_strand \
    --keep_var_ids 

    # write duplicated SV list
    bcftools query -f '%INTRASAMPLE_IDLIST\n' ${dir_dupcall}/${name}.dup_call.vcf \
    > ${dir_dupcall}/${name}.dup_call.txt
done

Then, we use represent to remove duplicated SVs by keeping the one with more supporting reads.

dir_dedup="output/dedup/"
[[ ! -d dir_dedup ]] && mkdir -p $dir_dedup

ls ${dir_harmonize}/*harmonized.vcf | while read vcf; do
    name=$(basename $vcf)
    name=${name%.vcf}

    harmonisv represent \
    -i $vcf \
    -o ${dir_dedup}/${name}.dedup.vcf \
    --merge ${dir_dupcall}/${name}.dup_call.txt \
    --by-max RE \
    --min-len-input 30 \
    --min-len-output 30
done

4. SV merging

After generating individual SV discovery call sets, we next identify non-redundant SVs across samples and methods by using jasmine. The results can be found in the test/sv_merge folder.

dir_merge="sv_merge/"
[[ ! -d $dir_merge ]] && mkdir -p $dir_merge

ls ${dir_harmonize}/*harmonized.vcf > ${dir_merge}/All_method.merge_vcf_list.txt

jasmine \
file_list=${dir_merge}/All_method.merge_vcf_list.txt \
out_file=${dir_merge}/All_method.merge.vcf \
--keep_var_ids \
--ignore_strand 

bcftools query -f '%IDLIST\n' ${dir_merge}/All_method.merge.vcf > ${dir_merge}/All_method.merge.txt

To select a representative SV from each group of redundant calls, we use the represent command to keep the one with the most frequent POS and SVLEN.

dir_represent="output/represent/"
[[ ! -d $dir_represent ]] && mkdir -p $dir_represent

harmonisv represent \
-f ${dir_merge}/All_method.merge_vcf_list.txt \
-o ${dir_represent}/All_method.representative.vcf \
--merge ${dir_merge}/All_method.merge.txt \
--by-freq \
--id-prefix All_method \
--save-id

5. SV re-genotyping

The VCF All_method.representative.vcf includes all non-redundant SVs discovered across methods and samples. To obtain a fully genotyped VCF for each sample, we re-genotype those SVs using sniffles and cuteSV. The results can be found in the test/force_call folder.

dir_force_call="force_call/"
[[ ! -d $dir_force_call ]] && mkdir -p $dir_force_call

# 1. re-genotyping
# sniffles
sniffles \
-i $bam \
-v "${dir_force_call}/HG002.minimap2.sniffles.vcf" \
--reference "hs37d5.fa" \
-t 40 \
--tandem-repeats "human_hs37d5.trf.bed" \
--genotype-vcf ${dir_represent}/All_method.representative.vcf

# cuteSV
cuteSV \
--max_cluster_bias_INS 100 \
--diff_ratio_merging_INS 0.3 \
--max_cluster_bias_DEL 200 \
--diff_ratio_merging_DEL 0.5 \
--genotype \
-t 40 \
-L -1 \
-Ivcf ${dir_represent}/All_method.representative.vcf \
$bam hs37d5.fa "${dir_force_call}/HG002.minimap2.cuteSV.vcf" $workdir

After re-genotyping, we need to harmonize them:

# 2. harmonize 
dir_force_call_harmonize="output/harmonize_force_call/"
# sniffles
ls ${dir_force_call}/*sniffles* | while read vcf; do
    name=$(basename $vcf)
    name=${name%.vcf}

    harmonisv harmonize \
    -i $vcf \
    -o ${dir_force_call_harmonize}/${name}.harmonized.vcf \
    --info SVTYPE,SVLEN,END \
    --format-to-info RE=DV \
    --format-to-info-sum DP=DR,DP=DV \
    --header $dir_header/harmonized_header.txt \
    --header-str 'STRANDS,1,String,Strand orientation of supporting reads'
done

# cuteSV
ls ${dir_force_call}/*cuteSV* | while read vcf; do
    name=$(basename $vcf)
    name=${name%.vcf}

    harmonisv harmonize \
    -i $vcf \
    -o ${dir_force_call_harmonize}/${name}.harmonized.vcf \
    --format-to-info-sum DP=DR,DP=DV \
    --header $dir_header/harmonized_header.txt \
    --header-str 'STRANDS,1,String,Strand orientation of supporting reads' \
    --keep-all
done

6. Combine per-sample SV genotyping results

Now we have 10 VCFs per sample: 4 discovery call sets and 6 re-genotyping call sets. We can use the genotype command to merge the results across all VCFs, producing the final per-sample VCF for downstream analysis.

dir_genotype="output/genotype/"
[[ ! -d $dir_genotype ]] && mkdir -p $dir_genotype

$harmonisv genotype \
-i ${dir_represent}/All_method.representative.vcf \
-f manifest_genotype.txt \
-o ${dir_genotype}/HG002.representative.genotyped.vcf \
--sample HG002

Particularly, the manifest_genotype.txt is a tab-separated file with the following columns:

• file: path to SV discovery and force calling VCF/BCF file
• sample: sample ID
• aligner: aligner used to generate the SV call set
• caller: SV caller used to generate the SV call set
• if_force_call: whether the SV call set is generated by force-calling (1: True, 0: False)
• info (optional): additional INFO tags beyond DP and RE to include in the output VCF

file	sample	aligner	caller	is_force_call	info
output/harmonize_force_call/HG002.NGMLR.cuteSV.harmonized.vcf	HG002	NGMLR	cuteSV	1	PRECISE,CIPOS,CILEN
output/harmonize_force_call/HG002.NGMLR.sniffles.harmonized.vcf	HG002	NGMLR	sniffles	1
output/harmonize_force_call/HG002.minimap2.cuteSV.harmonized.vcf	HG002	minimap2	cuteSV	1	PRECISE,CIPOS,CILEN
output/harmonize_force_call/HG002.minimap2.sniffles.harmonized.vcf	HG002	minimap2	sniffles	1
output/harmonize/HG002.NGMLR.cuteSV.harmonized.vcf	HG002	NGMLR	cuteSV	0
output/harmonize/HG002.NGMLR.sniffles.harmonized.vcf	HG002	NGMLR	sniffles	0
output/harmonize/HG002.NGMLR.svim.harmonized.vcf	HG002	NGMLR	svim	0
output/harmonize/HG002.minimap2.cuteSV.harmonized.vcf	HG002	minimap2	cuteSV	0
output/harmonize/HG002.minimap2.sniffles.harmonized.vcf	HG002	minimap2	sniffles	0
output/harmonize/HG002.minimap2.svim.harmonized.vcf	HG002	minimap2	svim	0

7. Random forest model for SV filtering

The per-sample VCF may contain false positive SVs. To filter out those false positives, we can use the filter to perform estimate the quality of each SV calls for all per-sample VCFs specified in --manifest. The model will be trained on known positive and negative SV sites in --train-sites and applied to all SV sites. If --bench-sites is provided, the model will also be evaluated on the benchmark sites (excluded from --train-sites). See the filter section for more details and test/filter for example input files.

Note: While the model only process one type of SV at a time, it includes all SVs in the output VCF. It is suggested to run the model on each SV type, do bcftools view -i and bcftools filter to extract SVs and assign filter status (e.g., by RF_SCORE), and then merge all SV types into a final VCF.

dir_filter="output/filter/"
[[ ! -d $dir_filter ]] && mkdir -p $dir_filter

harmonisv filter \
--manifest manifest_filter.txt \
--output ${dir_filter}/model_training \
--sv-type INS \
--feature "SVLEN,MEAN_VAF,STD_VAF,DP_MINIMAP2_CUTESV,VAF_MINIMAP2_CUTESV,VAF_MINIMAP2_SVIM,DP_NGMLR_CUTESV,VAF_NGMLR_CUTESV,VAF_NGMLR_SVIM" \
--train-sites filter/train_sites.txt \
--bench-sites filter/bench_sites.txt \
--bench-sample HG002 \
--train-size 0.8 \
--n-iter 10 \
--max-depth 7:9:1 \
--min-samples-leaf 1,5,10 \
--seed 42

harmonisv filter \
--manifest manifest_filter.txt \
--output ${dir_filter}/filtered_SV \
--apply-model ${dir_filter}/model_training.INS.model \
--sv-type INS \
--feature "SVLEN,MEAN_VAF,STD_VAF,DP_MINIMAP2_CUTESV,VAF_MINIMAP2_CUTESV,VAF_MINIMAP2_SVIM,DP_NGMLR_CUTESV,VAF_NGMLR_CUTESV,VAF_NGMLR_SVIM" \
--seed 42

8. Merge per-sample VCFs into a population VCF

After filtering, we have per-sample VCFs with confidence scores. The sample2pop command merge all per-sample VCFs specified in --manifest into a population VCF. Variants are mergedd across samples by their ID. INFO fields are merged using user-defined rules. In this example, SVLEN and SVTYPE use the first value (they should be identical across samples), sequencing depths (DP*, RE*) are summed, number of supporting callers/methods (SUPP_*) are averaged, MAX_RE takes the minimum, and the INFO/RF_SCORE takes the maximum. The --info-to-format option moves per-sample INFO/RF_SCORE into FORMAT/RF_SCORE so each sample retains its individual score. --filter-GT ensures variants that failed per-sample filters are set to ./.. See the sample2pop documentation for a full list of merge rules and options.

dir_pop="output/population/"
[[ ! -d $dir_pop ]] && mkdir -p $dir_pop

harmonisv sample2pop \
--manifest manifest_sample2pop.txt \
--outvcf ${dir_pop}/population.vcf.gz \
--filter-GT \
--info-first SVTYPE,SVLEN,END \
--info-sum "DP*,RE*" \
--info-avg "SUPP_*" \
--info-min MAX_RE \
--info-max RF_SCORE \
--info-to-format RF_SCORE \
--keep-format GT

9. Run time

The runtime of the above commands on a PC with an i7-14700 is shown below (excluding time for SV discovery and merging by other tools).

real    0m35.615s
user    1m27.436s
sys     1m37.989s