High chrY call rates for human female 30x WGS samples, what am I doing wrong?

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

gatk plink chromosomes sex callrate

• 257 views • link

updated 1 month ago by SeqBench 41 • | written 4 months ago by Joel Wallenius 22

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.

• 0 views • link

written 2 months ago by SeqBench 41

Thank you very much for pointing those two things out, I was steering by AI which evidently fails spectacularly sometimes. Perhaps I kept the --hwe and --maf options because I wanted to use both Y call rates and X homozygosity, without having to run PLINK once for each with different filters. I don't remember, honestly.

I switched to seGMM, it computes depth in various informative locations, ratios of aforementioned, and X homozygosity. Only the contaminated samples lack a clear sex distinction, with seGMM.

• 0 views • link

written 1 month ago by Joel Wallenius 22

Glad it sorted out. The contaminated ones staying ambiguous is expected rather than a seGMM failure - a female sample with male DNA in it genuinely sits between the two clusters, so the intermediate value is the contamination showing up rather than the sex check breaking. If you have a contamination estimate per sample (verifybamid or similar) it's worth plotting that against the seGMM score, they usually line up.

• 0 views • link

written 1 month ago by SeqBench 41