We have 4764 30x whole genome sequencing samples and I am attempting to import the individually genotyped files (.g.vcf.gz) with GLnexus on our hpc. These were created from GATK haplotypecaller. Each one is ~7 GB big. Just for chr21, the RocksDB building took a whopping 100 TB disk space and died after 10 days, the time limit on the cluster partition I used. It was about to finish. Offered 250 GB memory and that was fine, it used 180 GB. In a different job, I tried giving 3 TB memory and it only used ~200 GB again. It does not seem to be a memory issue. I think it is the sheer number of variants and samples that are inflating the db. Even if I try running it on our unlimited-time partition, it will take >10 days for chr21 only, and due to the disk space I could not run other contigs at the same time. Also, using local scratch did not help, it uses all 3 TB scratch available and dies. So I have been writing to regular disk space, which may be slowing it down.
Our overall goal is to run QTLtools and colocalization on all common noncoding variants and a metabolomics dataset to find variants and metabolites associated with our disease severity.
My colleagues and I are actually stumped what to try next. I also submitted an issue on the GLnexus github.
•
written
by
Claire Weber •
call the GVCFs files into chuncks of , say, 1Mb (overlap 100bp) , merge each chunk of 1MB with GLNexus, and at the end, 'bcftools concat' all the 1Mb vcfs.
EDIT: graphtyper is also good at calling thousand of genomes (by chuncks of 50kb)
GLNexus isn't a variant warehouse - it's just a joint genotyping tool. Since it's not even clear if you need joint genotyping (in the sense of changing genotype calls based on recurrent artefacts across samples) - it might be a waste of time anyway. It sounds like you just want a big genotype matrix.
See my discussion of variant warehouses below (I've updated it for you): Is there a file format better suited for the era of pangenomics than the .vcf? What are its attributes?
(Disclaimer: I am product manager for TileDB-VCF)
I am trying to get jointly called per-contig files. I think we do need to do joint calling for QTLtools. We plan to do it on common noncoding variants only MAF > 0.1. thank you
if we agree on what joint calling means in the gatk sense (removing artefacts) then wouldn't that be not so important for common variants? those probably aren't artefacts - they're just common.
Does MAF refer to global MAF or cohort MAF? If the former in theory you could use that information to filter your gVCFs as a first line (although I don't know how to filter gVCFs reference ranges and I do this for a living). Maybe convert them to VCFs first then filter, then get your matrix from PLINK.
that's a good point. well we are interested in the EUR MAF rarity. I did consider filtering first but not sure if it would make a difference for computational load. to update this post I am currently trying GATK import & joint calling.
i'll go on the record as saying i don't think you need joint calling for common variants. just get the 0/1/2 values for GT at common sites (from gnomAD) into a matrix and get on with your life.
I forgot to mention that QTLtools is not concerned with indels.
Guess what the main "rescued" variant type is for joint calling? Rare indels.
Log in to answer this question.
GitHub for
GLnexusmakes the above claim so in theory your dataset should be smaller than this though the statement above does not include any clue about what kind of hardware was used (the paper also does not seem to include this information). Considering that the last project update was in 2021,GLnexusmay either be stable software (or is no longer being actively maintained). Since many of the recent issues have not received responses, the latter may be the case.Exploring alternate software options may be the way to go.
True, it could be a software issue for this scale. What do you recommend then, just GATK for joint calling? thank you