A reference for navigating the data pipeline — from the GWAS Catalog through allele frequency computation, FST, Watterson's θ, and nucleotide diversity (π). The code is yours to write; this guide explains what each step needs to accomplish and what to watch out for along the way.
Your central question: do cancer-associated GWAS variants show different allele frequency patterns across global populations compared to random background variants? Three data sources connect through a shared key — the rsID — and you will build that connection step by step.
| Code | Label | Example sub-populations |
|---|---|---|
| AFR | African | YRI, LWK, GWD, MSL, ESN, ACB, ASW |
| AMR | Admixed American | MXL, PUR, CLM, PEL |
| EAS | East Asian | CHB, JPT, CHS, CDX, KHV |
| EUR | European | CEU, TSI, FIN, GBR, IBS |
| SAS | South Asian | GIH, PJL, BEB, STU, ITU |
Before writing any code, make sure you understand the structure of all three input sources. Misreading a column or index is the most common source of errors in this project.
Downloaded from ebi.ac.uk/gwas. Tab-separated with a header row. Key columns:
Three problems you will encounter and must fix.
(1) Allele-annotated rsIDs in the SNPS column — entries like rs16901979-A are rsID + risk allele joined by a dash. Strip the allele suffix before using the rsID for coordinate lookup.
(2) Multiple rsIDs in one row — entries like rs10187424;rs1859962 need to be split so each row represents exactly one SNP. After splitting, verify the number of rsIDs equals the number of positions — rows where they disagree are ambiguous and should be dropped.
(3) SNP_ID_CURRENT lacks the "rs" prefix — the catalog stores updated rsIDs as bare numbers (e.g., 4430796). Prepend "rs" before using this column, and fall back to SNPS when it is empty.
CHR_ID needs cleaning. The chromosome column can contain non-numeric entries: "X", "Y", "MT", combined values like "1;2", or entries with extra characters. The frq.count files only cover autosomes 1–22 — filter to integer chromosome IDs only, then convert to integer type for sorting.
No header row. Tab-separated. Large (~7 GB compressed.) The columns are at fixed zero-based positions. When you write your output file you will need to assign your own column names — we suggest the names below so they are easier to work with downstream. Only five columns matter for this project:
Use posEnd (col[3]), not posStart (col[2]). The frq.count files record SNP positions as 1-based coordinates that match posEnd. If you join on posStart you will lose the vast majority of your matches — posStart is 0-based and off by one. The chromosome column uses a "chr" prefix (e.g., chr1) — strip this when joining to frq.count files, which use bare integers.
One file per super-population × chromosome (e.g., AFR_chr6.frq.count). Tab-separated with a header row. The last two columns encode allele-count pairs.
| Column | Meaning | How you use it |
|---|---|---|
CHROM | Chromosome number (integer, no "chr" prefix) | Filter to your focal chromosome |
POS | Position (1-based) — matches posEnd from snp151 (col[3]) | Join key to your coordinate table |
N_CHR | Total chromosomes sampled in this population | Denominator for all frequency calculations |
{ALLELE:COUNT} | Space-separated allele:count pairs (e.g., T:1285 C:1) | Split on ":" to get allele identity and count |
Allele frequency formula. For each allele: frequency = count ÷ N_CHR. The minor allele frequency (MAF) is the smaller of the two frequencies — always between 0 and 0.5. Use MAF for all your distribution plots so comparisons across populations are not affected by which allele is labeled "reference."
Goal. Load all your GWAS Catalog TSV files, clean the messy fields, identify which chromosome has the most associated variants (your focal chromosome), and output two clean files: a full processed table and a plain list of rsIDs for the coordinate lookup.
Place all your downloaded TSVs in one folder. Read every file matching the naming pattern, stack them into a single data frame, and parse the cancer type from the filename. Report the total number of rows, unique studies (PubMed IDs), and unique disease traits — these numbers go in your Methods section.
Address the three problems described in the file format section above. After cleaning:
^rs[0-9]+$ — a bare rsID with no allele suffix.SNP_ID_CURRENT as your primary rsID (it reflects rsID merges), falling back to SNPS when SNP_ID_CURRENT is empty. Make sure to prepend "rs" since the catalog omits it.Check your counts at each filter step. Print how many rows remain after each operation. Unexpectedly large drops are a sign that a regex is too aggressive. Keep a note of your starting and ending SNP counts for the Methods section.
Count the number of unique SNPs per chromosome (de-duplicate rsIDs first so you don't inflate counts). The chromosome with the most associated variants is your focal chromosome — you will do all downstream allele frequency and FST analysis on this chromosome only. Make a bar chart of SNP counts per chromosome for your Methods or Results.
You need three outputs from this step:
Ancestry in INITIAL.SAMPLE.SIZE. The INITIAL.SAMPLE.SIZE column contains free-text descriptions of who was enrolled (e.g., "1,034 European ancestry individuals"). Extracting ancestry terms from this column is a good starting point for quantifying ancestry bias in your GWAS catalog download — useful for the Limitations section.
Goal. For every rsID from Step 1, retrieve its chromosome and 1-based position (posEnd, col[3]) from snp151.txt.gz. Output a table with your own header row — columns: rsID, chr, posStart, posEnd, alleles. This table bridges the GWAS rsIDs to the positions used in the allele count files.
snp151.txt.gz is ~7 GB compressed. Loading it entirely into memory will crash most laptops. Instead, open the file as a stream, read one line at a time, check whether the rsID in column 4 (zero-indexed) is in your target set, and write matching rows to an output file. Remove each found rsID from your search set as you go — this speeds up future lookups and lets you stop early once all rsIDs are found.
Write the output with a header line: rsID chr start end alleles.
Some rsIDs will not be found. A small fraction of GWAS Catalog rsIDs have been retired, merged into other IDs, or are not present in snp151. Losing roughly 0.1–1% is normal. Write a separate file listing the missing rsIDs so you can inspect them. Report how many were lost in your Methods section.
Load the output from 2a into R. Filter to rows where the chromosome matches your focal chromosome (remember snp151 uses the chr6 format — strip the prefix when needed). Remove any duplicate rsID entries, keeping one row per rsID. The posEnd column is the 1-based position you will join against POS in the frq.count files.
Goal. For each of the five super-populations, parse the frq.count file for your focal chromosome, join to your coordinate table on POS == posEnd, compute per-allele frequencies and MAF, and identify which frequency corresponds to the GWAS risk allele. Then sample 1,000 random background SNPs per population for comparison.
Write a single function that reads one frq.count file for one population. Inside it, skip the header, split the ALLELE:COUNT strings on ":" to get allele identity and raw count, then divide each count by N_CHR to get frequency. Compute MAF as the minimum of the two allele frequencies. Call this function in a loop over all five populations and stack the results.
Join on POS == posEnd. The frq.count POS column is 1-based and matches snp151's posEnd (col[3]). Join on this; joining on posStart (col[2]) will lose most matches because posStart is 0-based and off by one.
After joining, bring in the risk allele column from your cleaned GWAS table. For each SNP in each population, look up whether the risk allele matches the REF or ALT allele in the frq.count file and record the appropriate frequency. Rows where the risk allele matches neither (e.g., allele listed as "?" in the catalog) should be filtered out. This risk allele frequency — not MAF — is what you plot for associated variants in the allele frequency figures.
Risk allele frequency ≠ MAF. The risk allele may be common in one population and rare in another. MAF is always between 0 and 0.5; risk allele frequency can be anywhere from 0 to 1. Keep both — use risk allele frequency when you want to track the specific GWAS hit, and MAF when comparing distributions.
For each super-population, draw 1,000 SNPs at random from the full frq.count file, excluding any positions that are in your associated SNP set. Sample each population independently. Set a random seed at the top of your script and use the same seed throughout — report this number in your Methods so results are reproducible.
Sample from the frq.count file, not from a list of rsIDs. Most positions in the count file do not have rsIDs in the GWAS catalog — that is the point. You are drawing from the genomic background, not from another GWAS.
Goal. Quantify population differentiation using pairwise FST across all ten super-population pairs. Then, using your site frequency spectrum data, compute two summaries of within-population diversity — Watterson's θ and nucleotide diversity (π) — for both associated and random SNP sets. Compare these statistics across populations and between variant types.
FST at a single biallelic site summarizes how much more different two populations are from each other than each is internally. Use the following estimator, where p1 and p2 are the ALT allele frequencies and n1, n2 are the chromosome counts:
There are 10 pairwise combinations of 5 populations — use combn() in R to generate them automatically. Compute per-SNP FST for every pair, then summarize with distributions (boxplots or histograms) and mean FST heatmaps.
Negative FST values. The estimator can return small negative values due to sampling noise when within-population variance exceeds between-population variance. True FST is bounded between 0 and 1. Floor negative values at 0, report in your Methods that you did this, and note how much the mean FST changes before versus after flooring.
Watterson's θ estimates the population-scaled mutation rate from the number of segregating sites. For a set of S segregating sites observed in a sample of n chromosomes:
θW summarizes diversity in terms of how many variable sites you observe, normalized for sample size. A higher θW means more polymorphic sites. Compute this separately for your associated SNP set and your random SNP set within each super-population. Because your SNP sets are pre-selected (not a random genomic window), interpret θW as a relative measure comparing variant classes, not as an absolute population parameter.
S is just your count of SNPs in that set. Every site in your frq.count-matched set is by definition segregating (it has two alleles present). So S = number of sites in your set for that population. The tricky part is computing the harmonic number an — you can do this in R with sum(1/seq(1, n-1)) where n is N_CHR for that population.
Nucleotide diversity (π) is the average number of pairwise nucleotide differences per site between any two randomly chosen sequences. For a single biallelic site, using the alternate allele count (AC) and the total number of alleles (AN) directly from your frq.count data:
which simplifies to:
This counts all discordant pairs across the AN chromosomes and normalizes by the total number of possible pairs AN(AN − 1). Average πsite across all sites in your SNP set to get mean nucleotide diversity (π̄) for that set and population. This is directly computable from the AC and AN columns you already have from Step 3.
π vs 2pq. Using allele frequencies p = AC/AN and q = 1 − p, the formula above is equivalent to 2pq × AN/(AN−1). For large sample sizes (AN in the hundreds, as in 1000 Genomes) the correction factor AN/(AN−1) is very close to 1, so π ≈ 2pq. The true pairwise π is slightly preferred because it scales the count of discordant pairs by the exact number of possible pairs, rather than the approximation from random-mating expected heterozygosity. Both are sensitive to intermediate-frequency variants and both contrast naturally with θW, which is sensitive to rare variants.
Compute both θW and π for associated and random SNPs in all five super-populations. Key comparisons to make: (1) Are associated variants more or less diverse than the random background within each population? (2) Do θW and π tell the same story, or do they diverge — and what would that divergence mean given how each statistic is calculated? (3) Which super-population shows the highest nucleotide diversity for cancer-associated loci?
θW and π are related but not identical. θW depends only on how many sites are variable, while π depends on how common each variant is. A SNP set dominated by rare variants will have a higher θW relative to π; a set enriched for intermediate-frequency variants will have the reverse. The contrast between them is itself informative about the shape of the allele frequency spectrum — and connects directly to what you saw in Figures 1 and 2.
Below are schematic mockups showing what each figure should look like and what it should communicate. The exact aesthetic, color choices, and axis formatting are up to you. Every figure must have a descriptive title, labeled axes, and the name of the person who produced it in the subtitle or caption. Use ggsave() to export PDFs at 300 dpi.
Always normalize to proportions, not raw counts, before plotting MAF distributions. Divide bin counts by the group total so each bar represents the proportion of SNPs in that bin. This makes populations directly comparable regardless of how many SNPs each group has.
Show the distribution of minor allele frequencies for your cancer-associated SNPs, faceted by super-population. Bin MAF into 0.05-wide intervals from 0 to 0.5. Normalize counts to proportions within each population. Comment on whether the distributions are similar across populations or differ — particularly whether rare variants are more enriched in some populations than others.
Figure 1 — Binned MAF distributions for cancer-associated SNPs across all five super-populations. Normalize counts to proportions within each population; each bar spans a 0.05-wide MAF interval.
Compare allele frequency distributions of cancer-associated SNPs and randomly sampled background SNPs side-by-side, using a facet grid with one panel per super-population and dodged bars within each panel. Normalize both sets to proportions before plotting.
Random SNPs are strongly skewed toward the lowest MAF bin — this is the genomic null expectation under neutral drift. If your associated SNPs are less skewed or more evenly spread, this suggests GWAS preferentially identifies common variants. If the shapes are similar, the trait may be influenced by rare variants too.
Figure 2 — MAF distributions of cancer-associated (red) vs. randomly sampled (blue) SNPs, faceted by super-population. Compare the degree of skew and comment on what differences imply about variant discovery.
Compute Watterson's θW and mean nucleotide diversity (π) for associated and random SNPs separately within each super-population. Display both statistics as grouped bar charts in two side-by-side panels. Comment on: which population shows the highest diversity for cancer-associated loci? Do θW and π agree, or do they diverge — and what does divergence tell you about the allele frequency spectrum of your SNP sets?
Figure 3 — Panel A: Watterson's θW for associated and random SNPs in each super-population. Panel B: Mean nucleotide diversity (π = 2·AC·(AN−AC) / [AN·(AN−1)]) for the same sets. The contrasting patterns between the two panels reflect differences in the allele frequency spectrum — θW is sensitive to rare variants while π is maximized at intermediate frequency.
Compute pairwise FST for all 10 population pairs across your associated SNPs. Show the distribution of per-SNP FST values as boxplots (one box per pair) and summarize mean FST as a symmetric heatmap. Comment on which pairs are most and least differentiated and whether any pattern reflects known demographic history.
Figure 4 — Left: horizontal boxplot of per-SNP FST distributions for all ten pairwise population comparisons (population pair labels on y-axis, FST on x-axis). Right: symmetric mean FST heatmap. Values are illustrative — yours will differ.
Overlay or dodge the FST distributions of associated and random SNPs for each population pair. This lets you ask whether cancer-associated variants show more or less differentiation than the genomic background. Use a facet grid with one panel per population pair (all 10).
Figure 5 — Overlapping histograms of FST per population pair, comparing associated (red) and random (blue) variants. All 10 pairwise facets should appear in your final figure.