Clarification on VCF Output and Merging in PHG Workflow

[hyoj-seo]

Hello,
I am currently trying out your PHG software.

Thanks to your detailed tutorial, I was able to complete the imputation step successfully. However, I have a few questions before proceeding further.

(1) Two h.vcf files were generated after mapping with paired-end fastq files.

Here is the keyfile I used for phg map-reads:

sampleName      filename                            filename2  
2022-1          2022/1_S29_L001_R1_001.fastq.gz     2022/1_S29_L001_R2_001.fastq.gz  
2022-102        2022/102_S84_L003_R1_001.fastq.gz   2022/102_S84_L003_R2_001.fastq.gz

The mapping process completed successfully, but unlike the tutorial, I obtained the following files:

output/read_mappings/2022-1_1_readMapping.txt  
output/read_mappings/2022-1_2_readMapping.txt  
output/read_mappings/2022-102_1_readMapping.txt  
output/read_mappings/2022-102_2_readMapping.txt

For phg find-paths, I created a keyfile like this:

sampleName      filename  
2022-1          output/read_mappings/2022-1_1_readMapping.txt  
2022-1          output/read_mappings/2022-1_2_readMapping.txt  
2022-102        output/read_mappings/2022-102_1_readMapping.txt  
2022-102        output/read_mappings/2022-102_2_readMapping.txt

This resulted in one VCF file per sample.
Does this mean the VCFs were generated correctly?

(2) Is it okay to merge multiple imputed VCFs and use them with resequencing data?

I performed imputation using multiple samples, so I now have one h.vcf per sample.
I would like to combine these and use the merged result as a reference for future resequencing steps.

I found a command called merge-hvcfs in PHG.
Would it be appropriate to use this to generate a single merged.h.vcf file from the imputed VCFs?

Thank you in advance for your help.

Best regards,
HJ

[zrm22]

For 1: The reason why there are two read mapping files is that we use ropebwt3 and it is not paired end aware. So we align each file in the pair and put them in their own readMapping file. Then in the path keyfile you can merge them back together by using the same sampleName for the two read mappings. It looks like you have created your hVCFs correctly in this case.

2: merge-hvcfs should combine the multiple samples together which you could use for downstream tasks. I think that is the correct command for what it sounds like you are trying to do.

Thanks,
Zack

[hyoj-seo]

@zrm22
Thank you for your response!

Following your instructions, I attempted to create the composite FASTA file.

However, I encountered an issue where certain haplotype IDs present only in the merged VCF do not exist in the vcf_dbs. To resolve this, I tried to run phg load-vcf.

During that process, I received the following error:
[critical] Exception: Combined VCFs are currently not supported

According to the tutorial, it seems that I cannot proceed beyond this step without resolving the issue.
Is there any possible workaround?

Best regards,
HJ

(+)

I did some additional testing.
I tried running the command without specifying the DB path, but encountered an issue that forced me to specify the DB folder name.

[main] INFO net.maizegenetics.phgv2.cli.CliktLogging 2025-06-24 09:34:09,424: PHGv2 Command: create-fasta-from-hvcf --hvcf-file output/vcf_files_imputed_merged/imputed.h.vcf --fasta-type composite -o output/composite_assemblies/
    [--db-path] =
    [-o, --output-dir] = output/composite_assemblies/
    [--fasta-type] = composite
    [--hvcf-file] = HvcfFile(hvcfFile=output/vcf_files_imputed_merged/imputed.h.vcf)
    [--hvcf-dir] = null
    [--range-bedfile] = null
    [--conda-env-prefix] =
[main] INFO net.maizegenetics.phgv2.utils.VariantLoadingUtils 2025-06-24 09:34:09,428: TileDB datasets not found in folder /data3/03.GBS/02.2025_May_radish_PHG/01.PHG. Either send a valid dbPath folder variable, or run InitDB to create the datasets in the specified folder.
Exception in thread "main" java.lang.IllegalArgumentException: TileDB datasets not found in folder /data3/03.GBS/02.2025_May_radish_PHG/01.PHG. Please run InitDB to create the datasets.

[zrm22]

I am unsure of how you can get haplotypeIds which show up in your index without being in your vcfDB unless they were not loaded in prior to running the imputation. load-vcf is currently only designed to load in a single samples gvcf/hvcf files and as such is why you are seeing this error. We probably could work something out but it will likely take some time.

My suggestion would be to build a new tile-db instance then make sure all the samples you used in the index for imputation are loaded in.

[hyoj-seo]

@zrm22
Thank you for your reply.

As you recommended, I loaded all the imputed .h.vcf files into the tile-db using the following command:

phg load-vcf --vcf-dir output/vcf_files_imputed/gzip --db-path vcf_dbs --threads 196

After loading, I confirmed that {imputed_sample_name}.h.vcf.gz files were created inside ./vcf_dbs/hvcf_files/.
(Previously, only Ref.h.vcf and assembly.h.vcf existed.)

Then, I ran:

phg create-fasta-from-hvcf --db-path vcf_dbs --hvcf-file output/vcf_files_imputed_merged/gzip/imputed.h.vcf.gz --fasta-type composite -o output/composite_assemblies/

However, I encountered the following error:

[main] INFO net.maizegenetics.phgv2.cli.CliktLogging 2025-06-25 09:30:34,237: PHGv2 Command: create-fasta-from-hvcf --db-path vcf_dbs --hvcf-file output/vcf_files_imputed_merged/gzip/imputed.h.vcf.gz --fasta-type composite -o output/composite_assemblies/
    [--db-path] = vcf_dbs
    [-o, --output-dir] = output/composite_assemblies/
    [--fasta-type] = composite
    [--hvcf-file] = HvcfFile(hvcfFile=output/vcf_files_imputed_merged/gzip/imputed.h.vcf.gz)
    [--hvcf-dir] = null
    [--range-bedfile] = null
    [--conda-env-prefix] =
[main] INFO net.maizegenetics.phgv2.utils.VariantLoadingUtils 2025-06-25 09:30:34,242: begin Command:conda run -n phgv2-tiledb tiledbvcf stat --uri vcf_dbs/hvcf_dataset
[main] INFO net.maizegenetics.phgv2.utils.VariantLoadingUtils 2025-06-25 09:30:37,751: Using  TileDB datasets created in folder vcf_dbs.
Exception in thread "main" java.lang.IllegalStateException: Haplotype ID R000000_2022-1_G0 not found in ALT Header
        at net.maizegenetics.phgv2.cli.CreateFastaFromHvcf.createHaplotypeSequences(CreateFastaFromHvcf.kt:307)
        at net.maizegenetics.phgv2.cli.CreateFastaFromHvcf.processSingleHVCF(CreateFastaFromHvcf.kt:256)
        at net.maizegenetics.phgv2.cli.CreateFastaFromHvcf.buildFastaFromHVCF(CreateFastaFromHvcf.kt:136)
        at net.maizegenetics.phgv2.cli.CreateFastaFromHvcf.run(CreateFastaFromHvcf.kt:505)
        at com.github.ajalt.clikt.parsers.Parser.parse(Parser.kt:306)
        at com.github.ajalt.clikt.parsers.Parser.parse(Parser.kt:319)
        at com.github.ajalt.clikt.parsers.Parser.parse(Parser.kt:40)
        at com.github.ajalt.clikt.core.CliktCommand.parse(CliktCommand.kt:458)
        at com.github.ajalt.clikt.core.CliktCommand.parse$default(CliktCommand.kt:455)
        at com.github.ajalt.clikt.core.CliktCommand.main(CliktCommand.kt:475)
        at com.github.ajalt.clikt.core.CliktCommand.main(CliktCommand.kt:482)
        at net.maizegenetics.phgv2.cli.PhgKt.main(Phg.kt:43)

This haplotype ID R000000_2022-1_G0 appears only in the merged imputed VCF file:

#zgrep "R000000_2022-1_G0" *.h.vcf.gz
NC_018551.1     7923    .       A       <R000000_2022-1_G0>,<R000000_2022-158_G1>       .       .       END=212617      GT      1|1     1|1     1|1     1|1     1|1     1|1     1|1     1|1     1|1     1|1  1|1      1|1     1|1     1|1     1|1     1|1     1|1     1|1     1|1     1|1     1|1     1|2     1|2     1|1     1|1     1|1     1|1

In other single-sample imputed .h.vcf files, the same variant position looks like this:

#grep -P "NC_018551.1\t7923" 2023-197.h.vcf
NC_018551.1     7923    .       A       <d9c6e73ae940173fd174f9eb89387e9f>,<2cb4b159491f4ae722b08967a55072c6>   .       .       END=212617      GT      1|2

I suspect that phg merge-hvcfs is rewriting haplotype IDs during merging—am I correct?

Additionally, haplotype IDs like <d9c6e73ae940173fd174f9eb89387e9f> are defined in the header as:

##ALT=<ID=d9c6e73ae940173fd174f9eb89387e9f,Description=...>

Whereas <R000000_2022-1_G0> is not found in the header, or it appears mismatched like:

##ALT=<ID=R000000_RS05_G0,Description=...>

Regardless of these inconsistencies,
am I correct in understanding that if multi-sample load-vcf is not supported, then using a merged imputed VCF file as a reference for resequencing variant calling is not currently feasible?

Thank you for making your program available to use. I really appreciate it.