Hi everyone! I am trying to use nf-core/variantbenchmarking with GIAB data and rtg-tools option for benchmarking and I am getting this error: Error: An IO problem occurred: “Invalid file pointer: 14104788689 for HG001_GRCh37_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_highconf_PGandRTGphasetransfer_appchr_snps_cov.reheader.vcf.gz” I am not very familiar with nextflow so any help would be much appreciated! Also, are there any test files and test commands available for this pipeline? Apart from the samplesheet I could not find any working example in github. Thanks!

Hello Laura, Can you share your running configurations? I can try to help you. all the test cases are in conf/test directory. All of them are pulling the data automatically so that you can run some examples. One example is:

nextflow run nf-core/variantbecnhmarking -profile germline_small,docker

Hi Kübra, thanks so much for the quick reply. I was able to run some conf/test that you pointed out through the -profile param. So I guess I had some mispecified param or path in my running path. I also managed to run the test test_gagh.config test (downloading the test_gagh.config) and using the -c param (I think the test is not yet included). THis is my bash script: nextflow run nf-core/variantbenchmarking -r dev -profile docker -c ‘./test_ga4gh.config’ This is working now (ive downloaded the test_ga4gh.config from the varbenchmarking github). However, I think the output summary tables are incorrect. e.g. resutls/small/test1/benchmarks/happy/test1.HG002.strelka.summary.csv According to hap.py page, when running hap.py with rtg vcfeval engine, FP.AL column should be 0 and is non-zero in the output from test_ga4gh.config test (https://github.com/Illumina/hap.py/blob/master/doc/happy.md#comparison-engines). Do you know why? Has the vcfeval comparison strategy changed?

That sounds interesting and I have never tested to see myself. Let me investigate more.

Thank you! Ive tried to investigate myself but hap.py has not been updated since pre-covid as far as I can tell and thus I cannot build from binaries to see if that installation gives a different output

Hello Laura, I did experiment this comparing with and without rtg vcfeval engine. It is interesting since FP.AL changes when you use the engine but it is not being 0. I also run it without prepy normalization and also didnt change the result at all. To be honest, I am not sure why this weirdness happens. I open this thread https://github.com/Illumina/hap.py/issues/201 lets see if they have a posable explanation 😅

Hi Kübra, thanks so much for looking into this! Ive also tried to follow the commands in this paper supplementary info using the bioconda happy (https://doi.org/10.1038/s41588-018-0316-4): 2.3 Software, versions, and parameters for the variant calling benchmarks For evaluating the accuracy of the variant calls (against simulated data, GiaB truth-sets, and the SNP Array truth-sets) we used vcfeval version 3.6, followed by hap.py version (this is the same setup as used by the recent PrecisionFDA Truth Challenge: https://precision.fda.gov/challenges/truth/results).

rtg vcfeval -t $SDF -o $OUTPUT -T $THREADS --ref-overlap -m ga4gh -f "QUAL" -c $VCF -b $TRUTH --all-records

python qfy.py -o $OUTPUT --write-vcf --write-counts --verbose -t ga4gh --stratification $STRATS --roc- delta 20 --roc $QQ -r $FASTA $VCF

And still got some nonzero FP.AL values. It is a shame as I am trying to do the benchmarking for a diagnostics pipeline aso it is quite critical to understand the expected behaviour of happy/vcfeval. Also I am open to try other tools but as far as I know hap.py+vcfeval is the goldstandard toolset for benchmarking NGS pipelines (small variants). If you have any suggestion on what to try next please let me know! Also, thanks so much for opening an issue in hap.py repo but I am not very hopeful though as the hap.py github repository has not been updated in ages right? ANyway, thanks so much for looking into this I will keep investigating and let you know if I find out what is going on.

Thank you Laura for letting me know! As far as I know, yes happy+rtgtools is best practice. But you can also experiment the other tools (happy or rtgtools - standalone) to see the differences. Best.

Hi @Laura Mora Bitria, I'm also getting "Error: An IO problem occurred: “Invalid file pointer" errors when running rtg vcfeval with my own vcf. Were you able to find out more or fix this?

Hello @Felix Lenner, I would suggest to validate the input vcf files, and do you generate sdf files yourself or is it automatic? Also, if you are resuming your runs (with some errors) you might need to regenerate intermediate files (so maybe restart your run without --resume)

Hi, do you have a good suggestion for validating the input files? I ran rtgtools standalone and generated the sdf files myself.

I see, what I was suspecting, corrupted index or sdf file.

Ok, there's a mismatch between references, where I use

GRCh38_GIABv3_no_alt_analysis_set_maskedGRC_decoys_MAP2K3_KMT2C_KCNJ18.fasta

and the truth set uses

human_GRCh38_no_alt_analysis_set.fasta

. I only tried using the GIAB version as sdf, so I could try with the other one as well.

if you are only interested in main chromosomes, you can filter out decoys.

Probably my VCF that is a problem for rtg. Tried an unannotated version of it which works fine.

is it long read variant calling?

Yes!

I think I had similar error once trying..

Actually, now that running rtgtools outside variantbenchworking with my VCFs works, when running them with variantbenchmarking I get the same error, but on the truth set VCF instead. Any idea? 😅

Error: An IO problem occurred: "Invalid file pointer: 1825570489 for HG002_GRCh38_1_22_v4.2.1_benchmark.reheader.vcf.gz"

is it the same truth you were using?

Yeah, is the only thing happening to the truth within the pipeline reheadering?

if you are not using any preprocesses - yes