Entering edit mode
8 months ago
b.contreras.moreira
▴
310
Hi, after successfully running
phg create-ref-vcf \
--bed /data/phg/data/reference/ref_ranges.bed \
--reference-file ... \
--reference-name ... \
--db-path /data/phg
The next step failed, can't see any other log message that I can post here:
(phgv2-conda) $ phg create-maf-vcf --db-path /data/phg --bed /data/phg/data/reference/ref_ranges.bed --reference-file ...fa --maf-dir /data/phg/data/alignment_files -o /data/phg/data/vcf_files
[main] INFO net.maizegenetics.phgv2.utils.VariantLoadingUtils 2024-03-19 13:35:51,445: begin command:conda run -n phgv2-conda tiledbvcf stat --uri /data/phg/hvcf_dataset
[main] INFO net.maizegenetics.phgv2.utils.VariantLoadingUtils 2024-03-19 13:35:52,960: Using TileDB datasets created in folder /data/phg.
[main] INFO net.maizegenetics.phgv2.cli.CreateMafVcf 2024-03-19 13:35:53,370: CreateASMHvcfs: calling buildRefGenomeSeq
[main] INFO net.maizegenetics.phgv2.cli.CreateMafVcf 2024-03-19 13:36:38,240: CreateASMHvcfs: processing /data/phg/data/alignment_files/HOR13942.maf
getVariantContextsfromMAF: Processing a single genome
Exception in thread "main" java.lang.NullPointerException
at biokotlin.genome.MAFToGVCF.buildRefBlockVariantInfo(MAFToGVCF.kt:990)
at biokotlin.genome.MAFToGVCF.buildTempVariants(MAFToGVCF.kt:542)
at biokotlin.genome.MAFToGVCF.buildVariantsForAllAlignments(MAFToGVCF.kt:246)
at biokotlin.genome.MAFToGVCF.getVariantContextsfromMAF(MAFToGVCF.kt:167)
at biokotlin.genome.MAFToGVCF.getVariantContextsfromMAF$default(MAFToGVCF.kt:136)
at net.maizegenetics.phgv2.cli.CreateMafVcf.getGVCFVariantsFromMafFile(CreateMafVcf.kt:142)
at net.maizegenetics.phgv2.cli.CreateMafVcf.getGVCFVariantsFromMafFile$default(CreateMafVcf.kt:141)
at net.maizegenetics.phgv2.cli.CreateMafVcf.createASMHvcfs(CreateMafVcf.kt:80)
at net.maizegenetics.phgv2.cli.CreateMafVcf.createASMHvcfs$default(CreateMafVcf.kt:67)
at net.maizegenetics.phgv2.cli.CreateMafVcf.run(CreateMafVcf.kt:611)
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:45)
Any help appreciated, Bruno
WIll you check your reference fasta and the maf file to verify the chromosome names for the MAF reference matches the chromosome names used in the reference fasta?
Please see below lcj34 :
And:
Only thing I can see is that unplaced contigs in the reference (CAJHDD...) do not appear in the MAF file
In the second listing - I see they are sorted, so I cannot tell if the ref is chr1H or chr1H_LR890096.1. If the former, this would be the problem. (inconsistent ref naming). If the latter, then something else is going on . If the ref contig names are the same, the send me privately your reference file and MAF file and I'll try to debug.
Hi the MAF file in question looks like this:
I will email you the tarred files for your debugging, thanks for the offer, Bruno
I am re-running leaving out unplaced contigs to see if that fixes it, will let you know
Bruno - that won't fix it. I found a problem with the BioKotlin library that phgv2 calls. I have fixed the library (also Buckler Lab code) and tested with your data. Once my fix is approved I will create a new jar file, then we'll rebuild phgv2 with this new code and make that available.
On a related note, the latest version of the phgv2 documents (as of today) are clearer on the requirements of running the annotate-fastas command. This is necessary to keep the AGC compressed assemblies (sequence) and tiledb datasets (variants) in sync with the same sample names. I don't know if you have run that step, but it shoudl be run prior to running agc-compress or align-assemblies. See the documents for further details.
BTW - we appreciate you testing phgv2 and your patience as we work out issues.
As agreed by email, I have now updated the phg code to version 2.2.65.114. Now I can go past that stage, only to fail with bcftools index:
When I run that last command manualy it says:
You will need to sort the file before indexing.
Thanks GenoMax , I would have thought phg took care of that, right lcj34 ?
Edit: I tested sort followed by index and it worked, so I guess whatever the code is trying to index needs to be properly sorted first?
Yes, the gvcf created by phg should have been sorted. I am able to reproduce the problem using your file. Will look into this and let you know.
Bruno - I have fixed this issue in our library. The fix won't be available in phgv2 for a couple days as the person who creates our BioKotlin releases (phgv2 uses this library) is out early this week. It should be available by the end of the week - I'll let you know.
Thanks lcj34
FYI - phgv2 version Release 2.2.68.117 has the updated BioKotlin release with this fix.
Hi lcj34 , I downloaded the latest version and got a different error now:
Thanks, Bruno. I found a place with a missing sort call (working when BioKotlin was called directly, but when called from phgv2 a sort call was missing). I've fixed that - will let you know when the code has been reviewed and merge.
This fix is now available in phgv2 version 2.2.72.121
I confirm this now works with my test data, many thanks lcj34 !