############################################################
# Phenotypic QC practical - instructor/helper and stretch activities
# Instructor helper script: pre-clean the raw dataset
#
# Activities 4-12 are optional extensions that continue from the
# CORE worksheet Activities 1-3.
#
# Run this once on coriander_gwas_pheno_raw.csv before the tutorial.
# It performs the cleaning steps necessary to answer Activities 1-3 of the
# full 20-minute script (Phenotypic-CORE-worksheet.R) and writes
# a cleaned CSV that students load at the top of Activity 1.
#
#
# Inputs (expected in working directory):
#   coriander_gwas_pheno_raw.csv
#
# Outputs:
#   coriander_gwas_pheno_clean_for_20min_v1.1.csv  # cleaned student dataset
#   pheno_qc_prep_log.txt                     # what the script did + counts
#
# The cleaned dataset preserves the flag_problem_CT / flag_problem_CTS
# columns that Activity 3 in the student script expects to filter on.
############################################################

############################################################
# Variable guide for the mock coriander GWAS practical
############################################################

# Core ID and eligibility variables:
#
# participant_id:
#   Unique participant identifier.
#   Used to check for duplicate or repeated records and to match phenotype,
#   covariate and genotype data.
#
# ever_eaten_coriander:
#   Eligibility variable.
#   1 = participant has eaten coriander/cilantro.
#   0 = participant has never eaten coriander/cilantro.
#   Participants coded 0 should be excluded before phenotype analysis because
#   coriander taste phenotypes are not interpretable for people who have never
#   eaten coriander/cilantro.

# Phenotype variables:
#
# CT:
#   Binary coriander/cilantro soapy taste phenotype.
#   1 = participant reports that coriander/cilantro tastes soapy.
#   0 = participant reports that coriander/cilantro does not taste soapy.
#
# CTS:
#   Continuous strength of soapy coriander/cilantro taste.
#   Expected range: 0 to 10.
#   0 = not at all soapy.
#   10 = completely soapy.

# Minimum GWAS covariates:
#
# sex:
#   Participant sex.
#   1 = female.
#   0 = male.
#
# age:
#   Participant age in years.
#
# PC1 to PC10:
#   First 10 genetic principal components.
#   Included to adjust for genetic ancestry/population structure.
#
# genotyping_array:
#   Genotyping array or platform.
#   If there is more than one array, include k - 1 dummy variables or equivalent
#   array adjustment in the GWAS model.
#
# birth_year:
#   Participant birth year.
#   Used to check whether birth cohort adjustment is needed.
#
# birth_cohort:
#   Derived birth cohort group.
#   Used when participants span a wide range of birth years.

# Additional technical and teaching variables:
#
# batch:
#   Simplified genotyping batch variable.
#   Can be considered as an additional technical covariate.
#
# record_wave:
#   Indicates data collection wave or repeated-record version.
#   Used to decide which record to keep when a participant appears more than once.
#
# genetic_ancestry:
#   Broad genetic ancestry group.
#   Used to check sample size and phenotype availability by ancestry group.

# Main exclusion criterion:
#   Exclude participants who have never eaten coriander/cilantro.

# Minimum covariates for the mock GWAS:
#   sex, age, PC1 to PC10, genotyping_array, batch, and birth_cohort.

# Note:
#   participant_id is not a GWAS covariate. It is used for matching and checking.
#   record_wave is not a GWAS covariate. It is used for repeated-record handling.
#   genetic_ancestry is used here for sample checks; genetic PCs are the main
#   covariates used to adjust for population structure.
############################################################

############################################################
# Setup and analysis plan
############################################################

# Load packages used in this practical.
suppressPackageStartupMessages({
  library(dplyr)
  library(tidyr)
  library(ggplot2)
  library(mice)
})

# ---- Variable names ------------------------------------------------------
id_var       <- "participant_id"
eaten_var    <- "ever_eaten_coriander"
ct_var       <- "CT"
cts_var      <- "CTS"
sex_var      <- "sex"
age_var      <- "age"
pc_vars      <- paste0("PC", 1:10)
array_var    <- "genotyping_array"
birth_var    <- "birth_year"
cohort_var   <- "birth_cohort"
ancestry_var <- "genetic_ancestry"
wave_var     <- "record_wave"
batch_var    <- "batch"

# ---- Logging helper ------------------------------------------------------
log_lines <- character(0)
log_msg <- function(...) {
  msg <- paste0(format(Sys.time(), "%H:%M:%S"), "  ", paste0(..., collapse = ""))
  log_lines <<- c(log_lines, msg)
  message(msg)
}

# ---- Read raw data with import-time missing-value recoding -----------

common_missing_codes <- c(
  "", "NA", "N/A", "na", "n/a",
  "-9", "-99", "-999",
  "999", "9999",
  "missing", "Missing",
  "unknown", "Unknown",
  "prefer not to say"
)

dat_clean <- read.csv(
  file = "coriander_gwas_pheno_raw.csv",
  stringsAsFactors = FALSE,
  na.strings = common_missing_codes,
  strip.white = TRUE
)

log_msg("Read raw dataset: ", nrow(dat_clean), " rows, ", ncol(dat_clean), " columns.")

# ---- Activity 4: variable-specific invalid-value cleanup ------------
#  CTS must lie in [0, 10].
dat_clean$invalid_CTS <- dat_clean$CTS < 0 | dat_clean$CTS > 10
n_bad_cts <- sum(dat_clean$invalid_CTS, na.rm = TRUE)
dat_clean$CTS[dat_clean$invalid_CTS == TRUE & !is.na(dat_clean$invalid_CTS)] <- NA
log_msg("Set ", n_bad_cts, " out-of-range CTS values to NA.")

#  CT must be 0 or 1.
dat_clean$invalid_CT <- dat_clean$CT != 0 & dat_clean$CT != 1 & !is.na(dat_clean$CT)
n_bad_ct <- sum(dat_clean$invalid_CT, na.rm = TRUE)
dat_clean$CT[dat_clean$invalid_CT == TRUE] <- NA
log_msg("Set ", n_bad_ct, " invalid CT values to NA.")

#  Sex must be 0 or 1.
dat_clean$invalid_sex <- dat_clean$sex != 0 & dat_clean$sex != 1 & !is.na(dat_clean$sex)
n_bad_sex <- sum(dat_clean$invalid_sex, na.rm = TRUE)
dat_clean$sex[dat_clean$invalid_sex == TRUE] <- NA
log_msg("Set ", n_bad_sex, " invalid sex values to NA.")

# STOP AND ANSWER QUESTION 5 IN QUALTRICS BEFORE CONTINUING.
# Which coding or missing-value issues did you find, and how were they handled?

# ---- Activity 5: keep one row per participant -----------------------
#  Rule: keep the latest record_wave for each participant_id.
n_rows_before    <- nrow(dat_clean)
n_unique_ids     <- length(unique(dat_clean$participant_id))
n_repeat_rows    <- n_rows_before - n_unique_ids

dat_sorted <- dat_clean %>%
  arrange(participant_id, desc(record_wave)) %>%
  group_by(participant_id) %>%
  mutate(row_number_within_participant = row_number()) %>%
  ungroup()

dat_one_row <- dat_sorted %>%
  filter(row_number_within_participant == 1) %>%
  select(-row_number_within_participant)

log_msg(
  "Repeated-record handling: ", n_rows_before, " rows -> ", nrow(dat_one_row),
  " participants (", n_repeat_rows, " repeat rows dropped, latest record_wave kept)."
)

# STOP AND ANSWER QUESTION 6 IN QUALTRICS BEFORE CONTINUING.
# Explain whether duplicate/repeated records were present and how the script handled them.


# ---- Activity 6: apply coriander exposure exclusion -----------------
n_never_eaten <- sum(dat_one_row$ever_eaten_coriander == 0, na.rm = TRUE)
n_eaten_NA    <- sum(is.na(dat_one_row$ever_eaten_coriander))

dat_candidate <- dat_one_row %>%
  filter(ever_eaten_coriander == 1)

log_msg(
  "Exposure exclusion: dropped ", n_never_eaten,
  " never-eaten participants and ", n_eaten_NA, " with missing exposure status. ",
  "Candidate sample N = ", nrow(dat_candidate), "."
)

# STOP AND ANSWER QUESTION 7 IN QUALTRICS BEFORE CONTINUING.
# Report the candidate sample size and whether it meets the N > 100 criterion.

# QC note:
# Participants who have never eaten coriander will be excluded, even if CT or
# CTS values are present, because the taste phenotype is not interpretable for
# them. These inconsistencies should be addressed and reported.

# ---- Activity 7: implausible values + CT/CTS consistency flags -------
dat_candidate <- dat_candidate %>%
  mutate(
    flag_invalid_CT             = !is.na(CT)  & !(CT  %in% c(0, 1)),
    flag_invalid_CTS            = !is.na(CTS) & (CTS < 0 | CTS > 10),
    flag_invalid_sex            = !is.na(sex) & !(sex %in% c(0, 1)),
    flag_implausible_age        = !is.na(age) & (age < 0 | age > 120),
    flag_implausible_birth_year = !is.na(birth_year) & (birth_year < 1900 | birth_year > 2026),
    flag_CT0_high_CTS           = CT == 0 & !is.na(CTS) & CTS >= 6,
    flag_CT1_CTS0               = CT == 1 & !is.na(CTS) & CTS == 0,
    flag_inconsistent_CT_CTS    = (flag_CT0_high_CTS %in% TRUE) | (flag_CT1_CTS0 %in% TRUE),
    flag_problem_CT  = flag_invalid_CT  | flag_invalid_sex | flag_implausible_age |
      flag_implausible_birth_year | flag_inconsistent_CT_CTS,
    flag_problem_CTS = flag_invalid_CTS | flag_invalid_sex | flag_implausible_age |
      flag_implausible_birth_year | flag_inconsistent_CT_CTS
  )

log_msg(
  "Problem flags set on candidate sample: ",
  "flag_problem_CT = ",  sum(dat_candidate$flag_problem_CT,  na.rm = TRUE), ", ",
  "flag_problem_CTS = ", sum(dat_candidate$flag_problem_CTS, na.rm = TRUE), "."
)

# STOP AND ANSWER QUESTION 8 IN QUALTRICS BEFORE CONTINUING.
# Which suspicious values or inconsistencies remain, and what would you do with them?

############################################################
# Export dataset and log
# Note: This is required ONCE before running Phenotypic-CORE-worksheet.R
############################################################

# ---- Remove implausible age values ----
implausible_age <- dat_candidate %>% filter(!is.na(age) & (age < 0 | age > 120)) # Flag values

# Remove implausible values
dat_core <- dat_candidate %>%
  filter(!flag_implausible_age)

# ---- Write the cleaned dataset for the 20-min student script ---------
out_csv <- "coriander_gwas_pheno_clean_for_20min_v1.1.csv"
write.csv(dat_core, out_csv, row.names = FALSE)
log_msg("Wrote cleaned dataset: ", out_csv, " (", nrow(dat_core), " rows, ",
        ncol(dat_core), " columns).")

# ---- Write a small log file so we can audit what was done ------------
writeLines(log_lines, "pheno_qc_prep_log.txt")
message("\nDone. See pheno_qc_prep_log.txt for the full prep log.")


############################################################
# Additional stretch activities
# Use if there is extra time or as follow-up work
############################################################

# ---- Activity 8: Technical covariate dummy coding -------

# GWAS software often needs categorical covariates, such as array, batch and
# birth cohort, converted into numeric dummy variables.

# model.matrix() is a common base R function for this.

# Check the categories first.
table(dat_candidate$genotyping_array, useNA = "ifany")
table(dat_candidate$birth_cohort, useNA = "ifany")

if ("batch" %in% names(dat_candidate)) {
  table(dat_candidate$batch, useNA = "ifany")
}

# Create dummy variables for array, batch and birth cohort.
# One category is treated as the reference category automatically.
technical_dummies <- model.matrix(
  ~ genotyping_array + batch + birth_cohort,
  data = dat_candidate
)

# Inspect the dummy-coded covariates.
head(technical_dummies)

# Note 1:
# Dummy variables allow categorical covariates to be included in GWAS models.
# Always check your GWAS software documentation before using the file.

# Note 2:
# (Intercept) is automatically added by R/model.matrix() as part of the
# regression model structure. It is a standard modelling feature rather than
# a GWAS-specific covariate, and analysts usually do not need to interpret it.

# STOP AND ANSWER QUESTION 9 IN QUALTRICS BEFORE CONTINUING.
# Which categorical technical covariates were dummy coded, and why might GWAS software need dummy variables?


# ---- Activity 9: IQR outlier rule -------

# Outliers are unusual values, but they are not automatically errors.
# Here we use the IQR rule to flag CTS values for review.

# Calculate quartiles and IQR.
cts_q1 <- quantile(dat_candidate$CTS, 0.25, na.rm = TRUE)
cts_q3 <- quantile(dat_candidate$CTS, 0.75, na.rm = TRUE)
cts_iqr <- cts_q3 - cts_q1

# Calculate outlier thresholds.
cts_lower <- cts_q1 - 1.5 * cts_iqr
cts_upper <- cts_q3 + 1.5 * cts_iqr

cts_lower
cts_upper

# Flag CTS values outside the thresholds.
dat_candidate$CTS_iqr_outlier <- dat_candidate$CTS < cts_lower |
  dat_candidate$CTS > cts_upper

# Treat missing CTS as not flagged in this simple example.
dat_candidate$CTS_iqr_outlier[is.na(dat_candidate$CTS_iqr_outlier)] <- FALSE

# Count and inspect outliers.
table(dat_candidate$CTS_iqr_outlier, useNA = "ifany")

dat_candidate[
  dat_candidate$CTS_iqr_outlier == TRUE,
  c("participant_id", "CT", "CTS", "sex", "age")
]

# Note: IQR outliers should be reviewed. They should not be deleted automatically.

# STOP AND ANSWER QUESTION 10 IN QUALTRICS BEFORE CONTINUING.
# What did the IQR rule flag for CTS, and should those values be automatically removed?

# ---- Activity 10: Missingness patterns by age and sex -------

# This checks whether CTS missingness is related to observed variables.
# If missingness differs by sex or age, complete-case analysis may change the
# composition of the analysis sample.

# Create a TRUE/FALSE missingness indicator.
dat_candidate$missing_CTS <- is.na(dat_candidate$CTS)

# Count missing and observed CTS values.
table(dat_candidate$missing_CTS, useNA = "ifany")

# Check missingness by sex.
missing_by_sex <- table(dat_candidate$sex, dat_candidate$missing_CTS, useNA = "ifany")
missing_by_sex

# Convert to percentages within each sex group.
prop.table(missing_by_sex, margin = 1) * 100

# Compare age between observed and missing CTS groups.
aggregate(
  age ~ missing_CTS,
  data = dat_candidate,
  FUN = mean,
  na.rm = TRUE
)

# Plot age by CTS missingness.
ggplot(dat_candidate, aes(x = factor(missing_CTS), y = age)) +
  geom_boxplot(na.rm = TRUE) +
  labs(
    x = "CTS missing?",
    y = "Age",
    title = "Age by CTS missingness",
    subtitle = "FALSE = observed, TRUE = missing"
  ) +
  theme_minimal(base_size = 13)

# Note:
# This does not prove why data are missing. It only checks whether missingness is
# associated with observed variables.

# STOP AND ANSWER QUESTION 11 IN QUALTRICS BEFORE CONTINUING.
# Does CTS missingness appear related to sex or age? Why does this matter for GWAS?

# ---- Activity 11: Save cleaned files for downstream GWAS software -------

# Pre-activity setup

# CT and CTS are related phenotypes, but they are not the same.
# CT is binary:
#   0 = participant does not report coriander/cilantro tasting soapy
#   1 = participant reports coriander/cilantro tasting soapy
# CTS is continuous:
#   0  = not at all soapy
#   10 = completely soapy
# Because these two variables measure related information, we expect them to
# broadly agree. For example:
#   - CT = 0 should usually go with low CTS values.
#   - CT = 1 should usually go with CTS values above 0.
# CT = 0 with very high CTS, or CT = 1 with CTS = 0, may need review.

# Select variables needed for two hypothetical GWAS: CT or CTS as outcomes
required_for_CT_gwas <- c(ct_var, sex_var, age_var, pc_vars, array_var, batch_var, cohort_var)
required_for_CT_gwas <- intersect(required_for_CT_gwas, names(dat_candidate))

required_for_CTS_gwas <- c(cts_var, sex_var, age_var, pc_vars, array_var, batch_var, cohort_var)
required_for_CTS_gwas <- intersect(required_for_CTS_gwas, names(dat_candidate))

# Create complete-case dataset
dat_CT_complete <- dat_candidate %>%
  filter(if_all(all_of(required_for_CT_gwas), ~ !is.na(.x)))
dat_CTS_complete <- dat_candidate %>%
  filter(if_all(all_of(required_for_CTS_gwas), ~ !is.na(.x)))

# Check how row numbers differ
nrow(dat_candidate)
nrow(dat_CT_complete)
nrow(dat_CTS_complete)

# Create final CT analysis sample.
dat_final_CT <- dat_candidate %>%
  filter(
    ever_eaten_coriander == 1,
    !flag_problem_CT,
    if_all(all_of(required_for_CT_gwas), ~ !is.na(.x)),
    CT %in% c(0, 1)
  )

# Create final CTS analysis sample.
dat_final_CTS <- dat_candidate %>%
  filter(
    ever_eaten_coriander == 1,
    !flag_problem_CTS,
    if_all(all_of(required_for_CTS_gwas), ~ !is.na(.x)),
    CTS >= 0,
    CTS <= 10
  )

# Final sample sizes.
nrow(dat_final_CT)
nrow(dat_final_CTS)

# Save cleaned phenotype files.
write.csv(dat_final_CT, "final_CT_phenotype_for_GWAS.csv", row.names = FALSE)
write.csv(dat_final_CTS, "final_CTS_phenotype_for_GWAS.csv", row.names = FALSE)

# Create a simple CT covariate file.
# Keep participant_id so the file can be matched to genotype data.
covar_CT <- dat_final_CT[
  c(
    "participant_id",
    "sex",
    "age",
    "PC1", "PC2", "PC3", "PC4", "PC5",
    "PC6", "PC7", "PC8", "PC9", "PC10",
    "genotyping_array",
    "batch",
    "birth_cohort"
  )
]

# Convert categorical technical covariates to dummy variables.
technical_dummies <- model.matrix(
  ~ genotyping_array + batch + birth_cohort,
  data = covar_CT
)

# Keep numeric covariates.
numeric_covariates <- covar_CT[
  c(
    "participant_id",
    "sex",
    "age",
    "PC1", "PC2", "PC3", "PC4", "PC5",
    "PC6", "PC7", "PC8", "PC9", "PC10"
  )
]

# Combine numeric covariates and dummy variables.
covar_CT_out <- data.frame(
  numeric_covariates,
  technical_dummies,
  check.names = FALSE
)

head(covar_CT_out)

# Save covariate file.
write.csv(covar_CT_out, "CT_covariate_matrix.csv", row.names = FALSE)

# Note:
# Real GWAS software differs in how it expects phenotype and covariate files to
# be formatted. Always check the documentation before using these files.

# STOP AND ANSWER QUESTION 12 IN QUALTRICS BEFORE CONTINUING.
# Which files were saved, and why is participant_id kept in the covariate file?


# ---- Activity 12: Missing-data patterns with mice -------

# This is optional. It uses the mice package to show patterns of missingness.
# Each row in the output shows a different missing-data pattern.
#
# This can help us see whether missingness tends to occur in one variable only,
# or whether the same participants are missing data across several variables.

# Load mice only here because this is a stretch activity.
# If you do not have mice installed, run:
# install.packages("mice")
library(mice)

# Select a small number of key variables.
# md.pattern() is easier to read when we include only the most important
# phenotype and covariate variables.
#
# Note:
# We use all_of(batch_var) because batch_var stores the name of the batch column.
# In this dataset, batch_var is usually "batch".
missing_pattern_vars <- c(
  "CT",
  "CTS",
  "sex",
  "age",
  "genotyping_array",
  batch_var,
  "birth_cohort"
)

# Keep only variables that actually exist in dat_candidate.
missing_pattern_vars <- intersect(missing_pattern_vars, names(dat_candidate))

# Show missing-data patterns.
# In the output:
#   1 usually means the variable is observed.
#   0 usually means the variable is missing.
#
# The numbers at the left show how many participants have each pattern.
# The bottom row shows how many values are missing in each variable.
md.pattern(
  dat_candidate %>%
    select(all_of(missing_pattern_vars))
)

# Interpretation:
# A large complete-case row means many participants have all selected variables.
# Rows with several 0s show participants missing multiple variables.
# Variables with larger numbers in the bottom row contribute more to sample loss.

# STOP AND ANSWER QUESTION 13 IN QUALTRICS BEFORE CONTINUING.
# What does md.pattern() show, and how can this help with missing-data QC?


############################################################
# End of stretch activities
############################################################