# How to Perform Statistical Analysis in Python for GWAS Data By the end of this guide you will be able to take a genotype matrix and a phenotype file, run a full genome-wide association study (GWAS) in Python, and produce a results table with effect sizes, standard errors, p-values, corrected significance calls, principal-component covariates, and diagnostic plots. Every step includes real commands, working code, the expected output, and the checks you should run before trusting the result. This guide assumes you are comfortable with the command line and basic Python. Nothing here requires a cluster. A laptop with 16 GB of RAM handles the worked example, which uses a simulated dataset of 2,000 individuals and 100,000 single nucleotide polymorphisms (SNPs). ## What You Need Before You Start Install Python 3.10 or newer, then install these packages: ```bash conda create -n gwas python=3.11 -y conda activate gwas pip install numpy pandas scipy statsmodels scikit-learn matplotlib bed-reader ``` If you want to use PLINK 2 for quality control and large-scale regression, download the static binary from the PLINK 2 website at cog-genomics.org/plink/2.0. PLINK is a standalone command-line program, not a Python package, and it remains the fastest general-purpose tool for the heaviest GWAS steps. A common and defensible pattern in python for bioinformatics workflows is to run QC and [principal component analysis](/blog/data-analysis/principal-component-analysis-explained) (PCA) in PLINK, then do modeling, correction, and plotting in Python. Have these three files on hand before you start: 1. Genotype data in PLINK binary format (.bed, .bim, .fam) or in Variant Call Format (VCF). 2. A phenotype file, one row per sample, with an ID column matching the .fam file. 3. A covariate file with age, sex, batch, or other variables you want to adjust for. ## GWAS Concepts in Brief A GWAS tests, one variant at a time, whether a genetic variant is statistically associated with a trait. The core vocabulary you need: - **SNP**: a single nucleotide polymorphism, a position in the genome where one DNA letter varies between individuals. - **Genotype**: the pair of alleles an individual carries at that position. For a biallelic SNP with alleles A and G, genotypes are AA, AG, or GG. - **Allele dosage**: the count of one allele, coded 0, 1, or 2. Dosage 2 for the G allele means GG. This additive coding is the default in almost all GWAS regression models. - **Phenotype**: the trait you are studying. Binary phenotypes (case versus control) use [logistic regression](/blog/data-analysis/logistic-regression-definition-formula). Quantitative phenotypes (height, blood pressure, a lab measurement) use linear regression. - **Minor allele frequency (MAF)**: how often the less common allele appears in the population. Variants with very low MAF give unstable statistics and are usually filtered out. - **[Linkage disequilibrium](/knowledge/bioinformatics/linkage-disequilibrium-and-haplotype-mapping) (LD)**: nonrandom correlation between nearby variants. Two SNPs in high LD carry nearly the same information, which matters for PCA and for interpreting hits. - **Effect size**: the coefficient beta from the regression. For a quantitative trait, beta is the change in trait value per copy of the allele. For a binary trait, exponentiating the logistic coefficient gives an odds ratio. - **Covariates**: non-genetic variables included in the model, such as age, sex, batch, and ancestry principal components, so the SNP effect is not confounded by them. - **Population stratification**: systematic ancestry differences that create spurious associations. Controlling for ancestry is the single most important covariate decision in a GWAS. The conventional genome-wide significance threshold for a human GWAS with roughly one million independent common variants is 5 × 10⁻⁸, which corresponds to a Bonferroni correction for about one million tests at an alpha of 0.05. For organisms or arrays with different numbers of tests, compute the threshold yourself, as shown in Step 4. ## Setting Up Your Python Environment Verify the installation before writing any analysis code: ```bash python -c "import numpy, pandas, statsmodels, sklearn, matplotlib; print('ok')" ``` Expected output: ``` ok ``` Check the versions so your results are reproducible: ```python import numpy, pandas, statsmodels, sklearn print(numpy.__version__, pandas.__version__, statsmodels.__version__, sklearn.__version__) ``` Expected output on a current install: ``` 1.26.4 2.2.2 0.14.4 1.5.2 ``` Record these versions in a requirements.txt file committed alongside your code. A GWAS result that cannot be re-run with the same package versions cannot be fully checked later. ## Preparing and Loading GWAS Data ### Reading the PLINK map and sample files The .bim file lists variants. The .fam file lists samples. Both are whitespace-delimited text with no header row: ```python import pandas as pd bim = pd.read_csv( "data.bim", sep=r"\s+", header=None, names=["chr", "rsid", "cm", "pos", "a1", "a2"] ) fam = pd.read_csv( "data.fam", sep=r"\s+", header=None, names=["fid", "iid", "pat", "mat", "sex", "pheno"] ) print(bim.shape) print(fam.shape) print(bim.head()) ``` Expected output for the worked example: ``` (100000, 6) (2000, 6) chr rsid cm pos a1 a2 0 1 rs_sim_0001 0.0 10583 T C 1 1 rs_sim_0002 0.0 22417 A G 2 1 rs_sim_0003 0.0 31942 G A 3 1 rs_sim_0004 0.0 47120 C T 4 1 rs_sim_0005 0.0 60233 A G ``` In the .fam file, the standard PLINK convention encodes the phenotype column as 1 for control and 2 for case, with 0 or -9 meaning missing. Convert it early so you never mix up the coding: ```python fam["pheno"] = fam["pheno"].replace({1: 0, 2: 1, 0: None, -9: None}) ``` ### Reading the genotype matrix The .bed file is a compressed binary genotype matrix. The bed-reader package reads it directly into a NumPy array: ```python from bed_reader import open_bed bed = open_bed("data.bed") G = bed.read() # shape: n_samples x n_snps, float64, missing = NaN print(G.shape) print(G[:3, :5]) ``` Expected output: ``` (2000, 100000) [[0. 1. 2. 0. 1.] [1. 2. 0. 1. 0.] [0. 0. 1. 2. 2.]] ``` Three properties of this matrix matter for everything that follows: 1. Rows are samples, columns are SNPs, in exactly the order of the .fam and .bim files. Never sort one file and not the other. 2. Values are allele dosages, 0, 1, or 2, matching the a1 allele count. 3. Missing genotypes appear as NaN, and you must decide what to do with them before regression, because linear algebra functions raise errors on NaN input. If your data is in VCF, convert it to PLINK format first. This is simpler and faster than parsing VCF text in Python: ```bash plink2 --vcf cohort.vcf.gz --make-bed --out cohort ``` ### Loading phenotypes and covariates ```python pheno = pd.read_csv("phenotypes.csv") covs = pd.read_csv("covariates.csv") df = pheno.merge(covs, on="iid", how="inner").merge( fam[["iid", "pheno"]], on="iid", how="left" ) assert df["iid"].is_unique print(df.shape) ``` Expected output: ``` (2000, 6) ``` The `how="inner"` merge silently drops samples that lack covariates. Print how many samples you started with and how many survived the merge. If the number drops unexpectedly, the ID format differs between files, for example a family prefix present in one file and absent in the other. ## Step-by-Step Python GWAS Statistical Analysis Walkthrough ### Step 1. Quality control on variants and samples ![Flowchart of GWAS quality control: sample missingness, SNP missingness, MAF, HWE, then LD pruning before PCA.](/article-images/illustrations/python-gwas-statistical-analysis-qc-order.webp) GWAS quality control (QC) has a standard sequence. Applied to the .bim, .fam, and genotype data, the widely used thresholds are: | Filter | What it removes | Typical threshold | |--|--|--| | SNP missingness | Variants genotyped poorly in many samples | Missing rate above 2 percent | | Sample missingness | Samples with poorly genotyped DNA | Missing rate above 2 percent | | Minor allele frequency | Rare variants with unstable statistics | MAF below 0.01 to 0.05 | | Hardy-Weinberg equilibrium (HWE) | Genotyping artifacts, in controls for case-control data | Exact test p below 1 × 10⁻⁶ | | Heterozygosity outliers | Possible contamination or sample swaps | More than 3 to 4 standard deviations from the mean | Hardy-Weinberg equilibrium (HWE) is the expectation that genotype frequencies follow the proportions p², 2pq, q² for allele frequencies p and q in a randomly mating population. A severe deviation often signals a genotyping problem rather than real biology. Do QC in PLINK 2, which implements exact HWE tests and handles the bookkeeping: ```bash plink2 --bfile data \ --mind 0.02 \ --geno 0.02 \ --maf 0.01 \ --hwe 1e-6 \ --make-bed --out data_qc ``` Parameter notes: - `-mind 0.02` drops samples with more than 2 percent missing genotypes. Run it first, because sample removal changes SNP missing rates. - `-geno 0.02` drops SNPs missing in more than 2 percent of samples. - `-maf 0.01` keeps variants with MAF of at least 1 percent. For very large cohorts, many analysts raise this to 0.05. - `-hwe 1e-6` drops SNPs failing the HWE exact test. For case-control data, restrict this to controls by adding `midp` or running the filter on controls only, because true disease-associated loci can legitimately deviate from HWE among cases. Expected PLINK output at the end of the run: ``` Pruned 1247 variants; leaving 98753. Keys: --maf 98753 variants remaining after main filters. ``` The exact counts differ per dataset. What you should confirm is that the log reports how many variants each filter removed, and that the numbers are plausible. A filter that removes half the genome usually indicates a strand or allele-coding problem, not real QC signal. For PCA, additionally prune variants in LD so that correlated SNPs do not dominate the ancestry structure: ```bash plink2 --bfile data_qc --indep-pairwise 500 50 0.2 --out ld_prune plink2 --bfile data_qc --extract ld_prune.prune.in --pca 10 --out data_pca ``` `-indep-pairwise 500 50 0.2` walks the genome in windows of 500 SNPs, advancing 50 SNPs per step, and drops one of any pair with r² above 0.2. The `-pca 10` command writes the top 10 principal components to data_pca.eigenvec, one row per sample. ### Step 2. Load the genotype matrix and phenotypes After QC, reload everything and confirm the dimensions agree: ```python import numpy as np from bed_reader import open_bed bed = open_bed("data_qc.bed") fam_ids = pd.read_csv("data_qc.fam", sep=r"\s+", header=None, usecols=[1])[1].tolist() df = df.set_index("iid").loc[fam_ids].reset_index() G = bed.read() y = df["pheno"].to_numpy(dtype=float) print(G.shape, len(y)) ``` Expected output: ``` (2000, 98753) 2000 ``` The `.loc[fam_ids]` reindexing step is your guarantee that row 500 of the genotype matrix belongs to the same person as row 500 of the phenotype table. If any sample ID is missing, pandas raises `KeyError` there, which is exactly the failure you want rather than a silent misalignment. ### Step 3. Fit the association model #### Logistic regression for case-control phenotypes For each SNP, fit a logistic regression of case status on allele dosage plus covariates: ```python import statsmodels.api as sm cov_names = ["age", "sex", "pc1", "pc2", "pc3", "pc4", "pc5"] C = df[cov_names].to_numpy(dtype=float) def logistic_gwas_snp(g, y, C): """Return beta, se, p for one SNP. g is a dosage column, may contain NaN.""" keep = ~np.isnan(g) if keep.sum() < 10 or np.std(g[keep]) == 0: return np.nan, np.nan, np.nan X = np.column_stack([np.ones(keep.sum()), g[keep], C[keep]]) try: fit = sm.GLM(y[keep], X, family=sm.families.Binomial()).fit() return fit.params[1], fit.bse[1], fit.pvalues[1] except Exception: return np.nan, np.nan, np.nan beta, se, p = logistic_gwas_snp(G[:, 42], y, C) print(f"beta={beta:.4f} se={se:.4f} p={p:.3e}") ``` Expected output for a variant with a real but modest effect: ``` beta=0.1985 se=0.0621 p=1.368e-03 ``` Parameter notes: - The first column of X is the intercept, added explicitly with `np.ones`. Statsmodels does not add intercepts automatically when you pass a raw matrix. - `fit.params` is the SNP coefficient because the SNP sits in column 1, after the intercept. If you reorder X, reorder the index. - The odds ratio is `np.exp(beta)`. With beta = 0.1985, the odds ratio is 1.22, meaning each additional copy of the allele multiplies the odds of being a case by about 1.22. - Missing genotypes are handled by dropping that SNP's missing rows only, which is why the row count inside the function differs per SNP. Dropping rows per SNP is correct. Filling NaN with the SNP's mean dosage, called mean imputation, is also common and faster, and it keeps every regression on the full sample. Wrap the function in a loop over SNPs, storing results in a preallocated array rather than a Python list of tuples, and print progress every 10,000 variants: ```python n_snps = G.shape[1] out = np.full((n_snps, 3), np.nan) for j in range(n_snps): out[j, 0], out[j, 1], out[j, 2] = logistic_gwas_snp(G[:, j], y, C) if (j + 1) % 10000 == 0: print(f"{j + 1}/{n_snps} done") ``` #### Linear regression for quantitative phenotypes For a quantitative trait, replace the GLM with ordinary least squares (OLS): ```python def linear_gwas_snp(g, yq, C): keep = ~np.isnan(g) if keep.sum() < 10 or np.std(g[keep]) == 0: return np.nan, np.nan, np.nan X = np.column_stack([np.ones(keep.sum()), g[keep], C[keep]]) try: fit = sm.OLS(yq[keep], X).fit() return fit.params[1], fit.bse[1], fit.pvalues[1] except Exception: return np.nan, np.nan, np.nan ``` Before running OLS on a quantitative trait, inspect the phenotype distribution. A strongly skewed trait, such as a lipid concentration, is usually log-transformed or rank-normalized first, because OLS residuals should be approximately normal for the p-values to be trustworthy. #### The vectorized shortcut for linear models The per-SNP loop is clear but slow in pure Python. For linear models there is a fast algebraic trick: residualize the phenotype and genotypes on the covariates once, then every per-SNP regression reduces to a correlation. This is the standard trick behind fast GWAS software: ```python from scipy import stats # Residualize phenotype and genotype matrix on covariates (plus intercept) Xc = np.column_stack([np.ones(len(yq)), C]) Xc_inv = np.linalg.pinv(Xc) yq_resid = yq - Xc @ (Xc_inv @ yq) # Process genotypes in chunks of 10,000 SNPs n_snps = G.shape[1] betas = np.full(n_snps, np.nan) ses = np.full(n_snps, np.nan) pvals = np.full(n_snps, np.nan) chunk = 10000 for start in range(0, n_snps, chunk): g = G[:, start:start + chunk] g = np.where(np.isnan(g), np.nanmean(g, axis=0), g) # mean impute g_resid = g - Xc @ (Xc_inv @ g) n = g_resid.shape[0] g_ss = np.sum(g_resid**2, axis=0) y_ss = np.sum(yq_resid**2) cov_gy = g_resid.T @ yq_resid r = cov_gy / np.sqrt(g_ss * y_ss) t = r * np.sqrt((n - Xc.shape[1] - 1) / (1 - r**2 + 1e-12)) betas[start:start + chunk] = cov_gy / g_ss ses[start:start + chunk] = np.sqrt((1 - r**2) * y_ss / g_ss / (n - Xc.shape[1] - 1)) pvals[start:start + chunk] = 2 * stats.t.sf(np.abs(t), df=n - Xc.shape[1] - 1) ``` On the simulated dataset of 2,000 samples and about 100,000 SNPs, this completes in a few minutes on a laptop, versus tens of minutes for the per-SNP loop. The key parameter is `chunk`: 10,000 SNPs keeps each block around 16 MB in float64, which fits comfortably in cache and avoids loading the whole matrix twice. Note the mean imputation line. Imputing NaN with the SNP's mean dosage is an approximation that slightly shrinks effect sizes, and it is acceptable when missingness is below about 2 percent, which QC already guarantees. If your missingness is higher, go back and fix QC rather than relying on imputation. ### Step 4. Correct for multiple testing ![Worked calculation dividing 0.05 by 98,753 tests to get a Bonferroni threshold of 5.06e-07.](/article-images/illustrations/python-gwas-statistical-analysis-bonferroni-threshold.webp) You tested tens of thousands of variants, so some p-values below 0.05 appear by chance alone. Two corrections are standard: **Bonferroni** divides the alpha level by the number of tests. It is strict and simple: ```python M = np.sum(~np.isnan(pvals)) bonf_threshold = 0.05 / M print(f"Tests: {M}, Bonferroni threshold: {bonf_threshold:.2e}") ``` Expected output: ``` Tests: 98753, Bonferroni threshold: 5.06e-07 ``` For a human GWAS with around one million independent common variants, this same calculation gives the familiar 5 × 10⁻⁸ threshold. Compute it from your own variant count rather than copying 5 × 10⁻⁸ when your test count differs. **False discovery rate (FDR)** controls the expected proportion of false positives among your declared hits. It is more permissive and suits exploratory analyses: ```python from statsmodels.stats.multitest import multipletests valid = ~np.isnan(pvals) rej, qvals, _, _ = multipletests(pvals[valid], alpha=0.05, method="fdr_bh") results = pd.DataFrame({ "rsid": bim["rsid"].to_numpy()[valid], "beta": betas[valid], "se": ses[valid], "p": pvals[valid], "q": qvals, }) results["or"] = np.exp(results["beta"]) results.to_csv("gwas_results.txt", sep="\t", index=False, na_rep="NA") ``` Expected head of `gwas_results.txt`: ``` rsid beta se p q or rs_sim_0001 0.0121 0.0588 8.37e-01 9.21e-01 1.0122 rs_sim_0002 -0.0443 0.0601 4.61e-01 6.85e-01 0.9567 rs_sim_0003 0.2103 0.0597 4.30e-04 4.21e-03 1.2340 ``` Here `q` is the FDR-adjusted value. A variant with q below 0.05 is expected to be a false positive at most 5 percent of the time among all variants you call at that threshold. ### Step 5. Population stratification with PCA Even after including ancestry covariates in the model, check that they are doing their job. Load the PLINK PCA output and add the components to your covariate matrix: ```python pcs = pd.read_csv( "data_pca.eigenvec", sep=r"\s+", header=None, usecols=[1] + list(range(2, 12)), names=["iid"] + [f"pc{i}" for i in range(1, 11)] ) df = df.merge(pcs, on="iid", how="left") ``` Then inspect the relationship between PCs and the phenotype. If the phenotype differs strongly across PC-defined ancestry clusters and you did not include PCs, your top hits will cluster along ancestry gradients. A quick check on PC 1 versus PC 2: ```python import matplotlib.pyplot as plt fig, ax = plt.subplots(figsize=(5, 4)) ax.scatter(df["pc1"], df["pc2"], s=8, alpha=0.5) ax.set_xlabel("PC1") ax.set_ylabel("PC2") fig.savefig("pca_scatter.png", dpi=150) ``` If you see discrete clusters, your cohort mixes ancestries. Include enough PCs to absorb that structure. The common practice is 5 to 10 PCs for well-mixed human cohorts, with the number justified by how much phenotypic variance the PCs explain, rather than copied blindly. Also consider analyzing each ancestry cluster separately and meta-analyzing the results, which is the more rigorous solution when clusters are distinct. ### Step 6. Diagnostic plots and the genomic inflation factor #### Quantile-quantile plot A quantile-quantile (QQ) plot compares your observed p-values against what you would expect under the null hypothesis of no association. Systematic departure above the diagonal indicates inflation, usually from stratification or batch effects: ```python observed = np.sort(results["p"].dropna().to_numpy()) n = len(observed) expected = (np.arange(1, n + 1) - 0.5) / n fig, ax = plt.subplots(figsize=(5, 5)) ax.scatter(-np.log10(expected), -np.log10(observed), s=4) lims = [0, max(-np.log10(observed).max(), 6)] ax.plot(lims, lims, color="red", lw=1) ax.set_xlabel("Expected -log10(p)") ax.set_ylabel("Observed -log10(p)") fig.savefig("qq_plot.png", dpi=150) ``` Expected pattern: points hugging the red diagonal with a small number of points pulling upward at the tail. Those tail points are your candidate associations. #### Genomic inflation factor The genomic inflation factor lambda (λ) quantifies inflation using test statistics converted to one-degree-of-freedom chi-square values, which for a z or t statistic is simply the square. The expected median of a chi-square distribution with one degree of freedom is approximately 0.455: ```python z = results["beta"] / results["se"] lambda_gc = np.nanmedian(z**2) / 0.455 print(f"lambda_GC = {lambda_gc:.4f}") ``` Expected output for a well-controlled analysis: ``` lambda_GC = 1.0123 ``` A λ close to 1 means no systematic inflation. Values between 1.00 and about 1.05 are common even in clean analyses and can partly reflect genuine polygenic architecture. Values well above that signal stratification or batch effects, and the standard response is to add more PCs or re-examine QC rather than to accept the hits. The same statistic underlies the intercept of a linear mixed model analysis, which is the more refined approach when you have close relatives in the sample. #### Manhattan plot A Manhattan plot shows -log10(p) for every variant in genomic order: ```python results_plot = results.copy() results_plot["chr"] = bim["chr"].to_numpy()[valid] results_plot["pos"] = bim["pos"].to_numpy()[valid] fig, ax = plt.subplots(figsize=(12, 4)) chroms = results_plot["chr"].unique() colors = ["#4c72b0", "#dd8452"] x = 0 xticks, xlabels = [], [] for i, chrom in enumerate(chroms): sub = results_plot[results_plot["chr"] == chrom] ax.scatter(sub["pos"] + x, -np.log10(sub["p"]), s=3, color=colors[i % 2]) xticks.append(x + sub["pos"].mean()) xlabels.append(str(chrom)) x += sub["pos"].max() ax.axhline(-np.log10(bonf_threshold), color="red", lw=1, ls="--") ax.set_xticks(xticks) ax.set_xticklabels(xlabels) ax.set_xlabel("Chromosome") ax.set_ylabel("-log10(p)") fig.savefig("manhattan.png", dpi=150, bbox_inches="tight") ``` Points above the dashed red line meet the Bonferroni threshold. Alternating colors by chromosome make it easy to see whether hits cluster spuriously at chromosome edges, which often indicates a mapping or alignment problem. ## Key Options at a Glance | Decision | Options | When to use each | |--|--|--| | Regression model | Logistic (GLM Binomial), linear (OLS) | Logistic for case-control, linear for quantitative traits | | Genotype coding | Additive 0/1/2, dominant, recessive | Additive is the default and most powerful for common variants | | Missing genotypes | Drop rows per SNP, mean imputation | Drop when sample is large, impute for speed | | Multiple testing | Bonferroni, FDR (Benjamini-Hochberg) | Bonferroni for strict single-variant claims, FDR for discovery ranking | | Significance threshold | 0.05 divided by number of tests | Compute from your own test count, 5 × 10⁻⁸ for ~1M human tests | | Stratification control | Covariate PCs, linear mixed models | PCs for unrelated samples, mixed models when relatives are present | | QC engine | PLINK 2, Python loops | PLINK for full-genome QC and PCA, Python for modeling and plotting | ## Tips for Scaling Up A few practical tips that separate a working analysis from a painful one: 1. **Keep the heavy lifting outside Python when genomes are large.** For whole-genome human data with millions of variants, run the regression itself in PLINK 2 with `-glm`, then read the resulting summary statistics into Python for correction, annotation, and plotting. The Python statistics you compute above and the PLINK output should agree to within rounding. 2. **Use memory-mapped genotypes.** `bed_reader.open_bed()` supports reading a range of SNPs without loading the file, for example `bed.read(index=np.s_[:, start:stop])`, which keeps memory flat regardless of cohort size. 3. **Parallelize across chromosomes.** Split your variant list by chromosome, run independent Python processes, and concatenate the result tables. The per-SNP regressions share no state, so this parallelizes trivially. 4. **Annotate results with real gene positions.** This is where Biopython earns its place in a GWAS pipeline. Use `Bio.Entrez` to query NCBI databases or parse a reference annotation with `Bio.SeqIO`, then map lead SNPs to nearby genes by chromosomal position. A full Biopython tutorial is beyond this article's scope, but the mapping step needs only chromosome, position, and a gene interval table. 5. **Version everything.** Save the exact PLINK command lines, the requirements.txt file, the QC log files, and a random seed for any simulation. A GWAS you cannot re-run is a GWAS you cannot defend. ## Common Mistakes and Limitations These are the errors that most often invalidate a GWAS result, with the specific fix for each. **Skipping or misordering QC.** Running `-geno` before `-mind` double-counts the same bad data, because sample removal changes variant missingness. Follow the order: samples first, then variants, then MAF, then HWE. **Ignoring population stratification.** An ancestry-correlated trait analyzed without PC covariates produces a QQ plot lifted off the diagonal across its whole length, not just at the tail. If your λ is above roughly 1.05 and the QQ plot lifts uniformly, stop interpreting individual hits and fix the model. **Pruning before PCA is skipped.** Running PCA on all variants lets long regions of high LD, and even a single very large inversion, dominate PC1. Always LD-prune before PCA. A PC that separates samples along one chromosome rather than by ancestry is a red flag for exactly this problem. **Confusing statistical significance with effect size.** In a GWAS of 500,000 people, a variant can reach 5 × 10⁻⁹ with an odds ratio of 1.01, a biologically negligible effect. Always report beta or the odds ratio alongside the p-value, and never sort your candidate list by p-value alone. **Testing every phenotype and every model variant.** Each additional phenotype, covariate set, or transformation multiplies your multiple-testing burden. Decide your primary phenotype and primary model before looking at results, and treat everything else as exploratory with FDR control. **Misreading the p-value as the probability the association is real.** A p-value is the probability of seeing data at least this extreme if the null hypothesis were true. It is not the probability that the null is false. **Failing to check sample ID alignment.** The single most dangerous silent error in this workflow is a genotype row matched to the wrong phenotype row, which produces confident garbage. The explicit `.loc[fam_ids]` reindexing in Step 2 exists to make that failure loud instead of silent. **Interpreting unreplicated findings.** A single GWAS hit, however significant, is a hypothesis until it replicates in an independent cohort or has functional support. Plan replication before you start, not after you see the results. **Rare variants in a standard GWAS.** Variants with MAF below 1 percent lack power in single-variant tests and are filtered by the MAF threshold. They require dedicated burden orSKAT-style tests, which are a different analysis, not a looser threshold on the same one. One more limitation belongs in this list: this workflow tests one variant at a time and controls structure with covariates. Linear mixed model methods, which model the genetic relatedness matrix directly, are the standard for cohorts containing relatives or fine-scale structure, and implementing them is beyond what plain statsmodels offers. ## Frequently Asked Questions ### Do I need PLINK if I am doing GWAS analysis in Python? No, but it helps. Pure Python with pandas, NumPy, and statsmodels handles small to moderate studies completely, while PLINK 2 is faster for full-genome QC, PCA, and regression on very large cohorts. A hybrid workflow, PLINK for QC and Python for modeling and visualization, is the most common production pattern. ### What p-value threshold should I use for genome-wide significance? Divide 0.05 by your number of independent tests. With about one million independent common variants in a human GWAS this gives the conventional 5 × 10⁻⁸. For fewer tests, such as an exome array or a non-human genome, your threshold will be less strict, so compute it rather than copying the human convention. ### Logistic or linear regression: which do I use for my phenotype? Use logistic regression for binary case-control phenotypes and linear regression for quantitative traits measured on a continuous scale. If your trait is an ordered category with few levels, consider whether dichotomizing it or using a model for ordinal outcomes fits better than forcing it into either standard model. ### Can Python handle millions of SNPs in a single run? Pure Python loops cannot handle millions of SNPs efficiently, but memory-mapped readers and the vectorized residualization approach scale much further. For genuinely genome-scale human data, run the association scan in PLINK 2 and use Python for the downstream statistics, correction, and plotting. ### How many principal components should I include as covariates? Five to ten PCs is the common starting range for human cohorts, but the right number depends on your data. Increase the count until the genomic inflation factor stabilizes near 1 and the QQ plot no longer lifts uniformly across its length. Distinct ancestry clusters may warrant stratified analyses instead of more PCs. ### How should I handle missing genotype data in the regression? The two standard options are dropping the rows with missing values for each SNP, or imputing them with the SNP's mean dosage. Dropping is statistically cleanest, mean imputation is faster and acceptable once QC has already limited missingness to a few percent. Either way, verify that missingness never exceeds your QC threshold. ### What does the genomic inflation factor tell me? Lambda is the ratio of your median chi-square test statistic to its null expectation of about 0.455. A value near 1 means your test statistics behave as expected under the null. Values well above 1 indicate systematic inflation from stratification, batch effects, or cryptic relatedness, and the fix is better covariates or a mixed model, not a looser significance threshold. ### Where does Biopython fit into a GWAS workflow? Biopython handles the sequence and database side of python for bioinformatics, not the regression itself. Typical GWAS uses are retrieving reference sequences with Bio.Entrez, parsing genome annotation to map significant variants to genes, and validating allele assignments against a reference. For high-performance VCF parsing, dedicated libraries such as pysam or cyvcf2 are generally faster than Biopython's parsers.