Testing UserGWAS/GenomicSEM for Value with Three Moderately Genetically Correlated Chemokines

1 Analysis Summary

The analysis integrates Genomic Structural Equation Modeling (GenomicSEM) with multivariate summary statistics of three chemokines—CXCL3, CCL4, and CCL27—identified as moderately genetically correlated through linkage disequilibrium score regression (LDSC). By leveraging multivariate methods, the study highlights the advantages of capturing shared genetic architecture across traits while reducing noise compared to traditional univariate approaches.

Key findings include the identification of two genomic loci, on chromosomes 2 and 7, containing variants (e.g., rs1260326 near GCKR and rs7779873 near CD36) with significant associations and potential implications in metabolic pathways. Functional annotation using tools such as CADD/PHRED and FUMA underscores the biological relevance of these loci.

The value of GenomicSEM lies in its ability to provide a robust framework for testing the multivariate genetic architecture of traits, enhancing statistical power and uncovering pleiotropic effects that may remain undetected in pairwise or univariate analyses. By integrating functional annotations and mapping tools, this approach strengthens the interpretation of genetic findings and their translational potential in complex traits, such as inflammatory and metabolic disorders.

1.1 Bioinformatics Processing: Fetching rsID and Pvals

Show

#!/bin/bash

# Directories and variables
TAR_DIR="/Users/charleenadams/ukbppp/chemokine_list"
MAP_DIR="/Users/charleenadams/temp_BI/chemokine_rgs_olink/maps"
OUTPUT_BASE_DIR="/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list"
MERGED_DIR="$OUTPUT_BASE_DIR/merged_files"
LOG_FILE="$OUTPUT_BASE_DIR/processing_log.txt"

# Create required directories
mkdir -p "$OUTPUT_BASE_DIR"
mkdir -p "$MERGED_DIR"
echo "Processing started on $(date)" > "$LOG_FILE"

# Function to process each tar file
process_tar_file() {
    local tar_file="$1"
    local base_name=$(basename "$tar_file" .tar)
    local untarred_dir="$OUTPUT_BASE_DIR/$base_name"
    local output_dir="$OUTPUT_BASE_DIR/$base_name"

    echo "Processing tar file: $tar_file" | tee -a "$LOG_FILE"

    # Untar the file into the correct directory
    mkdir -p "$untarred_dir"
    tar -xf "$tar_file" -C "$untarred_dir" || {
        echo "Failed to untar: $tar_file" | tee -a "$LOG_FILE"
        return
    }

    # Process chromosome files within the untarred directory
    mkdir -p "$output_dir"
    local processed_files=()
    find "$untarred_dir" -type f -name "*.gz" | while read -r chr_file; do
        # Skip files that are already processed
        if [[ "$chr_file" == *"_rsid_added.gz" || "$chr_file" == *"_with_pval.gz" ]]; then
            echo "Skipping already processed file: $chr_file" | tee -a "$LOG_FILE"
            continue
        fi

        local chrom=$(basename "$chr_file" | sed -E 's/.*chr([0-9XY]+).*/\1/')
        local map_file="$MAP_DIR/olink_rsid_map_mac5_info03_b0_7_chr${chrom}_patched_v2.tsv.gz"
        local output_file="${chr_file%.gz}_rsid_added.gz"

        if [[ -n "$chrom" && -f "$map_file" ]]; then
            echo "Adding RSID to $chr_file (Chromosome $chrom)" | tee -a "$LOG_FILE"
            gunzip -c "$chr_file" | \
            awk 'BEGIN {OFS="\t"}
                 NR==FNR {map[$1]=$4; genpos37[$1]=$5; genpos38[$1]=$6; next}
                 ($3 in map) {print $0, map[$3], genpos37[$3], genpos38[$3]}
                 !($3 in map) {print $0, ".", ".", "."}' <(gunzip -c "$map_file") - | \
            gzip > "$output_file"

            if [[ -s "$output_file" ]]; then
                echo "RSID added to $output_file" | tee -a "$LOG_FILE"
                processed_files+=("$output_file")

                # Add P-values if not already added
                local pval_file="${output_file%.gz}_with_pval.gz"
                if [[ ! -f "$pval_file" ]]; then
                    gunzip -c "$output_file" | \
                    awk 'BEGIN {OFS="\t"} NR==1 {print $0, "P"} NR>1 {pval = (10 ^ -$13); print $0, pval}' | \
                    gzip > "$pval_file"
                    mv "$pval_file" "$output_dir/$(basename "$pval_file")"
                else
                    echo "P-values already added for $output_file" | tee -a "$LOG_FILE"
                fi
            else
                echo "Warning: RSID addition failed for $chr_file" | tee -a "$LOG_FILE"
            fi
        else
            echo "RSID map file for chromosome $chrom not found." | tee -a "$LOG_FILE"
        fi
    done

    # Merge chromosome files in the output directory
    if [[ ${#processed_files[@]} -gt 0 ]]; then
        echo "Merging chromosome files in $output_dir" | tee -a "$LOG_FILE"
        local merged_file="$MERGED_DIR/${base_name}_merged_rsid_added.gz"

        # Extract header from the first file
        gunzip -c "${processed_files[0]}" | head -n 1 > temp_header.txt

        # Process each file, excluding the header, and append to a temporary content file
        for file in "${processed_files[@]}"; do
            gunzip -c "$file" | tail -n +2 >> temp_content.txt
        done

        # Concatenate header and content into the final merged file
        cat temp_header.txt temp_content.txt | gzip -c > "$merged_file"

        # Clean up temporary files
        rm -f temp_header.txt temp_content.txt

        if [[ -s "$merged_file" ]]; then
            echo "Merged file created: $merged_file" | tee -a "$LOG_FILE"
        else
            echo "Warning: Merged file is empty for $output_dir" | tee -a "$LOG_FILE"
        fi
    else
        echo "No files to merge in $output_dir" | tee -a "$LOG_FILE"
    fi
}

# Export the function and variables for parallel processing
export -f process_tar_file
export MAP_DIR OUTPUT_BASE_DIR LOG_FILE MERGED_DIR

# Find all tar files and process them in parallel
find "$TAR_DIR" -type f -name "*.tar" | parallel -j15 process_tar_file

echo "Processing completed on $(date)" >> "$LOG_FILE"
echo "All results saved in $OUTPUT_BASE_DIR and merged files in $MERGED_DIR"

# tmux new -s my_session # didn't use tmux; keeping command to trigger memory
# nohup ./process_40_chemokines.sh > process_40_chemokines_log 2>&1 &
# pstree 91532
# ps -p 91532 -o stat,etime
# watch -n 1 ls -lt /Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list
# find /Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list -type f -name "*with_pval.gz" | wc -l

# remove the extra files
main_dir="/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list"

clean_directory() {
    local dir="$1"
    find "$dir" -type f -not -name "*_rsid_added_with_pval.gz" -not -name "*.txt" -delete
    echo "Processed: $dir"
}

find "$main_dir" -type d | while read -r dir; do
    clean_directory "$dir"
done

1.2 Munge by Chromsome

Show

library(GenomicSEM)
library(parallel)

# Set directories
base_dir <- "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list"
hm3_snplist <- "/Users/charleenadams/temp_BI/chemokine_rgs_olink/w_hm3.snplist"
output_directory <- "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40"

# Ensure the output directory exists
if (!dir.exists(output_directory)) {
  dir.create(output_directory, recursive = TRUE)
}

# Define column names as a named list
column_names <- list(
  SNP = "rsid",        # The column for SNP identifiers
  A1 = "ALLELE1",      # The column for the first allele
  A2 = "ALLELE0",      # The column for the second allele
  effect = "BETA",     # The column for the effect size
  INFO = "INFO",       # The column for the INFO score
  P = "P",             # The column for the P-value
  N = "N",             # The column for the sample size
  MAF = "A1FREQ",      # The column for the minor allele frequency
  Z = "CHISQ"          # The column for the Z-score or chi-squared value
)

# Function to extract `N` from the first row of a file
extract_N <- function(file, column) {
  header <- read.table(gzfile(file), nrows = 1, header = TRUE)
  if (column %in% colnames(header)) {
    N <- unique(header[[column]])
    if (length(N) == 1) {
      return(as.numeric(N))
    } else {
      stop(paste("N column has multiple values in file:", file))
    }
  } else {
    stop(paste("Column", column, "not found in file:", file))
  }
}

# Function to process a single subdirectory
process_subdirectory <- function(subdir) {
  tryCatch({
    # Get subdirectory name
    subdir_name <- basename(subdir)
    
    # Create output subdirectory
    subdir_output <- file.path(output_directory, subdir_name)
    if (!dir.exists(subdir_output)) {
      dir.create(subdir_output, recursive = TRUE)
    }
    
    # Get the list of chromosome-specific files for chromosomes 1–22
    chr_files <- list.files(subdir, full.names = TRUE, pattern = "chr[1-9]|chr1[0-9]|chr2[0-2]_.*rsid_added_with_pval\\.gz$")
    
    if (length(chr_files) == 0) {
      cat(paste("No files found in subdirectory:", subdir, "\n"))
      return(NULL)
    }
    
    # Extract `N` for each file
    N_values <- sapply(chr_files, extract_N, column = column_names$N)
    
    # Define trait names based on file names
    trait_names <- gsub("_rsid_added_with_pval\\.gz$", "", basename(chr_files))
    
    # Run the munge function with internal parallelization
    munge(
      files = chr_files,
      hm3 = hm3_snplist,
      trait.names = trait_names,
      N = N_values,
      info.filter = 0.9,       # Filter by INFO score
      maf.filter = 0.01,       # Filter by MAF
      column.names = column_names,
      parallel = TRUE,         # Enable internal parallelization
      cores = detectCores() - 1,  # Use all available cores minus one
      overwrite = TRUE         # Allow overwriting of existing files
    )
    
    cat(paste("Processing completed for subdirectory:", subdir, "\n"))
  }, error = function(e) {
    cat(paste("Error processing subdirectory:", subdir, "\nError message:", e$message, "\n"))
  })
}

# Get all subdirectories
subdirs <- list.dirs(base_dir, full.names = TRUE, recursive = FALSE)

# Sequentially process subdirectories
for (subdir in subdirs) {
  process_subdirectory(subdir)
}

1.3 Move to munged_40

Show

ls -lt /Users/charleenadams/ukbppp/*.sumstats.gz | grep "Dec  8" | wc -l
# Without notes:
directory="/Users/charleenadams/ukbppp"
target_date="Dec  8"
output_dir="/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40"
mkdir -p "$output_dir"

while IFS= read -r file; do
  cp "$file" "$output_dir"
done < <(ls -lt "$directory"/*.sumstats.gz | grep "$target_date" | awk '{print $NF}')
ls -lt /Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/*.sumstats.gz | grep "Dec  8" | wc -l

1.4 Merge the Munged Files for Each Protein

nohup Rscript merged_munged_40.R > merged_munged_40.log 2>&1 &

Show

library(parallel)

# Define directories
base_dir <- "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40"
results_dir <- file.path(base_dir, "merged_munged_protein_results")

# Ensure the results directory exists
if (!dir.exists(results_dir)) {
  dir.create(results_dir, recursive = TRUE)
}

# Function to merge chromosome files for a protein
merge_chromosome_files <- function(protein_name, chr_files, output_dir) {
  tryCatch({
    # Define output file path
    merged_file <- file.path(output_dir, paste0(protein_name, "_merged.sumstats.gz"))
    
    # Merge chromosome files
    all_data <- do.call(rbind, lapply(chr_files, function(f) {
      read.table(f, header = TRUE, sep = "\t", stringsAsFactors = FALSE)
    }))
    
    # Save the merged file
    write.table(
      all_data,
      file = gzfile(merged_file),
      quote = FALSE,
      sep = "\t",
      row.names = FALSE,
      col.names = TRUE
    )
    cat(paste("Merged protein:", protein_name, "to", merged_file, "\n"))
    return(merged_file)
  }, error = function(e) {
    cat(paste("Error merging protein:", protein_name, "\nError message:", e$message, "\n"))
    return(NULL)
  })
}

# Get list of all chromosome files
chr_files <- list.files(base_dir, pattern = "discovery_chr.*\\.sumstats\\.gz$", full.names = TRUE)

# Extract unique protein identifiers
proteins <- unique(gsub("discovery_chr[0-9]+_", "", gsub("\\.sumstats\\.gz$", "", basename(chr_files))))

# Function to process a single protein
process_protein <- function(protein) {
  protein_files <- grep(protein, chr_files, value = TRUE)
  merge_chromosome_files(protein, protein_files, output_dir = results_dir)
}

# Run merging in parallel for 12 proteins
num_cores <- detectCores() - 1  # Use all but one core
cl <- makeCluster(num_cores)
clusterExport(cl, c("chr_files", "proteins", "merge_chromosome_files", "results_dir"))
clusterEvalQ(cl, library(parallel))
merged_files <- parLapply(cl, proteins, process_protein)
stopCluster(cl)

# Print summary
cat("Merging completed for all proteins.\n")

# nohup Rscript merged_munged_40.R > merged_munged_40.log 2>&1 &

1.5 Add P Back for Sumstats

nohup Rscript add_p_sumstats.R > add_p_sumstats.log 2>&1 &

Show

library(data.table)

# Define the function to add a P column
add_p_column <- function(file, output_dir) {
  tryCatch({
    data <- fread(file)
    if (!"P" %in% colnames(data)) {
      if ("Z" %in% colnames(data)) {
        # Calculate P-values from Z-scores
        data[, P := 2 * pnorm(-abs(Z))]
        
        # Create output file path
        output_file <- file.path(output_dir, basename(file))
        
        # Save the updated data to the new file
        fwrite(data, output_file, sep = "\t")
        cat("Added P column to:", output_file, "\n")
      } else {
        stop("Z column missing. Cannot calculate P-values in file:", file)
      }
    } else {
      cat("P column already exists in:", file, "\n")
    }
  }, error = function(e) {
    cat("Error processing file:", file, "\nError message:", e$message, "\n")
  })
}

# Set input and output directories
input_dir <- "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results"
output_dir <- file.path(input_dir, "pvals_fix")

# Create the output directory if it doesn't exist
if (!dir.exists(output_dir)) {
  dir.create(output_dir)
}

# List all files matching the pattern in the input directory
files <- list.files(input_dir, pattern = "\\.sumstats\\.gz$", full.names = TRUE)

# Apply the function to all files
lapply(files, add_p_column, output_dir = output_dir)

1.6 Sumstats of Munged Data

nohup Rscript sumstats_on_sumstats.R > sumstats_on_sumstats.log 2>&1 &

Show

library(data.table)
library(GenomicSEM)

# Step 1: Define paths and variables
input_dir <- "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix"  # Input directory containing sumstats files
output_path <- "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/summary_stats.txt"      # Path to save the final summary stats
ref <- "/Users/charleenadams/temp_BI/chemokine_rgs_olink/reference.1000G.maf.0.005.txt"      # Reference file for SNP variance

# List all sumstats files
traits <- list.files(
  path = input_dir,
  pattern = "\\.sumstats\\.gz$",  # Match gzipped sumstats files
  full.names = TRUE
)

# Optional: Name your traits based on filenames without extensions
trait.names <- gsub("\\_.sumstats\\.gz$", "", basename(traits))

# Step 2: Extract sample sizes (N) from each trait file
# Assuming the column name for sample size is "N"
get_sample_size <- function(file) {
  data <- fread(file, select = "N", nrows = 1)  # Only read the first row to extract N
  return(data$N[1])
}

N <- sapply(traits, get_sample_size)
cat("Extracted sample sizes:\n")

# Step 3: Prepare the summary statistics with sumstats()
# Define arguments for sumstats()
se.logit <- rep(FALSE, length(traits))  # Continuous traits; SEs are not logistic
OLS <- rep(TRUE, length(traits))       # Ordinary least squares for continuous traits
linprob <- rep(FALSE, length(traits))  # Traits are not dichotomous; linear probability does not apply

# Run sumstats() to align and scale summary statistics
aligned_sumstats <- sumstats(
  files = traits,
  ref = ref,
  trait.names = trait.names,
  se.logit = se.logit,
  OLS = OLS,
  linprob = linprob,
  N = N,
  betas = NULL,
  #maf.filter = 0.01,
  keep.indel = FALSE,
  parallel = TRUE,
  cores = 15  # Adjust for available cores
)

# Save the aligned summary statistics
fwrite(aligned_sumstats, file = output_path, sep = "\t")
cat("Aligned summary statistics saved to:", output_path, "\n")

# Optional: Verify the saved file by re-loading it
summary_stats_read <- fread(output_path)

1.7 MV LDSC on 40 Chemokines

nohup Rscript ldsc40.R > ldsc40.log 2>&1 &

Show

# Load required library
library(GenomicSEM)
library(data.table)
library(parallel)
library(qqman)

# Define paths and variables
traits <- list.files("/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix", 
                     pattern = "\\.sumstats\\.gz$", full.names = TRUE)

# Set sample and population prevalences to NA to bypass liability scaling
sample.prev <- rep(NA, length(traits))  # Treat traits as continuous
population.prev <- rep(NA, length(traits))

ld <- "/Users/charleenadams/temp_BI/mvGWAS_chemokines/eur_w_ld_chr/"
wld <- "/Users/charleenadams/temp_BI/mvGWAS_chemokines/eur_w_ld_chr/"

# Optional: Name your traits based on filenames without extensions
trait.names <- gsub("\\.sumstats\\.gz$", "", basename(traits))
head(trait.names)

# Redirect console output to a log file
sink("/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/ldsc_40_rgs.log")

# Run multivariable LDSC
LDSCoutput <- ldsc(traits = traits,
                   sample.prev = sample.prev,
                   population.prev = population.prev,
                   ld = ld,
                   wld = wld,
                   trait.names = trait.names,
                   stand = TRUE)

# Stop capturing the console output
sink()

# Save the results
output_path <- "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/ldsc_results_rgs"
saveRDS(LDSCoutput, file = file.path(output_path, "LDSCoutput40.rds"))

# Print completion message
cat("Multivariable LDSC analysis complete and save in:\n", output_path, "\n")
</div>

1.8 Selected Three Moderately Genetically Correlated Chemokines to Format with sumstats: CXCL3 CCL4 CCL27

Pair	Chemokine 1	Chemokine 2	Genetic Correlation	Standard Error
CCL4_CXCL3	CCL4	CXCL3	0.3367	0.1022
CCL27_CXCL3	CCL27	CXCL3	0.4379	0.1042
CCL27_CCL4	CCL27	CCL4	0.2758	0.1001

• CXCL3:P19876:OID20788:v1:Inflammation_merged.sumstats.gz

• CCL4:P13236:OID20695:v1:Inflammation_merged.sumstats.gz

• CCL27:Q9Y4X3:OID20121:v1:Cardiometabolic_merged.sumstats.gz

Show

library(data.table)
library(GenomicSEM)

# Step 1: Define paths and variables
input_dir <- "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix"  # Input directory containing sumstats files
output_path <- "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/summary_stats_3_chemokines.txt"  # Path to save the final summary stats
ref <- "/Users/charleenadams/temp_BI/chemokine_rgs_olink/reference.1000G.maf.0.005.txt"  # Reference file for SNP variance

# Step 2: Specify the exact filenames of the desired chemokine files
selected_files <- c(
  "CXCL3:P19876:OID20788:v1:Inflammation_merged.sumstats.gz",
  "CCL4:P13236:OID20695:v1:Inflammation_merged.sumstats.gz",
  "CCL27:Q9Y4X3:OID20121:v1:Cardiometabolic_merged.sumstats.gz"
)
traits_filtered <- file.path(input_dir, selected_files)

# Check if all specified files exist
missing_files <- traits_filtered[!file.exists(traits_filtered)]
if (length(missing_files) > 0) {
  stop("The following files are missing:\n", paste(missing_files, collapse = "\n"))
}

# Optional: Name your traits based on filenames without extensions
trait.names <- gsub("\\.sumstats\\.gz$", "", basename(traits_filtered))

cat("Selected traits:\n", paste(trait.names, collapse = "\n"), "\n")

# Step 3: Extract sample sizes (N) from each trait file
get_sample_size <- function(file) {
  data <- fread(file, select = "N", nrows = 1)  # Only read the first row to extract N
  return(data$N[1])
}

N <- sapply(traits_filtered, get_sample_size)
cat("Extracted sample sizes:\n", paste(N, collapse = "\n"), "\n")

# Step 4: Prepare the summary statistics with sumstats()
# Define arguments for sumstats()
se.logit <- rep(FALSE, length(traits_filtered))  # Continuous traits; SEs are not logistic
OLS <- rep(TRUE, length(traits_filtered))       # Ordinary least squares for continuous traits
linprob <- rep(FALSE, length(traits_filtered))  # Traits are not dichotomous; linear probability does not apply

# Run sumstats() to align and scale summary statistics
aligned_sumstats <- sumstats(
  files = traits_filtered,
  ref = ref,
  trait.names = trait.names,
  se.logit = se.logit,
  OLS = OLS,
  linprob = linprob,
  N = N,
  betas = NULL,
  #maf.filter = 0.01,
  keep.indel = FALSE,
  parallel = TRUE,
  cores = 15  # Adjust for available cores
)

# Save the aligned summary statistics
fwrite(aligned_sumstats, file = output_path, sep = "\t")
cat("Aligned summary statistics saved to:", output_path, "\n")

# Optional: Verify the saved file by re-loading it
summary_stats_read <- fread(output_path)
print(head(summary_stats_read))

1.9 LDSC For Three Moderately Genetically Correlated Chemokines

Show

# Load required library
library(GenomicSEM)
library(data.table)
library(parallel)
library(qqman)

# Define paths and variables
input_dir <- "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix"
output_path <- "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/ldsc_results_rgs"
log_file <- file.path(output_path, "ldsc_3_chemokines.log")
selected_files <- c(
  "CXCL3:P19876:OID20788:v1:Inflammation_merged.sumstats.gz",
  "CCL4:P13236:OID20695:v1:Inflammation_merged.sumstats.gz",
  "CCL27:Q9Y4X3:OID20121:v1:Cardiometabolic_merged.sumstats.gz"
)
traits <- file.path(input_dir, selected_files)

# Check if files exist
missing_files <- traits[!file.exists(traits)]
if (length(missing_files) > 0) {
  stop("The following files are missing:\n", paste(missing_files, collapse = "\n"))
}

# Set sample and population prevalences to NA to bypass liability scaling
sample.prev <- rep(NA, length(traits))  # Treat traits as continuous
population.prev <- rep(NA, length(traits))

# Define LD and weight directories
ld <- "/Users/charleenadams/temp_BI/mvGWAS_chemokines/eur_w_ld_chr/"
wld <- "/Users/charleenadams/temp_BI/mvGWAS_chemokines/eur_w_ld_chr/"

# Optional: Name your traits based on filenames without extensions
trait.names <- gsub("\\.sumstats\\.gz$", "", basename(traits))

# Create output directory if it doesn't exist
if (!dir.exists(output_path)) {
  dir.create(output_path, recursive = TRUE)
}

# Redirect console output to a log file
sink(log_file)

# Run multivariable LDSC
LDSCoutput <- ldsc(traits = traits,
                   sample.prev = sample.prev,
                   population.prev = population.prev,
                   ld = ld,
                   wld = wld,
                   trait.names = trait.names,
                   stand = TRUE)

# Stop capturing the console output
sink()

# Save the results
saveRDS(LDSCoutput, file = file.path(output_path, "LDSCoutput_3_chemokines.rds"))

# Print completion message
cat("Multivariable LDSC analysis complete. Results saved in:\n", output_path, "\n")# Load required library

1.10 Parallel Analysis Based on Multivariate LDSC

1.10.1 Code for Parallel Analysis Based on Multivariate LDSC

Show

library(GenomicSEM)
library(knitr)
LDSCoutputfull <- readRDS("/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/ldsc_results_rgs/LDSCoutput_3_chemokines.rds")
#2. Genetic correlation matrix from LDSC output
S_Stand <- LDSCoutputfull$S_Stand
#3. Associated standardized multivariate sampling distribution matrix from LDSC output
V_Stand <- LDSCoutputfull$V_Stand

# Save the plot as a 600 DPI PNG
output_path <- "paLDSC_plot.png"
png(filename = output_path, width = 8, height = 6, units = "in", res = 600)

# Run the `paLDSC` function
paLDSC(S = S_Stand, V = V_Stand, r = 10)

## Running parallel analysis. Replication number:  1 
## Running parallel analysis. Replication number:  2 
## Running parallel analysis. Replication number:  3 
## Running parallel analysis. Replication number:  4 
## Running parallel analysis. Replication number:  5 
## Running parallel analysis. Replication number:  6 
## Running parallel analysis. Replication number:  7 
## Running parallel analysis. Replication number:  8 
## Running parallel analysis. Replication number:  9 
## Running parallel analysis. Replication number:  10

## ------------------------------------------------------------------------ 
## Parallel Analysis suggests extracting 1 components 
## ------------------------------------------------------------------------

# Close the graphics device
dev.off()

## quartz_off_screen 
##                 2

1.10.2 Parallel Analysis Based on Multivariate LDSC: Decision Figure

# Path to the saved plot
plot_path <- "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/ldsc_results_rgs/paLDSC_plot.png"

# Include the plot
include_graphics(plot_path)

1.11 Exploratory Factor Analysis

Show

# Load necessary libraries
library(psych)  # For EFA
library(GPArotation)  # For rotations

# Define the correlation matrix
# Extract the LDSC genetic correlation matrix
ldsc_results <- readRDS("/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/ldsc_results_rgs/LDSCoutput_3_chemokines.rds")
cor_matrix <- ldsc_results$S_Stand
write.csv(cor_matrix, "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/ldsc_results_rgs/EFA_cor_matrix3.csv")

cor_matrix <- cor.smooth(cor_matrix)
eigenvalues <- eigen(cor_matrix)$values
if (any(eigenvalues <= 0)) {
    cat("The smoothed matrix is still not positive definite. Smallest eigenvalue:", min(eigenvalues), "\n")
} else {
    cat("Matrix is now positive definite.\n")
}

## Matrix is now positive definite.

# Run EFA to determine the number of latent factors
efa_results1 <- fa(
  r = cor_matrix,       # Genetic correlation matrix
  nfactors = 1,         # Test up to 1 factors 
  rotate = "oblimin",   # Allow factors to correlate
  fm = "ml"             # Use maximum likelihood estimation
)

efa_results1

## Factor Analysis using method =  ml
## Call: fa(r = cor_matrix, nfactors = 1, rotate = "oblimin", fm = "ml")
## Standardized loadings (pattern matrix) based upon correlation matrix
##                                                  ML1   h2   u2 com
## CXCL3:P19876:OID20788:v1:Inflammation_merged    0.73 0.53 0.47   1
## CCL4:P13236:OID20695:v1:Inflammation_merged     0.46 0.21 0.79   1
## CCL27:Q9Y4X3:OID20121:v1:Cardiometabolic_merged 0.60 0.36 0.64   1
## 
##                 ML1
## SS loadings    1.11
## Proportion Var 0.37
## 
## Mean item complexity =  1
## Test of the hypothesis that 1 factor is sufficient.
## 
## df null model =  3  with the objective function =  0.36
## df of  the model are 0  and the objective function was  0 
## 
## The root mean square of the residuals (RMSR) is  0 
## The df corrected root mean square of the residuals is  NA 
## 
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy             
##                                                    ML1
## Correlation of (regression) scores with factors   0.81
## Multiple R square of scores with factors          0.66
## Minimum correlation of possible factor scores     0.33

# Run EFA to determine the number of latent factors
#efa_results2 <- fa(
#  r = cor_matrix,       # Genetic correlation matrix
#  nfactors = 2,         # Test up to 1 factors 
#  rotate = "oblimin",   # Allow factors to correlate
#  fm = "ml"             # Use maximum likelihood estimation
#)

1.12 One-Factor Multivariate GWAS (userGWAS)

1.12.1 Background Definitions

Structural Equation Modeling (SEM) is a statistical technique used to model and analyze complex relationships between observed (measured) and unobserved (latent) variables. It combines elements of multiple regression, factor analysis, and path analysis into a single, flexible framework.

Likewise, Genomic Structural Equation Modeling (GenomicSEM) is a statistical framework designed for SEM with GWAS summary statistics. It adapts SEM principles to work with genetic covariance matrices, enabling the study of complex relationships among traits, shared genetic architectures, and causal pathways. It uses two matrices:

1. S Matrix (Genetic Covariance Matrix):
   • The S matrix represents the genetic covariance among traits. It is derived from LDSC and provides unstandardized estimates of genetic covariances based on GWAS summary statistics.
   • For example, since we are analyzing three traits (e.g., CXCL3, CCL4, and CCL27), the S matrix would include pairwise genetic covariances between these traits.
2. V Matrix (Sampling Covariance Matrix):
   • The V matrix is the corresponding sampling covariance matrix for the genetic covariance estimates in the S matrix. It reflects the uncertainty (variance) and interdependence (covariance) of the estimates in the S matrix.
   • The V matrix is crucial for properly accounting for standard errors and dependencies in the modeling process.

---

• Modeling:

  • The S matrix is used to model the relationships between traits.
  • The V matrix ensures accurate statistical inference by incorporating uncertainties in the genetic covariance estimates.

---

1. General SEM Framework:

   • lavaan is a tool for specifying and estimating structural equation models. In GenomicSEM, the structural equation modeling principles are implemented in lavaan syntax. That is, the model is specified in lavaan.

2. Parameter Estimation:

   • lavaan relies on covariance matrices for parameter estimation in SEM.

---

nohup Rscript usergwas_3chemokines.R > usergwas_3chemokines.log 2>&1 &

Show

# Load necessary libraries
library(data.table)
library(GenomicSEM)
library(TwoSampleMR)
library(dplyr)

# Start timer to measure runtime
start_time <- Sys.time()

# Define paths and files
summary_stats_path <- "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/summary_stats_3_chemokines.txt"
ldsc_results_path <- "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/ldsc_results_rgs/LDSCoutput_3_chemokines.rds"

# Load summary statistics
summary_stats <- fread(summary_stats_path)

# Replace colons with underscores in column names for consistency
colnames(summary_stats) <- gsub(":", "_", colnames(summary_stats))

# Load LDSC results
ldsc_results <- readRDS(ldsc_results_path)

# Replace colons with underscores in the names of ldsc_results
if (!is.null(attr(ldsc_results$S, "dimnames"))) {
  attr(ldsc_results$S, "dimnames")[[2]] <- gsub(":", "_", attr(ldsc_results$S, "dimnames")[[2]])
}
if (!is.null(attr(ldsc_results$S_Stand, "dimnames"))) {
  attr(ldsc_results$S_Stand, "dimnames")[[2]] <- gsub(":", "_", attr(ldsc_results$S_Stand, "dimnames")[[2]])
}

# Check if names have been updated
print(attr(ldsc_results$S, "dimnames")[[2]])
print(attr(ldsc_results$S_Stand, "dimnames")[[2]])

# Define the model string
model_string <- paste0(
  "F1 =~ CXCL3_P19876_OID20788_v1_Inflammation_merged + CCL4_P13236_OID20695_v1_Inflammation_merged + CCL27_Q9Y4X3_OID20121_v1_Cardiometabolic_merged\n",
  "F1 ~ SNP\n"
)

# Step 6: Redefine SNPs after the renaming
SNPs2 <- summary_stats[, c("SNP", "CHR", "BP", "MAF", "A1", "A2", 
                           "beta.CXCL3_P19876_OID20788_v1_Inflammation_merged",
                           "se.CXCL3_P19876_OID20788_v1_Inflammation_merged",
                           "beta.CCL4_P13236_OID20695_v1_Inflammation_merged",
                           "se.CCL4_P13236_OID20695_v1_Inflammation_merged",
                           "beta.CCL27_Q9Y4X3_OID20121_v1_Cardiometabolic_merged",
                           "se.CCL27_Q9Y4X3_OID20121_v1_Cardiometabolic_merged")]

# Run userGWAS with SNPs2
CorrelatedFactors <- userGWAS(
  covstruc = ldsc_results, 
  SNPs = SNPs2, 
  model = model_string, 
  estimation = "DWLS", 
  printwarn = TRUE, 
  sub = c("F1~SNP"),
  toler = FALSE,
  SNPSE = FALSE,
  parallel = TRUE,
  GC = "standard",
  MPI = FALSE,
  smooth_check = TRUE,
  fix_measurement = TRUE,
  Q_SNP = TRUE
)

# Save userGWAS results
result_path <- "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/userGWAS_res/userGWAS_results_3_chemokines.rds"
saveRDS(CorrelatedFactors, file = result_path)
cat("GWAS completed. Results saved to:", result_path, "\n")

# Measure runtime
end_time <- Sys.time()
cat("Total runtime:", round(difftime(end_time, start_time, units = "mins"), 2), "minutes.\n")

CorrelatedFactors <- readRDS("/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/userGWAS_res/userGWAS_results_3_chemokines.rds")

# Optional: Filter low-quality SNPs
CF1_bad <- CorrelatedFactors[[1]][which(CorrelatedFactors[[1]]$Q_SNP_pval < 5e-8),]
write.table(CF1_bad$SNP, "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/userGWAS_res/poor_quality_snps.txt", row.names = FALSE, col.names = FALSE, quote = FALSE)

# Define low-quality SNPs and the full dataset
low_quality_snps <- CF1_bad  # Subset of SNPs with Q_SNP_pval < 5e-8
all_snps <- CorrelatedFactors[[1]]  # The full SNP dataset

# Parameters
ld_window_size <- 500000  # Proximity threshold in base pairs

# Filter low-quality SNPs for required columns
low_quality_snps <- low_quality_snps %>% select(SNP, CHR, BP)

# Exclude SNPs in LD with low-quality SNPs based on proximity
filtered_snps <- all_snps %>%
  filter(!sapply(1:n(), function(i) {
    any(abs(all_snps$BP[i] - low_quality_snps$BP[low_quality_snps$CHR == all_snps$CHR[i]]) <= ld_window_size)
  }))

# Save the filtered SNPs
write.table(filtered_snps, "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/userGWAS_res/F1_filtered_snps.tsv", row.names = FALSE, quote = FALSE, sep = "\t")

filtered_snps <- fread("/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/userGWAS_res/F1_filtered_snps.tsv")
filtered_snps <- filtered_snps[order(filtered_snps$CHR),]

# Output results before filtering
significant_snps1 <- filtered_snps[filtered_snps$Pval_Estimate < 5e-08, ]
low_heterogeneity1 <- significant_snps1[significant_snps1$Q_SNP_pval > 0.05, ]
snps_set <- significant_snps1$SNP
low_heterogeneity1_snps_set <- low_heterogeneity1$SNP

# Save the filtered SNPs
write.table(low_heterogeneity1, "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/userGWAS_res/low_heterogeneity1.tsv", row.names = FALSE, quote = FALSE, sep = "\t")

CXCL3_sig <- CXCL3[which(CXCL3$P<5E-8),]
latent_CXCL3_sig <- CXCL3_sig[which(CXCL3_sig$rsid%in%low_heterogeneity1_snps_set),]
CCL4_sig <- CCL4[which(CCL4$P<5E-8),]
latent_CCL4_sig <- CCL4_sig[which(CCL4_sig$rsid%in%low_heterogeneity1_snps_set),]
CCL27_sig <- CCL27[which(CCL27$P<5E-8),]
latent_CCL27_sig <- CCL27_sig[which(CCL27_sig$rsid%in%low_heterogeneity1_snps_set),]

`%notin%` <- Negate(`%in%`)
notin_CXCL3 <- low_heterogeneity1_snps_set[low_heterogeneity1_snps_set %notin% CXCL3_sig$rsid]
notin_CCL4 <- low_heterogeneity1_snps_set[low_heterogeneity1_snps_set %notin% CCL4_sig$rsid]
notin_CCL27 <- low_heterogeneity1_snps_set[low_heterogeneity1_snps_set %notin% CCL27_sig$rsid]
not_detected_in_CXCL3 <- CXCL3[CXCL3$rsid%in%notin_CXCL3]
not_detected_in_CCL4 <- CCL4[CCL4$rsid%in%notin_CCL4]
not_detected_in_CCL27 <- CCL27[CCL27$rsid%in%notin_CCL27]

1.13 Value

nohup Rscript /Users/charleenadams/ukbppp/comp_manhattan.R > comp_manhattan.log 2>&1 &

## CXCL3 Manhattan Plot

## CCL4 Manhattan Plot

## CCL27 Manhattan Plot

## Latent Factor Manhattan Plot (no filtering)

## Stringent QC: Labeled Latent Factor Manhattan Plot Manhattan Plot

1.13.1 Comparison with Original GWAS’s P-values

1.13.1.1 Code for P-values Comparison

Show

low_heterogeneity1 <- fread("/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/userGWAS_res/low_heterogeneity1.tsv")

# Table
low_heterogeneity1$SNP
low_heterogeneity1$rsid <- low_heterogeneity1$SNP
low_heterogeneity1$P <- low_heterogeneity1$Pval_Estimate
low_heterogeneity1$Dataset="multivariate"

CXCL3=fread("/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/fixed_files/header/fixed_CXCL3_P19876_OID20788_v1_Inflammation_merged_header_good.gz", fill = TRUE)
CCL4=fread("/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/fixed_files/header/fixed_CCL4_P13236_OID20695_v1_Inflammation_merged_header_good.gz", fill = TRUE)
CCL27=fread("/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/fixed_files/header/fixed_CCL27_Q9Y4X3_OID20121_v1_Cardiometabolic_merged_header_good.gz", fill = TRUE)

CXCL3$Dataset="CXCL3"
CCL4$Dataset="CCL4"
CCL27$Dataset="CCL27"

# Get the 7 SNPs from low_heterogeneity1
snp_list <- low_heterogeneity1$rsid

# Subset each dataset for the 7 SNPs and select the required columns
low_heterogeneity1_snps <- low_heterogeneity1[rsid %in% snp_list, .(Dataset, rsid, P)]
CXCL3_snps <- CXCL3[rsid %in% snp_list, .(Dataset, rsid, P)]
CCL4_snps <- CCL4[rsid %in% snp_list, .(Dataset, rsid, P)]
CCL27_snps <- CCL27[rsid %in% snp_list, .(Dataset, rsid, P)]

# Combine all datasets
combined_snps <- rbind(low_heterogeneity1_snps, CXCL3_snps, CCL4_snps, CCL27_snps)

# Save to a CSV file
write.csv(combined_snps, "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/ldsc_results_rgs/combined_snps_dataset.csv", row.names = FALSE)

1.13.2 P-value Figures

1.13.2.1 Code for P-value Figures

Show

# Load required libraries
library(ggplot2)
library(data.table)
library(gridExtra)

combined_snps <- fread("/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/ldsc_results_rgs/combined_snps_dataset.csv")

combined_snps[, neg_log10P := -log10(P)]

# Define genome-wide significance threshold
genomewide_threshold <- -log10(5e-8)

# Create Boxplot
boxplot_plot <- ggplot(combined_snps, aes(x = Dataset, y = neg_log10P, fill = Dataset)) +
  geom_boxplot() +
  geom_hline(yintercept = genomewide_threshold, linetype = "dotted", color = "red") +
  labs(
    title = "Comparison of -log10(P) Across Datasets",
    x = "Dataset",
    y = "-log10(P)"
  ) +
  theme_minimal()

# Create Scatterplot
scatterplot_plot <- ggplot(combined_snps, aes(x = rsid, y = neg_log10P, color = Dataset)) +
  geom_point(size = 3) +
  geom_hline(yintercept = genomewide_threshold, linetype = "dotted", color = "red") +
  labs(
    title = "Comparison of -log10(P) for SNPs Across Datasets",
    x = "SNP (rsid)",
    y = "-log10(P)"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Display the plots side by side
grid.arrange(boxplot_plot, scatterplot_plot, ncol = 2)

# Save the plots as 600 DPI
ggsave("boxplot_pvalues_across_datasets_with_threshold_600dpi.png", plot = boxplot_plot, width = 8, height = 6, dpi = 600)
ggsave("scatterplot_snp_pvalues_across_datasets_with_threshold_600dpi.png", plot = scatterplot_plot, width = 10, height = 6, dpi = 600)

1.14 Nearest Genes from Ensembl: CD36 and GCKR

1.14.1 Code for Nearest Genes from Ensembl

Show

library(biomaRt)

# Connect to Ensembl GRCh37 archive
ensembl_grch37 <- useMart(
  biomart = "ENSEMBL_MART_ENSEMBL",
  dataset = "hsapiens_gene_ensembl",
  host = "grch37.ensembl.org"
)

# Confirm the build
listDatasets(ensembl_grch37)

# Load your SNP dataset
filtered_snps <- fread("/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/userGWAS_res/F1_filtered_snps.tsv")
filtered_snps <- filtered_snps[order(filtered_snps$CHR), ]
filtered_snps <- filtered_snps[which(filtered_snps$Q_SNP_pval > 0.05), ]

# Query nearest genes using GRCh37
snps <- unique(filtered_snps[, .(CHR, BP, SNP)])

# Retrieve gene annotations
gene_annotations <- getBM(
  attributes = c("chromosome_name", "start_position", "end_position", 
                 "hgnc_symbol", "ensembl_gene_id"),
  filters = "chromosomal_region",
  values = paste0(snps$CHR, ":", snps$BP - 100000, ":", snps$BP + 100000),
  mart = ensembl_grch37
)

# Ensure filtered_snps and gene_annotations are data.tables
setDT(filtered_snps)
setDT(gene_annotations)

# Rename columns in gene_annotations to match filtered_snps for easier joining
setnames(gene_annotations, c("chromosome_name", "start_position", "end_position", "hgnc_symbol"), 
         c("CHR", "start", "end", "gene"))

# Find nearest genes for each SNP
filtered_snps[, nearest_gene := gene_annotations[.SD, 
                                                 on = .(CHR, start <= BP, end >= BP), 
                                                 mult = "first",  # Take the first match if multiple
                                                 x.gene]]

# Save the filtered SNPs
write.table(filtered_snps, "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/userGWAS_res/filtered_snps_q05_nearest_genes.tsv", row.names = FALSE, quote = FALSE, sep = "\t")

1.14.2 CD36: Cluster of Differentiation 36

• Location: Chromosome 7q21.11
  
• Function: Encodes a transmembrane glycoprotein involved in:
  • Fatty Acid Uptake: Critical for lipid metabolism.
    
  • Lipid Metabolism: Regulates lipid storage and utilization.
    
  • Immune Response: Contributes to inflammation and pathogen recognition.
    
  • Angiogenesis: Supports blood vessel formation.
    
• Clinical Relevance:
  • Implicated in insulin resistance, atherosclerosis, and cancer progression.
    
  • A key target for therapies addressing metabolic syndromes and cardiovascular diseases.

---

1.14.3 GCKR: Glucokinase Regulatory Protein

• Location: Chromosome 2p23.3
  
• Function: Encodes a protein that modulates glucokinase activity, playing a pivotal role in:
  • Glucose Metabolism: Regulates glucose levels in the liver and pancreas.
    
• Clinical Relevance:
  • Variants influence fasting glucose levels and triglyceride concentrations.
    
  • Associated with the risk of type 2 diabetes and non-alcoholic fatty liver disease (NAFLD).

---

1.14.4 Integrated Role in Metabolic Disorders

Together, Cluster of Differentiation 36 (CD36) and Glucokinase Regulatory Protein (GCKR) are critical to understanding the genetic architecture of metabolic disorders. These genes provide valuable insights for developing therapeutic strategies targeting:
- Diabetes
- Dyslipidemia
- Cardiovascular Diseases
- NAFLD

1.14.5 Code for Fetching Nearest Genes Bioinformatically

Show

library(biomaRt)

# Connect to Ensembl GRCh37 archive
ensembl_grch37 <- useMart(
  biomart = "ENSEMBL_MART_ENSEMBL",
  dataset = "hsapiens_gene_ensembl",
  host = "grch37.ensembl.org"
)

# Confirm the build
listDatasets(ensembl_grch37)

# Load filtered SNP dataset
filtered_snps <- fread("/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/userGWAS_res/F1_filtered_snps.tsv")
filtered_snps <- filtered_snps[order(filtered_snps$CHR), ]
filtered_snps <- filtered_snps[which(filtered_snps$Q_SNP_pval > 0.05), ]

# Query nearest genes using GRCh37
snps <- unique(filtered_snps[, .(CHR, BP, SNP)])

# Retrieve gene annotations
gene_annotations <- getBM(
  attributes = c("chromosome_name", "start_position", "end_position", 
                 "hgnc_symbol", "ensembl_gene_id"),
  filters = "chromosomal_region",
  values = paste0(snps$CHR, ":", snps$BP - 100000, ":", snps$BP + 100000),
  mart = ensembl_grch37
)

# Ensure filtered_snps and gene_annotations are data.tables
setDT(filtered_snps)
setDT(gene_annotations)

# Rename columns in gene_annotations to match filtered_snps for easier joining
setnames(gene_annotations, c("chromosome_name", "start_position", "end_position", "hgnc_symbol"), 
         c("CHR", "start", "end", "gene"))

# Find nearest genes for each SNP
filtered_snps[, nearest_gene := gene_annotations[.SD, 
                                                 on = .(CHR, start <= BP, end >= BP), 
                                                 mult = "first",  # Take the first match if multiple
                                                 x.gene]]

# Save the filtered SNPs
write.table(filtered_snps, "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/userGWAS_res/filtered_snps_q05_nearest_genes.tsv", row.names = FALSE, quote = FALSE, sep = "\t")


# Load SNPs to label
low_heterogeneity1 <- fread("/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/userGWAS_res/low_heterogeneity1.tsv")
label_snps <- low_heterogeneity1$SNP
tops <- filtered_snps[which(filtered_snps$SNP%in%label_snps),]

sigs <- filtered_snps[which(filtered_snps$Pval_Estimate<5E-8),]
sigs

1.15 CADD (Combined Annotation Dependent Depletion)

• Definition: CADD is a computational tool that scores the deleteriousness (harmfulness) of single nucleotide variants (SNVs) and small insertions/deletions (indels) in the genome.

• Purpose: It combines multiple annotations (e.g., conservation, regulatory impact, protein function, etc.) into a single score to predict how damaging a variant might be.

• Uses: Phred Quality Score (PHRED), which provides a scaled score indicating the relative deleteriousness of a variant.

SNP	CHR	BP	MAF	A1	A2	est	SE	Pval_Estimate	nearest_gene	Consequence	ConsDetail	GeneName	PHRED	PhastCons	PhyloP
rs1728918	2	27635463	0.272366	A	G	0.04235485	0.006950637	1.103550e-09	PPM1G	UPSTREAM	upstream	PPM1G	2.458	0.060	0.283
rs1260326	2	27730940	0.410537	T	C	0.04512821	0.006290023	7.253483e-13	GCKR	NON_SYNONYMOUS	splice,missense	GCKR	13.220	0.995	0.943
rs780094	2	27741237	0.410537	T	C	0.04275407	0.006290005	1.067132e-11	GCKR	INTRONIC	intron	GCKR	1.852	0.000	-0.872
rs780093	2	27742603	0.411531	T	C	0.04217555	0.006287704	1.978192e-11	GCKR	INTRONIC	intron	GCKR	1.541	0.000	-0.109
rs7779873	7	80211423	0.453280	G	A	-0.03564340	0.006215646	9.782121e-09	CD36	INTRONIC	intron	CD36	5.362	0.002	0.173
rs6961069	7	80218961	0.419483	C	T	-0.03646031	0.006270293	6.071746e-09	CD36	INTRONIC	intron	CD36	5.150	0.000	0.382
rs13236689	7	80236014	0.424453	T	G	-0.03506074	0.006260322	2.137726e-08	CD36	INTRONIC	intron	CD36	8.262	0.005	-0.008

1.15.1 Interpretations

• CADD-PHRED Score:
  • PHRED scores indicate the relative deleteriousness of variants.
  • Scores above 10: Variant is in the top 10% of deleterious variants.
  • Scores above 20: Variant is in the top 1% of deleterious variants.
• PhastCons:
  • Represents evolutionary conservation scores, with higher values indicating higher conservation across species.
  • \(0\) to \(1\)
  • 0: No conservation.
  • 1: Maximum conservation, indicating a highly conserved region across species.
• PhyloP:
  • Measures evolutionary conservation or acceleration at specific sites.
  • Positive scores indicate conservation, while negative scores indicate acceleration.
  • Negative to Positive (\(-\infty\) to \(+\infty\)).
  • Most scores typically fall within \(-3\) to \(+3\).
• Conclusion: The analyzed variants exhibit varying levels of functional significance based on CADD PHRED scores and conservation metrics. Variants such as rs1260326 (PHRED: 13.22) are of higher interest due to their non-synonymous coding impact and proximity to the gene GCKR. Variants within CD36 show low conservation but may play roles in regulatory mechanisms.

1.15.2 Protein Phosphatase, Mg²⁺/Mn²⁺ Dependent 1G (PPM1G)

PPM1G encodes a serine/threonine phosphatase that is part of the protein phosphatase 2C (PP2C) family. It plays a critical role in:
- Dephosphorylation: Removes phosphate groups from serine and threonine residues on proteins.
- Cell Cycle Regulation: Supports proper cell cycle progression.
- Pre-mRNA Splicing: Ensures accurate and efficient gene expression by regulating pre-mRNA splicing.
- Cell Stress Response: Dephosphorylates proteins involved in stress adaptation.
- Cancer Biology: Altered expression has been linked to tumorigenesis.
- Neurological Function: Plays a role in neuronal signaling and brain activity.

1.15.3 Code CADD/PHRED

Show

# CADD from UW website
library(data.table)

# Example dataset
sigs <- data.table(
  CHR = c(2, 2, 2, 2, 7, 7, 7),
  BP = c(27635463, 27730940, 27741237, 27742603, 80211423, 80218961, 80236014),
  ID = c("rs1728918", "rs1260326", "rs780094", "rs780093", "rs7779873", "rs6961069", "rs13236689"),
  A1 = c("A", "T", "T", "T", "G", "C", "T"),
  A2 = c("G", "C", "C", "C", "A", "T", "G")
)

# Rename columns to match VCF format
setnames(sigs, c("CHR", "BP", "ID","A2", "A1"), c("#CHROM", "POS", "ID","REF", "ALT"))

# Save as a VCF file (without a header, as CADD ignores it)
fwrite(sigs, "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/userGWAS_res/variants_for_cadd.vcf", sep = "\t", col.names = FALSE, quote = FALSE)

# Upload to https://cadd.gs.washington.edu/upload as .gz and select GRCh37-v1.4

# Download resulting gz and mv it to: GRCh37-v1.4_anno_2ff4228ba73e3002aa20b4141c1e5f31.tsv 
cadd <- fread("/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/userGWAS_res/GRCh37-v1.4_anno_2ff4228ba73e3002aa20b4141c1e5f31.tsv")

# Load SNPs to label
low_heterogeneity1 <- fread("/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/userGWAS_res/low_heterogeneity1.tsv")
label_snps <- low_heterogeneity1$SNP
tops <- filtered_snps[which(filtered_snps$SNP%in%label_snps),]

# Ensure both datasets are data.tables
setDT(cadd)
setDT(tops)

# Create a matching key in both datasets
tops[, cadd_key := paste(CHR, BP, sep = "-")]
cadd[, cadd_key := paste(`#Chrom`, Pos, sep = "-")]

# Merge on the generated key
merged_data <- merge(tops, cadd, by = "cadd_key", all.x = TRUE)

# Print the merged data to verify
print(merged_data, max.print=T)

# Save the merged data to a file for future reference
output_path <- "/Users/charleenadams/temp_BI/chemokine_rgs_olink/processed_ukbppp_chemokine_list/munged_40/merged_munged_protein_results/pvals_fix/userGWAS_res/merged_cadd_tops_results.csv"
fwrite(merged_data, output_path)

1.16 Functional Mapping and Annotation (FUMA) of Genome-Wide Association Studies

• FUMA is a platform that can be used to annotate, prioritize, visualize and interpret GWAS results. I plugged our multivariate summary statistics into FUMA and got various annotations. I requested GTEx v8 expression values for our two top genes.

1.16.1 GTEx Heatmap

## GTEx Heatmap

1.16.2 “Indepdendent” SNPs (r2 < 0.6) at the Chr 2 and Chr 7 Loci

FUMA’s default for independence is Plink’s clumping at r2 < 0.6. This is user-defined. I could have specified a stricter independence threshold, but these are “indepedent at r2 < 0.6.

Genomic Locus	uniqID	rsID	chr	pos
1	2:27635463:A:G	rs1728918	2	27635463
1	2:27730940:C:T	rs1260326	2	27730940
2	7:80211423:A:G	rs7779873	7	80211423
2	7:80218961:C:T	rs6961069	7	80218961

1.17 LD Pair Analysis from LDLink

I wanted to know the actual r2 between those identified with the r2 < 0.6 threshold in FUMA as “independent.” I used LD Pair:

Show

---

Query SNPs

• rs1728918 (chr2:27635463)
  
• rs1260326 (chr2:27730940)

EUR Haplotypes

	rs1260326
	C	T	Total
rs1728918
A	12	262	274 (0.272)
G	581	151	732 (0.728)
Total	593 (0.589)	413 (0.411)	1006

Haplotypes

• G_C: 581 (0.578)
  
• A_T: 262 (0.26)
  
• G_T: 151 (0.15)
  
• A_C: 12 (0.012)

Summary Statistics

• D’: 0.9257
  
• R²: 0.4606
  
• Chi-sq: 463.3242
  
• p-value: <0.0001

Observations

• rs1728918(A) allele is correlated with rs1260326(T) allele.
  
• rs1728918(G) allele is correlated with rs1260326(C) allele.

---

Query SNP

• rs7779873 (chr7:80211423)
  
• rs6961069 (chr7:80218961)

EUR Haplotypes

	rs6961069
	C	T	Total
rs7779873
A	72	384	456 (0.453)
G	512	38	550 (0.547)
Total	584 (0.581)	422 (0.419)	1006

Haplotypes

• G_C: 512 (0.509)
  
• A_T: 384 (0.382)
  
• A_C: 72 (0.072)
  
• G_T: 38 (0.038)

Summary Statistics

• D’: 0.8353
  
• R²: 0.6081
  
• Chi-sq: 611.7519
  
• p-value: <0.0001

Observations

• rs7779873(A) allele is correlated with rs6961069(T) allele.
  
• rs7779873(G) allele is correlated with rs6961069(C) allele.