You only need to write lines of code for each question. When answering questions that ask you to identify or interpret something, the length of your response doesn’t matter. For example, if the answer is just ‘yes,’ ‘no,’ or a number, you can just give that answer without adding anything else.
+
We will go through comparable code and concepts in the live learning session. If you run into trouble, start by using the help help() function in R, to get information about the datasets and function in question. The internet is also a great resource when coding (though note that no outside searches are required by the assignment!). If you do incorporate code from the internet, please cite the source within your code (providing a URL is sufficient).
+
Please bring questions that you cannot work out on your own to office hours, work periods or share with your peers on Slack. We will work with you through the issue.
+
You will need to install PLINK and run the analyses. Please follow the OS-specific setup guide in SETUP.md. PLINK is a free, open-source whole genome association analysis toolset, designed to perform a range of basic, large-scale analyses in a computationally efficient manner.
+
+
# Set up a consistent path for all chunks of code
+knitr::opts_knit$set(root.dir =normalizePath("../../"))
+
+library(data.table)
+library(ggplot2)
+library(seqminer)
+library(HardyWeinberg)
+
+
Loading required package: mice
+
+
+
Warning: package 'mice' was built under R version 4.5.2
+
+
+
+Attaching package: 'mice'
+
+
+
The following object is masked from 'package:stats':
+
+ filter
+
+
+
The following objects are masked from 'package:base':
+
+ cbind, rbind
+
+
+
Loading required package: nnet
+
+
+
Loading required package: Rsolnp
+
+
+
Loading required package: shape
+
+
library(dplyr)
+
+
+Attaching package: 'dplyr'
+
+
+
The following object is masked from 'package:HardyWeinberg':
+
+ recode
+
+
+
The following objects are masked from 'package:data.table':
+
+ between, first, last
+
+
+
The following objects are masked from 'package:stats':
+
+ filter, lag
+
+
+
The following objects are masked from 'package:base':
+
+ intersect, setdiff, setequal, union
+
+
+
+
Question 1: Data inspection
+
Before fitting any models, it is essential to understand the data. Use R or bash code to answer the following questions about the gwa.qc.A1.fam, gwa.qc.A1.bim, and gwa.qc.A1.bed files, available at the following Google Drive link: https://drive.google.com/drive/folders/11meVqGCY5yAyI1fh-fAlMEXQt0VmRGuz?usp=drive_link. Please download all three files from this link and place them in 02_activities/data/.
+
+
Read the .fam file. How many samples does the dataset contain?
Read the .bim file. How many SNPs does the dataset contain?
+
+
+
wc-l ./02_activities/data/gwa.qc.A1.bim
+
+#There are 101083 SNPs in the dataset.
+
+
101083 ./02_activities/data/gwa.qc.A1.bim
+
+
+
+
+
Question 2: Allele Frequency Estimation
+
+
Load the genotype matrix for SNPs rs1861, rs3813199, rs3128342, and rs11804831 using additive coding. What are the allele frequencies (AFs) for these four SNPs?
+
+
+
#Create SNP list file with 4 SNPs: rs1861, rs3813199, rs3128342, and rs11804831
+setopt interactivecomments
+printf"%s\n" rs1861 rs3813199 rs3128342 rs11804831 > ./02_activities/data/A1snplist.txt
+
+#Subset SNPs and output allele frequencies
+plink2--bfile ./02_activities/data/gwa.qc.A1 \
+--extract ./02_activities/data/A1snplist.txt \
+--recode A\
+--freq\
+--out ./02_activities/data/gwa_qc_A1_freq
+
+
bash: line 1: setopt: command not found
+PLINK v2.00a5.12 M1 (25 Jun 2024) www.cog-genomics.org/plink/2.0/
+(C) 2005-2024 Shaun Purcell, Christopher Chang GNU General Public License v3
+Logging to ./02_activities/data/gwa_qc_A1_freq.log.
+Options in effect:
+ --bfile ./02_activities/data/gwa.qc.A1
+ --export A
+ --extract ./02_activities/data/A1snplist.txt
+ --freq
+ --out ./02_activities/data/gwa_qc_A1_freq
+
+Start time: Wed Mar 18 12:26:12 2026
+16384 MiB RAM detected; reserving 8192 MiB for main workspace.
+Using up to 10 threads (change this with --threads).
+4000 samples (2000 females, 2000 males; 4000 founders) loaded from
+./02_activities/data/gwa.qc.A1.fam.
+101083 variants loaded from ./02_activities/data/gwa.qc.A1.bim.
+1 quantitative phenotype loaded (4000 values).
+--extract: 4 variants remaining.
+Calculating allele frequencies... 0%0%done.
+--freq: 0%0%25%50%75%
+--freq: Allele frequencies (founders only) written to
+./02_activities/data/gwa_qc_A1_freq.afreq .
+4 variants remaining after main filters.
+
+--export A pass 1/1: loading... 0%0%writing... 0%1%2%3%4%5%6%7%8%9%10%11%12%13%14%15%16%17%18%19%20%21%22%23%24%25%26%27%28%29%30%31%32%33%34%35%36%37%38%39%40%41%42%43%44%45%46%47%48%49%50%51%52%53%54%55%56%57%58%59%60%61%62%63%64%65%66%67%68%69%70%71%72%73%74%75%76%77%78%79%80%81%82%83%84%85%86%87%88%89%90%91%92%93%94%95%96%97%98%99% done.
+--export A: ./02_activities/data/gwa_qc_A1_freq.raw written.
+End time: Wed Mar 18 12:26:12 2026
#CHROM ID REF ALT PROVISIONAL_REF? ALT_FREQS OBS_CT
+ <int> <char> <char> <char> <char> <num> <int>
+1: 1 rs3813199 G A Y 0.0569126 7942
+2: 1 rs11804831 T C Y 0.1543410 7924
+3: 1 rs3128342 C A Y 0.3051210 7928
+4: 1 rs1861 C A Y 0.0539859 7928
+
+
#The allele frequencies for the alternative allele of rs1861, rs3813199, rs3128342, and rs11804831 are as follows: 0.0539859, 0.0569126, 0.3051210, and 0.1543410.
+
+A1freq$REF_FREQ <-1- A1freq$ALT_FREQS
+A1freq
+
+
#CHROM ID REF ALT PROVISIONAL_REF? ALT_FREQS OBS_CT REF_FREQ
+ <int> <char> <char> <char> <char> <num> <int> <num>
+1: 1 rs3813199 G A Y 0.0569126 7942 0.9430874
+2: 1 rs11804831 T C Y 0.1543410 7924 0.8456590
+3: 1 rs3128342 C A Y 0.3051210 7928 0.6948790
+4: 1 rs1861 C A Y 0.0539859 7928 0.9460141
+
+
#The allele frequencies for the reference allele of rs1861, rs3813199, rs3128342, and rs11804831 are as follows: 0.9460141, 0.9430874, 0.6948790, and 0.8456590.
+
+
+
What are the minor allele frequencies (MAFs) for these four SNPs?
#The minor allele frequencies for rs1861, rs3813199, rs3128342, and rs11804831 are as follows: 0.0539859, 0.0569126, 0.3051210, and 0.1543410.
+
+
+
+
Question 3: Hardy–Weinberg Equilibrium (HWE) Test
+
+
Conduct the Hardy–Weinberg Equilibrium (HWE) test for all SNPs in the .bim file. Then, load the file containing the HWE p-value results and display the first few rows of the resulting data frame.
geno_rs1861 <- geno_A1_subset %>%
+select(FID, IID, PHENOTYPE, rs1861_C) %>%
+na.omit()
+
+#Conduct a linear regression model between rs1861 and the phenotype
+lm_rs1861 <-lm(PHENOTYPE ~ rs1861_C, data = geno_rs1861)
+summary(lm_rs1861)
+
+
+Call:
+lm(formula = PHENOTYPE ~ rs1861_C, data = geno_rs1861)
+
+Residuals:
+ Min 1Q Median 3Q Max
+-3.5439 -0.6850 0.0021 0.6993 3.3268
+
+Coefficients:
+ Estimate Std. Error t value Pr(>|t|)
+(Intercept) -0.94762 0.09486 -9.989 <2e-16 ***
+rs1861_C 0.97382 0.04943 19.703 <2e-16 ***
+---
+Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+
+Residual standard error: 1.003 on 3962 degrees of freedom
+Multiple R-squared: 0.08924, Adjusted R-squared: 0.08901
+F-statistic: 388.2 on 1 and 3962 DF, p-value: < 2.2e-16
+
+
#The p-value is less than 2 × 10⁻¹⁶. There is a statistically significant relationship between rs1861 and the phenotype.
+
+
+
How would you interpret the beta coefficient from this regression?
+
+
We performed a linear regression to test the association between genotype rs1861 and the continuous phenotype. The beta-coefficient is 0.97382. This means that for every additional C allele an individual has for rs1861, their phenotype increases by 0.97382 units.
+
+
Plot the scatterplot of phenotype versus the genotype of SNP rs1861. Add the regression line to the plot.
+
+
+
# Create a new data frame containing a sequence of genotype values
+ geno_range <-data.frame(
+rs1861_C =seq(0, 2, length.out =100)
+ )
+# Use the fitted linear regression model to predict phenotype probability
+# for each genotype value in geno_range
+ geno_range$predicted_prob <-predict(lm_rs1861, newdata = geno_range, type ="response")
+
+# Start from the observed genotype dataset
+count_data <- geno_A1_subset %>%
+# Remove rows with missing genotype or missing phenotype
+filter(!is.na(rs1861_C), !is.na(PHENOTYPE)) %>%
+# Group data by genotype value and phenotype category
+group_by(rs1861_C, PHENOTYPE) %>%
+# Count how many observations fall into each genotype-phenotype combination
+summarise(count =n(), .groups ="drop")
+
+# Start building the plot
+p1=ggplot() +
+# Add points showing observed genotype-phenotype combinations
+geom_point(
+# Use the summarized count_data for the points
+data = count_data,
+# Map genotype to x-axis, phenotype to y-axis, and count to point size
+aes(x = rs1861_C, y = PHENOTYPE, size = count),
+color ="blue"
+ ) +
+# Add a regression line showing predicted probability across genotype values
+geom_line(# Use the prediction data frame for the line
+data = geno_range,
+aes(x = rs1861_C, y = predicted_prob),
+color ="red",
+# Set line thickness
+size =1.2
+ ) +
+# Control the size range of the points
+scale_size_continuous(range =c(2, 10)) +
+labs(
+title ="Linear Regression: PHENOTYPE ~ Genotype",
+x ="Genotype (0/1/2)",
+y ="Phenotype ",
+size ="Count"
+ ) +
+# Restrict the visible x- and y-axis range
+coord_cartesian(xlim =c(-0.5, 2.5)) +
+# Apply a minimal theme for a clean appearance
+theme_minimal()
+
+
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
+ℹ Please use `linewidth` instead.
+
+
print(p1)
+
+
+
+
+
+
+
+
+
+
Convert the genotype coding for rs1861 to recessive coding.
+
+
+
geno_A1_subset_rec <- geno_A1_subset # make a copy
+
+#change into recessive coding
+geno_A1_subset_rec[, 7:ncol(geno_A1_subset_rec)] <-as.data.frame(
+lapply(geno_A1_subset[, 7:ncol(geno_A1_subset)], function(x) ifelse(x ==2, 1, 0))
+)
+
+#Subset to just rs1861
+geno_rs1861_rec <- geno_A1_subset_rec %>%
+select(FID, IID, PHENOTYPE, rs1861_C) %>%
+na.omit()
+
+
+
Conduct a linear regression to test the association between the recessive-coded rs1861 and the phenotype. What is the p-value?
+
+
+
#Conduct a linear regression model between rs1861 and the phenotype under recessive coding
+lm_rs1861_rec <-lm(PHENOTYPE ~ rs1861_C, data = geno_rs1861_rec)
+summary(lm_rs1861_rec)
+
+
+Call:
+lm(formula = PHENOTYPE ~ rs1861_C, data = geno_rs1861_rec)
+
+Residuals:
+ Min 1Q Median 3Q Max
+-3.5437 -0.6892 0.0015 0.7016 3.3270
+
+Coefficients:
+ Estimate Std. Error t value Pr(>|t|)
+(Intercept) -0.007686 0.049446 -0.155 0.876
+rs1861_C 1.007543 0.052242 19.286 <2e-16 ***
+---
+Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+
+Residual standard error: 1.005 on 3962 degrees of freedom
+Multiple R-squared: 0.08582, Adjusted R-squared: 0.08559
+F-statistic: 372 on 1 and 3962 DF, p-value: < 2.2e-16
+
+
#The p-value is less than 2 × 10⁻¹⁶. There is a statistically significant relationship between rs1861 and the phenotype.
+
+
+
Plot the scatterplot of phenotype versus the recessive-coded genotype of rs1861. Add the regression line to the plot.
+
+
+
# Create a new data frame containing a sequence of genotype values
+ geno_range_rec <-data.frame(
+rs1861_C =seq(0, 2, length.out =100)
+ )
+# Use the fitted linear regression model to predict phenotype probability
+# for each genotype value in geno_range
+ geno_range_rec$predicted_prob <-predict(lm_rs1861_rec, newdata = geno_range_rec, type ="response")
+
+# Start from the observed genotype dataset
+count_data_rec <- geno_A1_subset_rec %>%
+# Remove rows with missing genotype or missing phenotype
+filter(!is.na(rs1861_C), !is.na(PHENOTYPE)) %>%
+# Group data by genotype value and phenotype category
+group_by(rs1861_C, PHENOTYPE) %>%
+# Count how many observations fall into each genotype-phenotype combination
+summarise(count =n(), .groups ="drop")
+
+# Start building the plot
+p2=ggplot() +
+# Add points showing observed genotype-phenotype combinations
+geom_point(
+# Use the summarized count_data for the points
+data = count_data_rec,
+# Map genotype to x-axis, phenotype to y-axis, and count to point size
+aes(x = rs1861_C, y = PHENOTYPE, size = count),
+color ="blue"
+ ) +
+# Add a regression line showing predicted probability across genotype values
+geom_line(# Use the prediction data frame for the line
+data = geno_range_rec,
+aes(x = rs1861_C, y = predicted_prob),
+color ="red",
+# Set line thickness
+size =1.2
+ ) +
+# Control the size range of the points
+scale_size_continuous(range =c(2, 10)) +
+labs(
+title ="Linear Regression: PHENOTYPE ~ Genotype",
+x ="Genotype (0/1)",
+y ="Phenotype ",
+size ="Count") +
+# Restrict the visible x- and y-axis range
+coord_cartesian(xlim =c(-0.5,1.5)) +
+scale_x_continuous(breaks =c(0, 1))+
+# Apply a minimal theme for a clean appearance
+theme_minimal()
+
+print(p2)
# The AIC is higher in the recessive model (AIC=11291.74), indicating that the additive model (AIC=11276.91) may fit better slightly better overall by capturing a more gradual allele-phenotype effect.
+
+# However, it is also important to note that there are not a lot of individuals with the homozygous variant genotype in this dataset. As a result, that may be why the recessive model and additive models fit similarly.
+
+
+
Criteria
+
+
+
+
+
+
+
+
+
Criteria
+
Complete
+
Incomplete
+
+
+
+
+
Data Inspection
+
Correct sample/SNP counts and variable type identified.
+
Missing or incorrect counts or variable type.
+
+
+
Allele Frequency Estimation
+
Correct allele and minor allele frequencies computed.
+
Frequencies missing or wrong.
+
+
+
Hardy–Weinberg Equilibrium Test
+
Correct PLINK command and p-value extraction in R.
+
PLINK command or extraction incorrect/missing.
+
+
+
Genetic Association Test
+
Correct regressions, plots, coding, and interpretation.
+
Regression, plots, or interpretation missing/incomplete.
+
+
+
+
+
+
Submission Information
+
📌 Please review our Assignment Submission Guide for detailed instructions on how to format, branch, and submit your work. Following these guidelines is crucial for your submissions to be evaluated correctly.
+
+
Note:
+
If you like, you may collaborate with others in the cohort. If you choose to do so, please indicate with whom you have worked with in your pull request by tagging their GitHub username. Separate submissions are required.
+
+
+
+
Submission Parameters
+
+
Submission Due Date: 11:59 PM – 16/03/2026
+
Branch name for your repo should be: assignment-1
+
What to submit for this assignment:
+
+
Populate this Quarto document (assignment_1.qmd).
+
Render the document with Quarto: quarto render assignment_1.qmd.
+
Submit both assignment_1.qmd and the rendered HTML file assignment_1.html in your pull request.
+
+
What the pull request link should look like for this assignment: https://github.com/<your_github_username>/gen_data/pull/<pr_id>
+
+
Open a private window in your browser. Copy and paste the link to your pull request into the address bar. Make sure you can see your pull request properly. This helps the technical facilitator and learning support team review your submission easily.
+
+
+
+
Checklist:
+
+
Created a branch with the correct naming convention.
+
Ensured that the repository is public.
+
Reviewed the PR description guidelines and adhered to them.
+
Verified that the link is accessible in a private browser window.
+
Confirmed that both assignment_1.qmd and assignment_1.html are included in the pull request.
+
+
If you encounter any difficulties or have questions, please don’t hesitate to reach out to our team via our Slack help channel. Our technical facilitators and learning support team are here to help you navigate any challenges.
+
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/02_activities/assignments/assignment_1.qmd b/02_activities/assignments/assignment_1.qmd
index 550af3d..064ee7e 100644
--- a/02_activities/assignments/assignment_1.qmd
+++ b/02_activities/assignments/assignment_1.qmd
@@ -13,98 +13,298 @@ Please bring questions that you cannot work out on your own to office hours, wor
You will need to install PLINK and run the analyses. Please follow the OS-specific setup guide in [`SETUP.md`](../../SETUP.md). PLINK is a free, open-source whole genome association analysis toolset, designed to perform a range of basic, large-scale analyses in a computationally efficient manner.
+```{r setup}
+# Set up a consistent path for all chunks of code
+knitr::opts_knit$set(root.dir = normalizePath("../../"))
+
+library(data.table)
+library(ggplot2)
+library(seqminer)
+library(HardyWeinberg)
+library(dplyr)
+```
+
+
#### Question 1: Data inspection
Before fitting any models, it is essential to understand the data. Use R or bash code to answer the following questions about the `gwa.qc.A1.fam`, `gwa.qc.A1.bim`, and `gwa.qc.A1.bed` files, available at the following Google Drive link: . Please download all three files from this link and place them in `02_activities/data/`.
(i) Read the .fam file. How many samples does the dataset contain?
-```
-# Your answer here...
+```{bash}
+#Read the .fam file
+wc -l ./02_activities/data/gwa.qc.A1.fam
+
+#The dataset contains 4000 samples.
```
(ii) What is the 'variable type' of the response variable (i.e.Continuous or binary)?
-```
-# Your answer here...
+```{bash}
+head ./02_activities/data/gwa.qc.A1.fam
+
+#The response variable is continuous.
```
(iii) Read the .bim file. How many SNPs does the dataset contain?
-```
-# Your answer here...
+```{bash}
+wc -l ./02_activities/data/gwa.qc.A1.bim
+
+#There are 101083 SNPs in the dataset.
```
#### Question 2: Allele Frequency Estimation
(i) Load the genotype matrix for SNPs rs1861, rs3813199, rs3128342, and rs11804831 using additive coding. What are the allele frequencies (AFs) for these four SNPs?
-```
-# Your code here...
+```{bash}
+#Create SNP list file with 4 SNPs: rs1861, rs3813199, rs3128342, and rs11804831
+setopt interactivecomments
+printf "%s\n" rs1861 rs3813199 rs3128342 rs11804831 > ./02_activities/data/A1snplist.txt
+
+#Subset SNPs and output allele frequencies
+plink2 --bfile ./02_activities/data/gwa.qc.A1 \
+ --extract ./02_activities/data/A1snplist.txt \
+ --recode A\
+ --freq \
+ --out ./02_activities/data/gwa_qc_A1_freq
+```
+```{r}
+#View allele frequencies
+A1freq <- fread("./02_activities/data/gwa_qc_A1_freq.afreq")
+head(A1freq)
+
+#The allele frequencies for the alternative allele of rs1861, rs3813199, rs3128342, and rs11804831 are as follows: 0.0539859, 0.0569126, 0.3051210, and 0.1543410.
+
+A1freq$REF_FREQ <- 1 - A1freq$ALT_FREQS
+A1freq
+
+#The allele frequencies for the reference allele of rs1861, rs3813199, rs3128342, and rs11804831 are as follows: 0.9460141, 0.9430874, 0.6948790, and 0.8456590.
```
(ii) What are the minor allele frequencies (MAFs) for these four SNPs?
-```
-# Your code here...
+```{r}
+A1freq$MAF <- pmin(A1freq$ALT_FREQS, 1 - A1freq$ALT_FREQS)
+
+A1freq[, .(ID, AF = ALT_FREQS, MAF)]
+
+#The minor allele frequencies for rs1861, rs3813199, rs3128342, and rs11804831 are as follows: 0.0539859, 0.0569126, 0.3051210, and 0.1543410.
```
+
#### Question 3: Hardy–Weinberg Equilibrium (HWE) Test
(i) Conduct the Hardy–Weinberg Equilibrium (HWE) test for all SNPs in the .bim file. Then, load the file containing the HWE p-value results and display the first few rows of the resulting data frame.
-```
-# Your code here...
+```{bash}
+plink2 --bfile ./02_activities/data/gwa.qc.A1 --hardy --out ./02_activities/data/gwa_qc_A1_hwe
+```
+
+```{r}
+A1hwe <- fread("./02_activities/data/gwa_qc_A1_hwe.hardy")
+head(A1hwe)
```
(ii) What are the HWE p-values for SNPs rs1861, rs3813199, rs3128342, and rs11804831?
-```
-# Your code here...
+```{bash}
+#Subset SNPs and output allele frequencies
+plink2 --bfile ./02_activities/data/gwa.qc.A1 \
+ --extract ./02_activities/data/A1snplist.txt \
+ --hardy \
+ --out ./02_activities/data/gwa_qc_A1_hwe_snplist
+```
+
+```{r}
+#View allele frequencies
+A1hwesnplist <- fread("./02_activities/data/gwa_qc_A1_hwe_snplist.hardy")
+head(A1hwesnplist)
+
+#The HWE p-values for rs1861, rs3813199, rs3128342, and rs11804831 are the following: 0.274719, 1.000000, 0.330273, and 0.113354.
```
#### Question 4: Genetic Association Test
(i) Conduct a linear regression to test the association between SNP rs1861 and the phenotype. What is the p-value?
-```
-# Your code here...
+```{bash}
+plink2 --bfile ./02_activities/data/gwa.qc.A1 \
+ --extract ./02_activities/data/A1snplist.txt \
+ --recode A \
+ --out ./02_activities/data/geno.A1.additive
+```
+
+```{r}
+# Load genotype data
+geno_A1_subset <- fread("./02_activities/data/geno.A1.additive.raw")
+geno_A1_subset=geno_A1_subset[!is.na(rs1861_C), ]
+geno_A1_subset$PHENOTYPE=geno_A1_subset$PHENOTYPE-1
+head(geno_A1_subset)
+
+geno_rs1861 <- geno_A1_subset %>%
+ select(FID, IID, PHENOTYPE, rs1861_C) %>%
+ na.omit()
+
+#Conduct a linear regression model between rs1861 and the phenotype
+lm_rs1861 <- lm(PHENOTYPE ~ rs1861_C, data = geno_rs1861)
+summary(lm_rs1861)
+
+#The p-value is less than 2 × 10⁻¹⁶. There is a statistically significant relationship between rs1861 and the phenotype.
```
(ii) How would you interpret the beta coefficient from this regression?
```
-# Your answer here...
+We performed a linear regression to test the association between genotype rs1861 and the continuous phenotype. The beta-coefficient is 0.97382. This means that for every additional C allele an individual has for rs1861, their phenotype increases by 0.97382 units.
```
(iii) Plot the scatterplot of phenotype versus the genotype of SNP rs1861. Add the regression line to the plot.
-```
-# Your code here...
+```{r}
+# Create a new data frame containing a sequence of genotype values
+ geno_range <- data.frame(
+ rs1861_C = seq(0, 2, length.out = 100)
+ )
+# Use the fitted linear regression model to predict phenotype probability
+# for each genotype value in geno_range
+ geno_range$predicted_prob <- predict(lm_rs1861, newdata = geno_range, type = "response")
+
+# Start from the observed genotype dataset
+count_data <- geno_A1_subset %>%
+# Remove rows with missing genotype or missing phenotype
+ filter(!is.na(rs1861_C), !is.na(PHENOTYPE)) %>%
+ # Group data by genotype value and phenotype category
+ group_by(rs1861_C, PHENOTYPE) %>%
+# Count how many observations fall into each genotype-phenotype combination
+ summarise(count = n(), .groups = "drop")
+
+# Start building the plot
+p1=ggplot() +
+# Add points showing observed genotype-phenotype combinations
+ geom_point(
+ # Use the summarized count_data for the points
+ data = count_data,
+ # Map genotype to x-axis, phenotype to y-axis, and count to point size
+ aes(x = rs1861_C, y = PHENOTYPE, size = count),
+ color = "blue"
+ ) +
+ # Add a regression line showing predicted probability across genotype values
+ geom_line(# Use the prediction data frame for the line
+ data = geno_range,
+ aes(x = rs1861_C, y = predicted_prob),
+ color = "red",
+ # Set line thickness
+ size = 1.2
+ ) +
+ # Control the size range of the points
+ scale_size_continuous(range = c(2, 10)) +
+ labs(
+ title = "Linear Regression: PHENOTYPE ~ Genotype",
+ x = "Genotype (0/1/2)",
+ y = "Phenotype ",
+ size = "Count"
+ ) +
+ # Restrict the visible x- and y-axis range
+ coord_cartesian(xlim = c(-0.5, 2.5)) +
+ # Apply a minimal theme for a clean appearance
+ theme_minimal()
+
+print(p1)
```
(iv) Convert the genotype coding for rs1861 to recessive coding.
-```
-# Your code here...
+``` {r}
+geno_A1_subset_rec <- geno_A1_subset # make a copy
+
+#change into recessive coding
+geno_A1_subset_rec[, 7:ncol(geno_A1_subset_rec)] <- as.data.frame(
+ lapply(geno_A1_subset[, 7:ncol(geno_A1_subset)], function(x) ifelse(x == 2, 1, 0))
+)
+
+#Subset to just rs1861
+geno_rs1861_rec <- geno_A1_subset_rec %>%
+ select(FID, IID, PHENOTYPE, rs1861_C) %>%
+ na.omit()
+
```
(v) Conduct a linear regression to test the association between the recessive-coded rs1861 and the phenotype. What is the p-value?
-```
-# Your code here...
+``` {r}
+#Conduct a linear regression model between rs1861 and the phenotype under recessive coding
+lm_rs1861_rec <- lm(PHENOTYPE ~ rs1861_C, data = geno_rs1861_rec)
+summary(lm_rs1861_rec)
+
+#The p-value is less than 2 × 10⁻¹⁶. There is a statistically significant relationship between rs1861 and the phenotype.
```
(vi) Plot the scatterplot of phenotype versus the recessive-coded genotype of rs1861. Add the regression line to the plot.
-```
-# Your code here...
+```{r}
+# Create a new data frame containing a sequence of genotype values
+ geno_range_rec <- data.frame(
+ rs1861_C = seq(0, 2, length.out = 100)
+ )
+# Use the fitted linear regression model to predict phenotype probability
+# for each genotype value in geno_range
+ geno_range_rec$predicted_prob <- predict(lm_rs1861_rec, newdata = geno_range_rec, type = "response")
+
+# Start from the observed genotype dataset
+count_data_rec <- geno_A1_subset_rec %>%
+# Remove rows with missing genotype or missing phenotype
+ filter(!is.na(rs1861_C), !is.na(PHENOTYPE)) %>%
+ # Group data by genotype value and phenotype category
+ group_by(rs1861_C, PHENOTYPE) %>%
+# Count how many observations fall into each genotype-phenotype combination
+ summarise(count = n(), .groups = "drop")
+
+# Start building the plot
+p2=ggplot() +
+# Add points showing observed genotype-phenotype combinations
+ geom_point(
+ # Use the summarized count_data for the points
+ data = count_data_rec,
+ # Map genotype to x-axis, phenotype to y-axis, and count to point size
+ aes(x = rs1861_C, y = PHENOTYPE, size = count),
+ color = "blue"
+ ) +
+ # Add a regression line showing predicted probability across genotype values
+ geom_line(# Use the prediction data frame for the line
+ data = geno_range_rec,
+ aes(x = rs1861_C, y = predicted_prob),
+ color = "red",
+ # Set line thickness
+ size = 1.2
+ ) +
+ # Control the size range of the points
+ scale_size_continuous(range = c(2, 10)) +
+ labs(
+ title = "Linear Regression: PHENOTYPE ~ Genotype",
+ x = "Genotype (0/1)",
+ y = "Phenotype ",
+ size = "Count") +
+ # Restrict the visible x- and y-axis range
+ coord_cartesian(xlim = c(-0.5,1.5)) +
+ scale_x_continuous(breaks = c(0, 1))+
+ # Apply a minimal theme for a clean appearance
+ theme_minimal()
+
+print(p2)
```
(vii) Which model fits better? Justify your answer.
+```{r}
+AIC(lm_rs1861, lm_rs1861_rec)
+```
+
```
-# Your answer here...
+# The AIC is higher in the recessive model (AIC=11291.74), indicating that the additive model (AIC=11276.91) may fit better slightly better overall by capturing a more gradual allele-phenotype effect.
+
+# However, it is also important to note that there are not a lot of individuals with the homozygous variant genotype in this dataset. As a result, that may be why the recessive model and additive models fit similarly.
```
### Criteria
diff --git a/02_activities/assignments/assignment_2.html b/02_activities/assignments/assignment_2.html
new file mode 100644
index 0000000..057a085
--- /dev/null
+++ b/02_activities/assignments/assignment_2.html
@@ -0,0 +1,1018 @@
+
+
+
+
+
+
+
+
+
+Assignment #2
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Assignment #2
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
You will need PLINK2 installed and available in your PATH. Please follow the OS-specific setup guide in SETUP.md. The dataset for this assignment consists of the following binary PLINK files: gwa.A2.bed, gwa.A2.bim, gwa.A2.fam , available at the following Google Drive link: https://drive.google.com/drive/folders/1rHoy3z52Yyj985ukjjtLhfIBxchpyYtZ?usp=drive_link. Please download all three files and save them in 02_activities/data/.
+
+
Question 1: Data inspection
+
Before you run any models, first get familiar with the dataset. You may find data.table::fread() in R helpful for reading .bim and .fam files.
+
+
Read the .fam file. How many samples does the dataset contain?
+
+
+
#Read the .fam file
+
+fam <-fread("./02_activities/data/gwa.A2.fam", header =FALSE)
+nrow(fam)
+
+
[1] 4000
+
+
#The dataset contains 4000 samples.
+
+
+
Read the .bim file. How many SNPs does the dataset contain?
+
+
+
#Read the .bim file
+bim <-fread("./02_activities/data/gwa.A2.bim", header =FALSE)
+nrow(bim)
+
+
[1] 306102
+
+
#The dataset contains 306102 SNPs.
+
+
+
Make a histogram (or density plot) of the phenotype. Does it look roughly Gaussian?
#Create a histogram
+hist(fam$PHENOTYPE,
+main ="Histogram of Phenotype",
+xlab ="PHENOTYPE",
+breaks =30)
+
+
+
+
+
+
+
+
#Yes, the distribution of the phenotype looks roughly Gaussian.
+
+
+
+
Question 2: Quality control (QC)
+
Now we will perform QC using PLINK2 for the genotype files in gwa.A2.
+
+
Using PLINK2 from the command line (bash), perform basic QC with the following filters: MAF ≥ 0.05, SNP missingness (--geno) ≤ 0.01, individual missingness (--mind) ≤ 0.10, and HWE p-value ≥ 0.00005, and output the QC’ed dataset as gwa.qc.A2.
PLINK v2.00a5.12 M1 (25 Jun 2024) www.cog-genomics.org/plink/2.0/
+(C) 2005-2024 Shaun Purcell, Christopher Chang GNU General Public License v3
+Logging to ./02_activities/data/gwa.qc.A2.log.
+Options in effect:
+ --bfile ./02_activities/data/gwa.A2
+ --geno 0.01
+ --hwe 0.00005
+ --maf 0.05
+ --make-bed
+ --mind 0.1
+ --out ./02_activities/data/gwa.qc.A2
+
+Start time: Tue Mar 31 14:11:15 2026
+16384 MiB RAM detected; reserving 8192 MiB for main workspace.
+Using up to 10 threads (change this with --threads).
+4000 samples (2000 females, 2000 males; 4000 founders) loaded from
+./02_activities/data/gwa.A2.fam.
+306102 variants loaded from ./02_activities/data/gwa.A2.bim.
+1 quantitative phenotype loaded (4000 values).
+Calculating sample missingness rates... 0%21%42%64%85%done.
+0 samples removed due to missing genotype data (--mind).
+4000 samples (2000 females, 2000 males; 4000 founders) remaining after main
+filters.
+4000 quantitative phenotype values remaining after main filters.
+Calculating allele frequencies... 0%21%42%64%85%done.
+--geno: 196578 variants removed due to missing genotype data.
+--hwe: 6 variants removed due to Hardy-Weinberg exact test (founders only).
+8435 variants removed due to allele frequency threshold(s)
+(--maf/--max-maf/--mac/--max-mac).
+101083 variants remaining after main filters.
+Writing ./02_activities/data/gwa.qc.A2.fam ... done.
+Writing ./02_activities/data/gwa.qc.A2.bim ... done.
+Writing ./02_activities/data/gwa.qc.A2.bed ... 0%21%43%65%87%done.
+End time: Tue Mar 31 14:11:16 2026
+
+
+
+
+
Question 3: Relatedness
+
In this question, you will use PLINK2’s built-in KING-robust kinship (--king-cutoff) detect and remove related individuals.
+
+
Perform LD pruning on gwa.qc.A2 using PLINK2 with the following parameters: --indep-pairwise 500 50 0.05, and then generate a new dataset containing only the pruned SNPs.
+
+
+
#Create a subset of approximately independent SNPs
+plink2--bfile ./02_activities/data/gwa.qc.A2 --indep-pairwise 500 50 0.05 --out ./02_activities/data/gwa.qc.A2_pruned
+
+#Build a new dataset using only pruned SNPs
+plink2--bfile ./02_activities/data/gwa.qc.A2 --extract ./02_activities/data/gwa.qc.A2_pruned.prune.in --make-bed--out ./02_activities/data/gwa.qc.A2_pruned
+
+
PLINK v2.00a5.12 M1 (25 Jun 2024) www.cog-genomics.org/plink/2.0/
+(C) 2005-2024 Shaun Purcell, Christopher Chang GNU General Public License v3
+Logging to ./02_activities/data/gwa.qc.A2_pruned.log.
+Options in effect:
+ --bfile ./02_activities/data/gwa.qc.A2
+ --indep-pairwise 500 50 0.05
+ --out ./02_activities/data/gwa.qc.A2_pruned
+
+Start time: Tue Mar 31 14:11:16 2026
+16384 MiB RAM detected; reserving 8192 MiB for main workspace.
+Using up to 10 threads (change this with --threads).
+4000 samples (2000 females, 2000 males; 4000 founders) loaded from
+./02_activities/data/gwa.qc.A2.fam.
+101083 variants loaded from ./02_activities/data/gwa.qc.A2.bim.
+1 quantitative phenotype loaded (4000 values).
+Calculating allele frequencies... 0%64%done.
+--indep-pairwise (9 compute threads): 0%50%79169/101083 variants removed.
+Writing...
+Variant lists written to ./02_activities/data/gwa.qc.A2_pruned.prune.in and
+./02_activities/data/gwa.qc.A2_pruned.prune.out .
+End time: Tue Mar 31 14:11:16 2026
+PLINK v2.00a5.12 M1 (25 Jun 2024) www.cog-genomics.org/plink/2.0/
+(C) 2005-2024 Shaun Purcell, Christopher Chang GNU General Public License v3
+Logging to ./02_activities/data/gwa.qc.A2_pruned.log.
+Options in effect:
+ --bfile ./02_activities/data/gwa.qc.A2
+ --extract ./02_activities/data/gwa.qc.A2_pruned.prune.in
+ --make-bed
+ --out ./02_activities/data/gwa.qc.A2_pruned
+
+Start time: Tue Mar 31 14:11:16 2026
+16384 MiB RAM detected; reserving 8192 MiB for main workspace.
+Using up to 10 threads (change this with --threads).
+4000 samples (2000 females, 2000 males; 4000 founders) loaded from
+./02_activities/data/gwa.qc.A2.fam.
+101083 variants loaded from ./02_activities/data/gwa.qc.A2.bim.
+1 quantitative phenotype loaded (4000 values).
+--extract: 21914 variants remaining.
+21914 variants remaining after main filters.
+Writing ./02_activities/data/gwa.qc.A2_pruned.fam ... done.
+Writing ./02_activities/data/gwa.qc.A2_pruned.bim ... done.
+Writing ./02_activities/data/gwa.qc.A2_pruned.bed ... 0%61%done.
+End time: Tue Mar 31 14:11:16 2026
+
+
+
+
Use PLINK2 on the LD-pruned dataset to identify a set of unrelated individuals up to (approximately) 2nd-degree relatives (use a kinship cutoff of 0.0884).
PLINK v2.00a5.12 M1 (25 Jun 2024) www.cog-genomics.org/plink/2.0/
+(C) 2005-2024 Shaun Purcell, Christopher Chang GNU General Public License v3
+Logging to ./02_activities/data/gwa.qc_assoc.A2.log.
+Options in effect:
+ --adjust
+ --bfile ./02_activities/data/gwa.qc.A2
+ --ci 0.95
+ --covar ./02_activities/data/gwa_unrel_PCA.A2.eigenvec
+ --covar-name PC1,PC2,PC3
+ --glm
+ --out ./02_activities/data/gwa.qc_assoc.A2
+
+Start time: Tue Mar 31 14:11:22 2026
+16384 MiB RAM detected; reserving 8192 MiB for main workspace.
+Using up to 10 threads (change this with --threads).
+4000 samples (2000 females, 2000 males; 4000 founders) loaded from
+./02_activities/data/gwa.qc.A2.fam.
+101083 variants loaded from ./02_activities/data/gwa.qc.A2.bim.
+1 quantitative phenotype loaded (4000 values).
+3 covariates loaded from ./02_activities/data/gwa_unrel_PCA.A2.eigenvec.
+Calculating allele frequencies... 0%64%done.
+--glm linear regression on phenotype 'PHENO1': 0%64%done.
+Results written to ./02_activities/data/gwa.qc_assoc.A2.PHENO1.glm.linear .
+--adjust: Genomic inflation est. lambda (based on median chisq) = 1.01779.
+--adjust values (101083 tests) written to
+./02_activities/data/gwa.qc_assoc.A2.PHENO1.glm.linear.adjusted .
+End time: Tue Mar 31 14:11:23 2026
+
+
+
+
Create a Manhattan plot of the GWAS results.
+
+
+
# Load PLINK2 linear results
+assocA2 <- data.table::fread("./02_activities/data/gwa.qc_assoc.A2.PHENO1.glm.linear")
+
+head(assocA2)
+
+
#CHROM POS ID REF ALT PROVISIONAL_REF? A1 OMITTED
+ <int> <int> <char> <char> <char> <char> <char> <char>
+1: 1 1011278 rs3737728 G A Y A G
+2: 1 1011278 rs3737728 G A Y A G
+3: 1 1011278 rs3737728 G A Y A G
+4: 1 1011278 rs3737728 G A Y A G
+5: 1 1109721 rs1320565 C T Y T C
+6: 1 1109721 rs1320565 C T Y T C
+ A1_FREQ TEST OBS_CT BETA SE L95 U95 T_STAT
+ <num> <char> <int> <num> <num> <num> <num> <num>
+1: 0.3386490 ADD 3982 0.0182086 0.0240795 -0.0289864 0.0654035 0.756187
+2: 0.3386490 PC1 3982 -0.3779240 1.0015500 -2.3409200 1.5850800 -0.377340
+3: 0.3386490 PC2 3982 0.9818050 1.0015100 -0.9811150 2.9447300 0.980326
+4: 0.3386490 PC3 3982 0.4832040 1.0015300 -1.4797500 2.4461600 0.482467
+5: 0.0788481 ADD 3976 -0.0153297 0.0420641 -0.0977739 0.0671145 -0.364436
+6: 0.0788481 PC1 3976 -0.4762890 1.0037500 -2.4436000 1.4910200 -0.474511
+ P ERRCODE
+ <num> <char>
+1: 0.449582 .
+2: 0.705941 .
+3: 0.326985 .
+4: 0.629501 .
+5: 0.715552 .
+6: 0.635162 .
+
+
# Keep only the additive test (one row per SNP)
+assoc_addA2 <- assocA2[TEST =="ADD"]
+
+# qqman expects columns named CHR, BP, SNP, and P.
+# In PLINK2 output they are #CHROM, POS, ID, and P.
+setnames(assoc_addA2, c("#CHROM","POS","ID"), c("CHR","BP","SNP"))
+
+png("./02_activities/assignments/manhattan_plot.png", width =1200, height =800, res =150)
+# Manhattan plot
+manhattan(assoc_addA2,
+chr ="CHR",
+bp ="BP",
+snp ="SNP",
+p ="P",
+xlab ="Chromosome",
+ylab ="-log10(P)",
+suggestiveline =FALSE,
+cex.axis =1.5,
+col =c("lightblue", "lightslateblue"))
+dev.off()
+
+
quartz_off_screen
+ 2
+
+
+
+
Create a QQ plot of the GWAS p-values.
+
+
+
png("./02_activities/assignments/QQ_plot.png", width =1200, height =800, res =150)
+# Q-Q plot
+qq(assoc_addA2$P, main ="Q-Q plot of GWAS p-values")
+dev.off()
+
+
quartz_off_screen
+ 2
+
+
+
+
+
Criteria
+
+
+
+
+
+
+
+
+
Criteria
+
Complete
+
Incomplete
+
+
+
+
+
Data inspection
+
Correct sample and SNP counts; phenotype plot with a brief comment on Gaussianity.
+
Counts or phenotype description/plot missing or incorrect.
+
+
+
QC & LD pruning
+
Correct PLINK2 QC command and thresholds.
+
QC/pruning commands, thresholds, or output datasets missing or incorrect.
+
+
+
Relatedness & PCA
+
Correct use of PLINK2 command to obtain unrelated samples and PCA run on pruned SNPs.
+
Relatedness step, unrelated dataset, or PCA analysis missing or incorrect.
+
+
+
GWAS & visualisation
+
Linear regression GWAS with PCs as covariates; Manhattan and QQ plots produced.
+
GWAS command, or Manhattan/QQ plots missing or clearly incorrect.
+
+
+
+
+
+
Submission Information
+
🚨 Please review our Assignment Submission Guide 🚨 for detailed instructions on how to format, branch, and submit your work. Following these guidelines is crucial for your submissions to be evaluated correctly.
+
+
+
Submission Parameters
+
+
Submission Due Date: 11:59 PM – 01/04/2026
+
The branch name for your repo should be: assignment-2
+
What to submit for this assignment:
+
+
Populate this Quarto document (assignment_2.qmd).
+
Render the document with Quarto: quarto render assignment_2.qmd.
+
Submit both assignment_2.qmd and the rendered HTML file assignment_2.html in your pull request.
+
+
What the pull request link should look like for this assignment: https://github.com/<your_github_username>/gen_data/pull/<pr_id>
+
+
Open a private window in your browser. Copy and paste the link to your pull request into the address bar. Make sure you can see your pull request properly. This helps the technical facilitator and learning support team review your submission easily.
Verify that your link is accessible in a private browser window.
+
Confirm that both assignment_2.qmd and assignment_2.html are included in the pull request.
+
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/02_activities/assignments/assignment_2.qmd b/02_activities/assignments/assignment_2.qmd
index fd93722..a5399bf 100644
--- a/02_activities/assignments/assignment_2.qmd
+++ b/02_activities/assignments/assignment_2.qmd
@@ -5,6 +5,14 @@ format: html
You will need **PLINK2** installed and available in your PATH. Please follow the OS-specific setup guide in [`SETUP.md`](../../SETUP.md). The dataset for this assignment consists of the following binary PLINK files: `gwa.A2.bed`, `gwa.A2.bim`, `gwa.A2.fam` , available at the following Google Drive link: . Please download all three files and save them in `02_activities/data/`.
+```{r setup, include=FALSE}
+# Input files are expected in ./02_activities/data/.
+#knitr::opts_chunk$set(echo = TRUE)
+knitr::opts_knit$set(root.dir = normalizePath("../../"))
+library(data.table)
+library(qqman)
+```
+
#### Question 1: Data inspection
Before you run any models, first get familiar with the dataset. You may find `data.table::fread()` in R helpful for reading `.bim` and `.fam` files.
@@ -12,19 +20,45 @@ Before you run any models, first get familiar with the dataset. You may find `da
(i) Read the `.fam` file. How many samples does the dataset contain?
```{r}
-# Your code here...
+#Read the .fam file
+
+fam <- fread("./02_activities/data/gwa.A2.fam", header = FALSE)
+nrow(fam)
+
+#The dataset contains 4000 samples.
```
(ii) Read the `.bim` file. How many SNPs does the dataset contain?
```{r}
-# Your code here...
+#Read the .bim file
+bim <- fread("./02_activities/data/gwa.A2.bim", header = FALSE)
+nrow(bim)
+
+#The dataset contains 306102 SNPs.
```
(iii) Make a histogram (or density plot) of the phenotype. Does it look roughly Gaussian?
```{r}
-# Your code here...
+#Check the header
+head(fam)
+
+#Name the columns
+colnames(fam) <- c("FID","IID","PID","MID","SEX","PHENOTYPE")
+
+#Double check that worked and it corresponds
+names(fam)
+head(fam)
+
+#Create a histogram
+hist(fam$PHENOTYPE,
+ main = "Histogram of Phenotype",
+ xlab = "PHENOTYPE",
+ breaks = 30)
+
+#Yes, the distribution of the phenotype looks roughly Gaussian.
+
```
#### Question 2: Quality control (QC)
@@ -34,7 +68,7 @@ Now we will perform QC using PLINK2 for the genotype files in `gwa.A2`.
(i) Using PLINK2 from the command line (bash), perform basic QC with the following filters: MAF ≥ 0.05, SNP missingness (`--geno`) ≤ 0.01, individual missingness (`--mind`) ≤ 0.10, and HWE p-value ≥ 0.00005, and output the QC’ed dataset as `gwa.qc.A2`.
```{bash}
-# Your code here...
+plink2 --bfile ./02_activities/data/gwa.A2 --maf 0.05 --geno 0.01 --mind 0.1 --hwe 0.00005 --make-bed --out ./02_activities/data/gwa.qc.A2
```
#### Question 3: Relatedness
@@ -44,19 +78,25 @@ In this question, you will use **PLINK2’s built-in KING-robust kinship** (`--k
i) Perform LD pruning on `gwa.qc.A2` using PLINK2 with the following parameters: `--indep-pairwise 500 50 0.05`, and then generate a new dataset containing only the pruned SNPs.
```{bash}
-# Your code here...
+#Create a subset of approximately independent SNPs
+plink2 --bfile ./02_activities/data/gwa.qc.A2 --indep-pairwise 500 50 0.05 --out ./02_activities/data/gwa.qc.A2_pruned
+
+#Build a new dataset using only pruned SNPs
+plink2 --bfile ./02_activities/data/gwa.qc.A2 --extract ./02_activities/data/gwa.qc.A2_pruned.prune.in --make-bed --out ./02_activities/data/gwa.qc.A2_pruned
```
(ii) Use PLINK2 on the LD-pruned dataset to identify a set of unrelated individuals up to (approximately) 2nd-degree relatives (use a kinship cutoff of 0.0884).
```{bash}
-# Your code here...
+
+plink2 --bfile ./02_activities/data/gwa.qc.A2_pruned --king-cutoff 0.0884 --out ./02_activities/data/gwa.A2_king
+
```
(iii) Using the unrelated individual list from part (ii), create a dataset containing only unrelated individuals and name it `gwa.qc_unrelated.A2`.
```{bash}
-# Your code here...
+plink2 --bfile ./02_activities/data/gwa.qc.A2 --keep ./02_activities/data/gwa.A2_king.king.cutoff.in.id --make-bed --out ./02_activities/data/gwa.qc_unrelated.A2
```
#### Question 4: principal component analysis (PCA)
@@ -66,7 +106,9 @@ i) Perform LD pruning on `gwa.qc.A2` using PLINK2 with the following parameters
(Hint: Refer to the PCA section in the GWAS Tutorial II. You should use a `*.prune.in` file to select a set of independent SNPs before performing PCA.)
```{bash}
-# Your code here...
+plink2 --bfile ./02_activities/data/gwa.qc_unrelated.A2 --extract ./02_activities/data/gwa.qc.A2_pruned.prune.in --make-bed --out ./02_activities/data/gwa.qc_unrelated_pruned.A2
+
+plink2 --bfile ./02_activities/data/gwa.qc_unrelated_pruned.A2 --pca 20 --out ./02_activities/data/gwa_unrel_PCA.A2
```
#### Question 5: GWAS analyses
@@ -82,19 +124,46 @@ Assume that:
(i) Using PLINK2, run a linear regression GWAS using `gwa.qc_unrelated.A2` as the input dataset. Adjust for PCs 1–3 as covariates.
```{bash}
-# Your code here...
+plink2 --bfile ./02_activities/data/gwa.qc.A2 --ci 0.95 --logistic --covar ./02_activities/data/gwa_unrel_PCA.A2.eigenvec --covar-name PC1,PC2,PC3 --adjust --out ./02_activities/data/gwa.qc_assoc.A2
```
(ii) Create a Manhattan plot of the GWAS results.
```{r}
-# Your code here...
+# Load PLINK2 linear results
+assocA2 <- data.table::fread("./02_activities/data/gwa.qc_assoc.A2.PHENO1.glm.linear")
+
+head(assocA2)
+
+# Keep only the additive test (one row per SNP)
+assoc_addA2 <- assocA2[TEST == "ADD"]
+
+# qqman expects columns named CHR, BP, SNP, and P.
+# In PLINK2 output they are #CHROM, POS, ID, and P.
+setnames(assoc_addA2, c("#CHROM","POS","ID"), c("CHR","BP","SNP"))
+
+png("./02_activities/assignments/manhattan_plot.png", width = 1200, height = 800, res = 150)
+# Manhattan plot
+manhattan(assoc_addA2,
+ chr = "CHR",
+ bp = "BP",
+ snp = "SNP",
+ p = "P",
+ xlab = "Chromosome",
+ ylab = "-log10(P)",
+ suggestiveline = FALSE,
+ cex.axis = 1.5,
+ col = c("lightblue", "lightslateblue"))
+dev.off()
```
(iii) Create a QQ plot of the GWAS p-values.
```{r}
-# Your code here...
+png("./02_activities/assignments/QQ_plot.png", width = 1200, height = 800, res = 150)
+# Q-Q plot
+qq(assoc_addA2$P, main = "Q-Q plot of GWAS p-values")
+dev.off()
```
### Criteria
diff --git a/02_activities/assignments/manhattan_plot.png b/02_activities/assignments/manhattan_plot.png
new file mode 100644
index 0000000..71b3b30
Binary files /dev/null and b/02_activities/assignments/manhattan_plot.png differ
