Hello,
I'm pulling my hair! For days now I've been trying different filters so that I can clearly differentiate female from male samples in a histogram of chrY call rates. The best I can do is end up with males all having a chrY call rate of exactly 1.0, while the females are dispersed between 0.8 and 0.99.
I tried adding filters in PLINK2 like --min-vcf-GQ and --min-vcf-DP, but that just makes it worse because the males instead get call rates around 0.98.
The pipeline is roughly:
GATK calls non-PAR regions that are mappable:
chrY:1-10000
chrY:2781480-56887902
chrY:57217416-57227415
but excluding the unmappable bits:
chrY:4343800-4345800
chrY:10246200-11041200
chrY:11072100-11335300
chrY:11486600-11757800
chrY:26637300-57227400
BCFtools sets all calls with DP=0 to missing.
PLINK2 imports and ultimately does --impute-sex.
$plink2 --vcf "$vcf"\
--psam "input/samples_with_sex.psam"\
--set-all-var-ids '@_#_$r_$a'\
--make-pgen\
--snps-only\
--max-alleles 2\
--hwe 1e-6\
--maf 0.01\
--geno 0.05\
--threads 48\
--memory 240000\
--out output/sexcheck/pfiles/sexcheck_MP1
There are 1233 variants left after these filters. If I add the filters I mentioned above, the results are worse.
What am I missing here? Certainly 114 WGS samples can't all be poorly sequenced or something... I used a no-ALTs reference of hg38.
Sincerely, Joel
•
written
by
Joel Wallenius
1 answer
Two problems, and the second one is doing most of the damage. Call rate on chrY genotypes is a weak statistic to start with, because females still pick up reads across non-PAR chrY from the X-transposed and ampliconic regions, so masking DP=0 won't clear them.
The bigger issue is that your filters are actively selecting for exactly those sites. --hwe is meaningless on a haploid chromosome, and --maf 0.01 across a mixed-sex cohort keeps precisely the positions where females carry non-reference calls, which are the X-Y mismapping artefacts. So your surviving 1233 variants are close to the worst possible set for telling the sexes apart.
I'd drop genotypes for this and use depth instead. Mean coverage over the non-PAR chrY intervals you already defined, normalised to the autosomes, comes out around 0.5 for males and near zero for females, and it doesn't care about mismapping the way variant calls do. mosdepth over those intervals is a few minutes of work. If you want a genotype-based cross-check, chrX heterozygosity is far more trustworthy than anything on Y.
Log in to answer this question.