Hi,

First off, that pattern in your Manhattan plot --gaps around the centromeres-- is pretty common and usually not a red flag for your analysis. Technically, it's reasonable because centromeric regions are notorious for being tough nuts to crack in sequencing: they're highly repetitive, have weird GC content, and reads map poorly there (or not at all). Your hard filters (especially MQ < 60 and the rank sum tests) are spot-on for GATK best practices, but they naturally cull SNPs in those areas since mapping quality tanks and strand bias creeps in from the repeats. Biologically, it makes sense too—recombination is suppressed near centromeres, so polymorphism is often lower anyway, especially if your species has that heterochromatin setup (assuming this isn't human; if it is, same deal).

Your summary stats look solid overall: median MAF at 0.09 is decent for discovery, missingness under 0.25 is fine, and those QD/FS values aren't screaming problems. The mean FS at 57.63 jumps out a bit (even with your >5 filter), but if it's driven by a few outliers post-filter, it might just be noise from dense regions elsewhere.

To double-check reliability, here's a stepwise plan—focus on diagnostics first, then poke at the data:

1. QQ plot and genomic control lambda: This is your quick gut-check for model fit. In rMVP, after x<-MVP() , grab the p-values and plot:

   qqplot(-log10(1:p), -log10(gwas$p.values), main="QQ Plot")
   abline(0,1)
   lambda <- median(qchisq(1-gwas$p.values,1), na.rm=T) / qchisq(0.5,1)
   

   Lambda around 1-1.05 is golden; >1.1 might mean unaccounted structure (try adding more PCs if you have them).

2. SNP density by chromosomal position: Slice your VCF or genofile by position and histogram it per chromosome to visualize those gaps. Load via vcfR (it's stable, no updates needed):

   library(vcfR)
   vcf <- read.vcfR("your_filtered.vcf.gz")
   pos <- getFIX(vcf)$POS
   chrom <- getCHROM(vcf)
   # For chr1 example
   hist(pos[chrom=="1"], breaks=1000, main="SNP Density Chr1")
   # Overlay known centromere coords if you have 'em (e.g., from UCSC)
   

   If gaps align with centromeres (grab coords from your genome browser), you're good—it's not your filters eating signal.

3. Mappability track: Pull a mappability score (e.g., from ENCODE or your aligner logs) and correlate it with your SNP positions. Low mappability (<0.8) should match your gaps. If not, loosen MQ to 50 and re-run filtering to see if centromeric SNPs pop up (but watch for false positives).

4. Cross-validate the model: Re-run MLM without kinship (if you're using it) or try FarmCPU in rMVP for comparison—MVP.Config(method="FarmCPU"). If peaks shift but centromere gaps stay, it's technical. Also, subsample your SNPs (e.g., LD-prune harder) and re-plot to rule out over-pruning.

5. External benchmark: Compare your SNP set to a public panel for your species (e.g., dbSNP or a ref GWAS). Overlap and see if your missing centromeric hits are also absent there.

If lambda's off or densities look wonky beyond centromeres, dig into population structure (PCA on your genotypes) or batch effects. Otherwise, this is likely just the data behaving as expected—plenty of published plant GWAS (guessing from rMVP?) show those barren zones.

Kevin

I don't know how this is usually handled in GWAS studies, but the peri-centromeric regions of chromosomes (as well as the short arms of acrocentric chromosomes) are highly repetative. The centromeric regions themsevles are simple arrays of microsatelites, and as such, it is generally not possible to align reads or call SNPs in these regions.

What I don't really know of the top of my head is how large you'd expect these uncalled regions to be.