GWAS Analysis in R: A Working Guide

By the end of this guide you will have a quality-controlled genotype matrix on disk and principal components that capture population structure. You will also have a table of per-variant association statistics with a genomic inflation factor you can defend, a Manhattan and QQ plot, and a polygenic score computed on your own genome from published summary statistics.

The requirements are modest by genomics standards: R 4.3+, PLINK 2.0, bcftools, roughly 16 GB of RAM and 100 GB of disk for a cohort of a few thousand samples, and the R packages bigsnpr, bigstatsr, data.table, and ggplot2.

One point is worth settling before you start. A genome-wide association study is a cohort method, and with a single sample you cannot estimate an association at all. What you can do with one genome is consume other people’s GWAS results, and that is the subject of the second half of this guide.

1. Decide which of the two jobs you are doing

People use the phrase “GWAS analysis in R” for two workflows that share tooling but answer different questions. Knowing which one you are doing determines everything that follows, so it is the first decision to make.

The first workflow is a discovery GWAS. You have genotypes and a phenotype for thousands of people, and you regress the phenotype on each variant separately. Typical scale is 500k to 1M directly genotyped variants on an array, or 8 to 15 million variants after imputation or from whole-genome sequencing once you filter to minor allele frequency above 0.5 percent. A recent example gives a sense of the shape of a real study: whole-genome sequencing of 6,218 Qatari individuals tested against 45 clinically relevant traits, which recovered known loci and also found population-specific ones ¹.

The second workflow applies existing results to a genome you already have. You take public summary statistics, match them to your variants, and compute a polygenic score or look up individual hits. This is the workflow that applies to a Personal Molecular Profile, and it is where most of the interpretive traps live.

Because both workflows use the same file formats and much of the same R code, we will build the machinery once and use it twice. If you want to practice the discovery half without a cohort of your own, use 1000 Genomes Phase 3 with a simulated phenotype. It is the fastest way to learn the failure modes without waiting on real data.

2. Get from VCF to PLINK binary format

Your starting point is almost always a VCF, the Variant Call Format text file that lists each variant and each sample’s genotype at it. The goal of this section is to convert that into PLINK’s compact binary format, which is what every downstream step expects.

Resist the temptation to do quality control in R. PLINK 2 performs the same operations two orders of magnitude faster, and its flags are the de facto standard, so anyone can reproduce your work from the command line alone. R comes in once you have a clean matrix.

Normalize the VCF first. Multiallelic sites and indels that have not been left-aligned are the most common source of silent allele mismatches downstream, and normalization removes both problems before they can propagate.

bcftools norm -m -any -f GRCh38.fa cohort.vcf.gz -Ou \
  | bcftools norm -d exact -Oz -o cohort.norm.vcf.gz
bcftools index cohort.norm.vcf.gz

plink2 --vcf cohort.norm.vcf.gz \
  --max-alleles 2 --snps-only just-acgt \
  --set-all-var-ids '@:#:$r:$a' --new-id-max-allele-len 60 truncate \
  --rm-dup exclude-all \
  --make-bed --out cohort.raw

The flag --set-all-var-ids '@:#:$r:$a' gives every variant a chromosome:position:ref:alt identifier. This matters because rsIDs are unstable across dbSNP builds and missing entirely for many sequenced variants, whereas position-based IDs let you join tables later without consulting a lookup service. Record your genome build now and write it into the filename, since a build mismatch between your genotypes and downstream summary statistics is the single most common failure in this whole pipeline.

3. Sample and variant QC

Quality control decides which samples and which variants you can trust. The core filters run in a single PLINK command:

plink2 --bfile cohort.raw \
  --geno 0.02 --mind 0.02 --maf 0.01 --hwe 1e-10 midp \
  --make-bed --out cohort.qc

Order matters here in a way that is easy to miss. PLINK applies variant filters before sample filters within a single run, so if you have a few catastrophically bad samples, run --mind alone first and then re-run the variant filters on the survivors. Otherwise those bad samples inflate missingness at every variant and you discard good sites along with the bad.

There are two further checks that PLINK will not perform unless you ask for them, and both catch problems the standard filters miss:

# heterozygosity outliers (contamination or sample mixture)
plink2 --bfile cohort.qc --het --out cohort
# relatedness
plink2 --bfile cohort.qc --king-cutoff 0.0884 --out cohort.unrel

A KING coefficient above 0.0884 corresponds roughly to third-degree relatives or closer. Either remove one member of each such pair or move to a mixed model that handles relatedness explicitly. Running ordinary least squares on a sample with cryptic relatedness and hoping the principal components absorb it is not a workable strategy.

With a clean dataset in hand, load it into R once as a file-backed object:

library(bigsnpr)
snp_readBed("cohort.qc.bed")            # writes cohort.qc.rds + .bk
obj <- snp_attach("cohort.qc.rds")
G   <- obj$genotypes                     # FBM.code256, not in RAM
CHR <- obj$map$chromosome
POS <- obj$map$physical.pos
dim(G)   # samples x variants

bigsnpr keeps the genotype matrix memory-mapped in a .bk backing file, which means a 5,000 x 8,000,000 matrix costs you 40 GB of disk and almost no RAM. That property is why we prefer it to snpStats or GWASTools for anything at sequencing scale.

4. Population structure

Ancestry stratification is the classic way to produce associations that are real, reproducible, and completely spurious. If allele frequencies and phenotype means both vary with ancestry, every variant that differs between groups will look associated. The standard correction is to compute principal components on LD-pruned variants, with long-range LD regions removed first.

excl <- snp_indLRLDR(infos.chr = CHR, infos.pos = POS)
svd  <- snp_autoSVD(G, CHR, POS, ind.excl = excl,
                    k = 20, ncores = nb_cores())
plot(svd, type = "scores", scores = 1:6)

snp_autoSVD works iteratively. It prunes on linkage disequilibrium and runs a truncated singular value decomposition. It then detects variants with outlying loadings, usually an unflagged inversion or a leftover LD block, removes them and repeats. Inspect both the scree plot and the loadings: if PC7 is driven by a few dozen variants in one 2 Mb window, that is an artifact rather than ancestry.

Keep as many PCs as show genuine structure, which is typically 10 to 20. In an ancestrally homogeneous cohort, 4 is often enough. If your cohort is genuinely admixed, PCs alone will not fully fix the problem and you should move to a linear mixed model. Use either GENESIS::assocTestSingle in R or GCTA, REGENIE, or SAIGE outside it.

5. Run the association test

Now the actual regression. For a quantitative trait with covariates, the call is short because all the heavy lifting happens on the file-backed matrix:

y <- obj$fam$affection            # or read your phenotype file
covar <- cbind(svd$u[, 1:10], age, sex)

gwas <- big_univLinReg(G, y.train = y,
                       covar.train = covar,
                       ncores = nb_cores())
head(gwas)   # estimate, std.err, score (t-stat)

For a binary trait, big_univLogReg does the same thing with logistic regression. It is slower, and it will fail to converge on rare variants in unbalanced case-control designs. If your case fraction is below about 5 percent, or you are testing singleton-level rare variants, use SAIGE or REGENIE instead of a plain logistic model. Score-test-based methods with saddlepoint approximation exist precisely because the Wald test breaks down in that regime.

Before you look at any hits, look at inflation. The genomic inflation factor λ_GC compares your observed test statistics to what you would expect under the null, and it is the quickest way to find out whether your results are trustworthy:

snp_qq(gwas)                  # prints lambdaGC in the subtitle
ggplot2::last_plot()

A λ_GC near 1.00 to 1.05 is fine. A λ_GC of 1.2 in a study of 3,000 people means you have residual structure, relatedness, or a batch effect. Every peak on your Manhattan plot is then suspect. λ_GC also rises legitimately with sample size for highly polygenic traits, which is why the LD score regression intercept becomes the better diagnostic once you are above roughly 20,000 samples.

With inflation under control, plot the results:

snp_manhattan(gwas, CHR, POS, npoints = 50000) +
  ggplot2::geom_hline(yintercept = -log10(5e-8), linetype = 2)

The 5×10⁻⁸ threshold is 0.05 divided by roughly one million independent common-variant tests in European-ancestry LD structure. It is a convention rather than a law. Dense whole-genome sequencing data in African-ancestry cohorts contains more independent tests and justifies a stricter line.

6. Clump, then interpret conservatively

Genome-wide significant variants arrive in correlated blocks, so a raw hit list overstates how many findings you have. Clumping reduces those blocks to independent signals:

plink2 --bfile cohort.qc --clump gwas.tsv \
  --clump-p1 5e-8 --clump-r2 0.1 --clump-kb 1000 \
  --out gwas.clumped

Interpretation then needs care in three respects. The lead variant is almost never the causal variant; it is the best-tagged marker in an LD block. The nearest gene is a guess. And effect sizes for common variants on complex traits are small for structural reasons: selection keeps large-effect alleles rare, so the common variants that GWAS is powered to find sit at the low end of the effect distribution ². Fifteen years of this work has produced tens of thousands of trait-associated loci, with the bulk of them in noncoding regulatory sequence rather than protein-coding exons ³.

If you want to move from a locus to a mechanism, the productive next step is functional data rather than more association testing. Matching GWAS signals to expression QTLs across tissues is how a substantial fraction of noncoding associations get assigned to a gene, and whole-genome sequencing combined with RNA-seq in the same individuals improves causal variant identification over imputed genotypes alone ⁴. If you have your own RNA-seq, this is the join that makes a locus mean something.

A coarser but still informative view comes from aggregating to genes or pathways. sumFREGAT runs gene-based tests directly on summary statistics plus an LD reference, which is useful when you cannot share individual-level genotypes ⁵. Pathway-level analysis on GWAS or WGS output brings its own set of method choices and multiple-testing corrections, which are worth reading about before you run one ⁶.

7. Apply published summary statistics to your own genome

This is the part that works with a single sample. The approach is to download summary statistics from the GWAS Catalog in the harmonized format and then use LDpred2 in bigsnpr to turn them into weights you can apply to your own genotypes.

library(bigsnpr)
sumstats <- bigreadr::fread2("GCST90000001_buildGRCh37.tsv.gz")
names(sumstats) <- c("chr","pos","a1","a0","freq","beta","beta_se","p","n_eff")

# your own genotypes, same build, HapMap3+ variants
map <- obj$map[, c("chromosome","physical.pos","allele1","allele2")]
names(map) <- c("chr","pos","a0","a1")

info_snp <- snp_match(sumstats, map, strand_flip = TRUE)

snp_match handles reverse-complement and allele-swap cases and reports how many variants survived the match. If it keeps less than about 70 percent of the HapMap3+ set, you have a build mismatch or a column ordering error, so check that before proceeding.

The next step is to estimate per-chromosome LD from a reference panel, run LD score regression to get a heritability estimate, and fit the model:

ldsc <- snp_ldsc2(corr, df_beta)     # h2 estimate + intercept
h2   <- ldsc[["h2"]]
auto <- snp_ldpred2_auto(corr, df_beta, h2_init = h2,
                         vec_p_init = seq_log(1e-4, 0.2, 30),
                         ncores = nb_cores())
beta_auto <- rowMeans(sapply(auto, function(a) a$beta_est))
pred <- big_prodVec(G, beta_auto, ind.col = info_snp$`_NUM_ID_`)

snp_ldpred2_auto learns the polygenicity parameter from the data rather than requiring a validation cohort, which is exactly what you want when you have one sample. Drop any chains whose scale is an outlier relative to the median before averaging.

Three things are worth understanding about the number that comes out the other end. First, a raw polygenic risk score is meaningless in isolation. It only means something as a percentile within an ancestry-matched reference distribution, so score a few thousand 1000 Genomes samples from the matched superpopulation with the same weights and place yourself in that distribution. Second, scores derived in European-ancestry cohorts lose a large fraction of their predictive accuracy when applied to other ancestries, because LD patterns and allele frequencies differ; replication and portability have been the persistent weak points of the field ⁷. Third, a high percentile is a statement about a population distribution rather than about you. It does not diagnose anything, it does not tell you what will happen, and any question about what to do with the result belongs with a clinician or a genetic counselor rather than with your R session.

Common problems

The failures in this pipeline tend to be quiet rather than loud, producing plausible output from broken inputs. The list below covers the ones you are most likely to meet, with what to do about each.

Genome build mismatch. GRCh37 and GRCh38 coordinates differ for essentially every variant. snp_match will return a plausible-looking but tiny overlap, or worse, a moderate overlap built from coincidences. Always confirm the build of both inputs from the file header rather than the filename. Lift over with liftOver or CrossMap if needed, and accept that a few percent of variants will fail to map.

Palindromic SNPs. A/T and C/G variants cannot be strand-resolved by allele codes alone. snp_match(strand_flip = TRUE) removes ambiguous ones by default, and you should not override this to raise your overlap count. Frequency-based resolution only works when the minor allele frequency is far from 0.5 and both files report frequencies on the same allele.

Odds ratios logged, or not. GWAS Catalog files sometimes report an odds ratio and sometimes a beta. If a column named beta has a median near 1.0, it is an odds ratio and needs log(). Feeding odds ratios to LDpred2 produces a score that correlates with nothing.

λ_GC above 1.1 in a small cohort. Add PCs, check relatedness with KING, and check for a batch effect by regressing PC1 on genotyping plate. If a single genomic region dominates the QQ plot tail, the cause is probably an unmasked long-range LD region or an inversion.

Out of memory in R. If you find yourself calling as.matrix() on a genotype object, stop. Everything in bigsnpr operates on the file-backed matrix through ind.row and ind.col index vectors. Check your disk as well, because the .bk file is not cleaned up automatically and accumulates across runs.

Sharing individual-level genotypes. You often cannot, for good reason. Summary-statistic methods cover most of what you need, and there is an established literature on secure multi-party computation for association testing when raw genotypes must stay in place ⁸. Newer work surveys where machine learning helps with the translation gap between association signals and usable interpretation ⁹. Before you design a study at all, read a plain treatment of when GWAS is the right tool, because for monogenic or family-structured questions it usually is not ¹⁰.

Woolf Software builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

1. Gaurav Thareja, Yasser Al-Sarraj, Aziz Belkadi, et al. Whole genome sequencing in the Middle Eastern Qatari population identifies genetic associations with 45 clinically relevant traits. Nature Communications, 2021. https://doi.org/10.1038/s41467-021-21381-3 ↩

2. Yuval B. Simons, Kevin Bullaughey, Richard R. Hudson, et al. A population genetic interpretation of GWAS findings for human quantitative traits. PLOS Biology, 2018. https://doi.org/10.1371/journal.pbio.2002985 ↩

3. Abdel Abdellaoui, Loic Yengo, Karin J.H. Verweij, et al. 15 years of GWAS discovery: Realizing the promise. The American Journal of Human Genetics, 2023. https://doi.org/10.1016/j.ajhg.2022.12.011 ↩

4. Andrew Anand Brown, Ana Viñuela, Olivier Delaneau, et al. Predicting causal variants affecting expression by using whole-genome sequencing and RNA-seq from multiple human tissues. Nature Genetics, 2017. https://doi.org/10.1038/ng.3979 ↩

5. Gulnara R Svishcheva, Nadezhda M Belonogova, Irina V Zorkoltseva, et al. Gene-based association tests using GWAS summary statistics. Bioinformatics, 2019. https://doi.org/10.1093/bioinformatics/btz172 ↩

6. Marquitta J. White, Brian L. Yaspan, Olivia J. Veatch, et al. Strategies for Pathway Analysis Using GWAS and WGS Data. Current Protocols in Human Genetics, 2018. https://doi.org/10.1002/cphg.79 ↩

7. Urko M. Marigorta, Juan Antonio Rodríguez, Greg Gibson, et al. Replicability and Prediction: Lessons and Challenges from GWAS. Trends in Genetics, 2018. https://doi.org/10.1016/j.tig.2018.03.005 ↩

8. Yihua Zhang, Marina Blanton, Ghada Almashaqbeh. Secure distributed genome analysis for GWAS and sequence comparison computation. BMC Medical Informatics and Decision Making, 2015. https://doi.org/10.1186/1472-6947-15-s5-s4 ↩

9. Jie Huang, Gary R. McLean, Andre Franke. Twenty years of genome-wide association studies: Health translation challenges and AI opportunities. European Journal of Human Genetics, 2025. https://doi.org/10.1038/s41431-025-01951-5 ↩

10. Fan Wang. Genome-wide association studies (GWAS): What are they, when to use them?. Rigor and Reproducibility in Genetics and Genomics, 2024. https://doi.org/10.1016/b978-0-12-817218-6.00004-8 ↩