Performance evaluation
To examine the impact of known SVs, we evaluate the performance of SVision on HG002 with two different training sets:
- Original training SVs
- SVs outside of HG002 passed SVs. The passed SVs are used because the Truvari evaluates with
--passonlyoption.
There are 12,745 passed HG002 SVs, and they are based on hg19 coordinate system. Since our training set is obtained on GRCh38, we lifted over the passed HG002 SVs and obtained 12,683 SVs. We then use bedtools to test the overlaps between our training set and 12,683 HG002 SVs.
The performance of using different training sets are shown below:
| Training set | Precision | Recall | F1 |
|---|---|---|---|
| All SVs | 0.9469 | 0.8594 | 0.9011 |
| SVs outside HG002 | 0.9465 | 0.8594 | 0.9009 |
Besides visual confirmation of CSVs in NA12878, we further applied dipcall and PAV to detect SVs from phased assembly.
The phased assembly is obtained from HGSVC recent release. We then used default parameters to run dipcall and PAV, and their calls were compared to 18 visual confirmed CSVs (results download link) used in performance evaluation.
Except the validation introduced in the paper, we further added trio analysis to show the Mendelian error of 80 CSVs detected by SVision.
The HiFi data of HG00732 and HG00731 could be downloaded from HGSVC ftp site.
We applied the same parameters as we did for HG00733 to align and call SVs from HG00732 (VCF download link) and HG00731 (VCF download link). This analysis (trio detection summary download link) suggests a Mendelian error rate of 1.25% for the trio.
| HG00731 | HG00732 | HG00733 | Count |
|---|---|---|---|
| 1 | 1 | 1 | 40 |
| 1 | 0 | 1 | 18 |
| 0 | 1 | 1 | 21 |
| 0 | 0 | 1 | 1 |
We followed the procedure used in CuteSV Supplementary Note to evaluate the detection of simple SVs from HG002.
pbmm2 index reference.fa reference.mmi
pbmm2 align reference.mmi reads.fa HG002.pbmm2.bam --sort --rg '@RG\tID:HG002' --sample HG002
samtools view -Sb -s ratio HG002.pbmm2.bam > HG002.pbmm2.subset.bam
samtools index HG002.pbmm2.subset.bam
# For ONT reads
cuteSV HG002.pbmm2.ont.bam HG002.cutesv.ont.vcf ./workdir -l 50 --max_cluster_bias_INS 100 --diff_ratio_merging_INS 0.3 --max_cluster_bias_DEL 100 --diff_ratio_fitering_DEL 0.3
# For HiFi reads
cuteSV HG002.pbmm2.hifi.bam HG002.cutesv.hifi.vcf ./workdir -l 50 --max_cluster_bias_INS 1000 --diff_ratio_merging_INS 0.9 --max_cluster_bias_DEL 1000 --diff_ratio_fitering_DEL 0.5
sniffles -l 50 -m HG002.pbmm2.ont.bam -v ./HG002.sniffles.ont.vcf --genotype
sniffles -l 50 -m HG002.pbmm2.hifi.bam -v ./HG002.sniffles.hifi.vcf --genotype
pbsv discover HG002.pbmm2.bam ./sigs.svsig.gz
pbsv call referece.fa ./sigs.svsig.gz HG002.pbsv.vcf
svim alignment --min_sv_size 50 ./workdir HG002.pbmm2.bam reference.fa
cat {workdir/final_results.vcf} | awk '{ if($1 ~ /^#/) { print $0 } else { if($5=="<DEL>" || $5=="<INS>") { print $0 }}}' | grep -v 'SUPPORT=1;\|SUPPORT=2;\|SUPPORT=3;\|SUPPORT=4;\|SUPPORT=5;\|SUPPORT=6;\|SUPPORT=7;\|SUPPORT=8;\|SUPPORT=9;' | sed 's/DUP:INT/INS/g' | sed 's/DUP:TANDEM/INS/g' > HG002.svim.vcf
bgzip HG002.tools.vcf && tabix HG002.tools.vcf.gz
truvari -f reference.fa -b HG002_SVs_Tier1_v0.6.vcf.gz --includebed HG002_SVs_Tier1_v0.6.bed -o tools_truvari --passonly --giabreport -r 1000 -p 0.00 -c HG002.tools.vcf.gz
The performance is shown in the figure below.
SVision is directly applicable to assembled genomes due to image coded sequence differences. Here we show the demo about running and benchmarking with the most recently published PAV calls made by PAV caller.
-
Download the HiFi reads
-
Download the phased assemblies
-
Download the PAV calls: we are using the freeze2 calls published by HGSVC in the performance benchmarking for NA12878 and HG00733.
Contig alignment with minimap2, key parameters are shown below:
-x asm20 -m 10000 -z 10000,50 -r 50000 --end-bonus=100 --secondary=no -a --eqx -Y -O 5,56 -E 4,1 -B 5
HiFi reads are also aligned with minimap2.
- We only selected autosome PAV calls to benchmark the performance, where all insertions and deletions are included.
- PAV calls are used as haploid benchmark and diploid benchmarks. The haploid benchmark is used to evaluate reads alignment based detection, while diploid benchmarks (i.e., H1 and H2) are used to evaluate contig alignment based detection.
- For haploid evaluation, PAV calls are separated by SV length and saved to different BED files, including [50,100), [100,300), [300,1000) and [1000,6000).
- Contig mode is activated by setting
--contig - The output from reads and contigs are directly used without any post-processing.
- SVision read- and contig-calls are separate by SV length.
- Evaluating SVision calls with 50% reciprocal overlaps.
bedtools intersect -c -a pav.bed -b svision.bed -f 0.5 -r
This evaluates the detection of haploid PAV benchmark.
This evaluates the detection of haploid PAV benchmark.
We also tested SVision scalability to long-read assemblies for SVs of different length. The performance was under our expectation, where longer read input would improve the detection performance. The number on x-axis is the maximum length of SVs considered in the evaluation, 'All' indicates all events detected by SVision is used in evaluation.
For the diploid comparison, support read number is one for SVision, CuteSV, Sniffles and pbsv, while SVIM-ASM is specifically designed for contig SV detection and run with default settings.
