A reference for navigating the data pipeline — from the GWAS Catalog through allele frequency computation, FST, Watterson's θ, and expected heterozygosity. 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 expected heterozygosity (He) — 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.
Expected heterozygosity at a single biallelic site is the probability that two randomly drawn chromosomes carry different alleles. With allele frequencies p and q = (1−p):
Average He,site across all sites in your SNP set to get mean expected heterozygosity for that set and population. This is directly computable from the frequency columns you already have from Step 3. Unlike θW, He is sensitive to allele frequency — sites at intermediate frequency contribute most.
Compute both θW and He 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 He 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 diversity for cancer-associated loci?
θW and He are related but not identical. Both estimate diversity, but θW depends only on how many sites are variable, while He depends on how common each variant is. A SNP set dominated by rare variants will have a higher θW relative to He; a set enriched for intermediate-frequency variants will have the reverse. The contrast between them is itself informative about the shape of the SFS — and connects directly to what you saw in Figure 5.
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 (top). Figure 2 — Example panel showing associated vs. random SNPs side-by-side in one population (bottom); your plot should have all five populations as facets.
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 3 — Boxplot of per-SNP FST distributions for all ten pairwise comparisons (left) and mean FST heatmap (right) for cancer-associated variants. Numbers in the mockup 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 4 — Overlapping histograms of FST per population pair, comparing associated (red) and random (blue) variants. All 10 pairwise facets should appear in your final figure.
Plot the distribution of MAF values for your associated SNPs, one histogram panel per super-population. Normalize to proportions within each population. The SFS is a classical population genetics summary: an excess of rare variants suggests population growth or purifying selection; a flat or intermediate-frequency-enriched SFS suggests common variation and GWAS ascertainment. Comment on the shape within each population and how populations compare.
Figure 5 — Site frequency spectrum of cancer-associated variants across all five super-populations. Bar heights represent the proportion of associated SNPs in each 0.05-wide MAF interval.
Compute Watterson's θW and mean expected heterozygosity (He) for associated and random SNPs separately within each super-population. Display both statistics as grouped bar charts, either as two side-by-side panels or as a single faceted figure. Comment on: which population shows the highest diversity for cancer-associated loci? Do θW and He agree, or do they diverge — and what does divergence tell you about the allele frequency spectrum of your SNP sets?
Figure 6 — Panel A: Watterson's θW for associated and random SNPs in each super-population. Panel B: Mean expected heterozygosity (He = average 2pq) 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 He is maximized at intermediate frequency.
Figure assignment summary. Each group member is responsible for at least one figure. Write your name in the subtitle or caption of every figure you produce. The grader uses this to assign individual scores for the Results section.
The final paper must be ≤ 12 double-spaced pages (excluding references). Sections are graded individually — write your name next to every section you are responsible for.
Write the abstract last. It summarizes everything — you cannot write it accurately until the paper is done.
Structure as a funnel: broad background → specific gap → your study.
A reader should be able to reproduce your analysis exactly from your Methods section. Address:
One paragraph per figure. For every result:
Upload all code and processed data files to a shared Google Drive folder and include the link. Your folder should contain every input and output file plus your full analysis scripts so that another person could reproduce your results from scratch.
QBIO 475 · Statistical and Evolutionary Genetics · University of Southern California · mooney-lab.github.io · Paper due December 15 at 1:00 pm