inst/complete_workflow_example.R

In CoxMK: A Model-X Knockoff Method for Genome-Wide Survival Association Analysis

# CoxMK Package Complete Workflow Example
library(CoxMK)
library(survival)
library(Matrix)

cat("=== CoxMK Package Complete Workflow Example ===\n")
cat("Note: This example uses temporary directories to comply with CRAN policies.\n\n")

# Create temporary directory for outputs
temp_dir <- tempdir()
cat("Working in temporary directory:", temp_dir, "\n\n")

# 1. Load data files
cat("1. Loading PLINK and phenotype data...\n")
extdata_path <- system.file("extdata", package = "CoxMK")

# Read phenotype and covariate data directly
phenotype_data <- read.table(file.path(extdata_path, "tte_phenotype.txt"), 
                            header = TRUE, stringsAsFactors = FALSE)
covariate_data <- read.table(file.path(extdata_path, "covariates.txt"), 
                            header = TRUE, stringsAsFactors = FALSE)

cat("   - Phenotypes:", nrow(phenotype_data), "samples\n")
cat("   - Covariates:", nrow(covariate_data), "samples\n")
cat("   - Events observed:", sum(phenotype_data$status), "out of", nrow(phenotype_data), "samples\n\n")

# 2. Use cox_knockoff_analysis for complete workflow
cat("2. Running complete Cox knockoff analysis...\n")

# Prepare covariate matrix (exclude ID columns)
covariates <- covariate_data[, c("age", "sex", "bmi", "smoking")]

result <- cox_knockoff_analysis(
  plink_prefix = file.path(extdata_path, "sample"),
  time = phenotype_data$time,
  status = phenotype_data$status,
  covariates = covariates,
  sample_ids = phenotype_data$IID,
  M = 5,
  fdr = 0.05,
  output_dir = temp_dir
)

cat("   - Analysis completed successfully\n")
cat("   - Variables selected:", length(result$selected_vars), "\n")
cat("   - W statistics calculated for all variants\n\n")

# 3. Save results
cat("3. Saving analysis results...\n")

# Save results to RDS file
results_file <- file.path(temp_dir, "analysis_results.rds")
saveRDS(result, results_file)
cat("   - Analysis results saved to:", basename(results_file), "\n")

# Save selected SNP information if any selected
if (length(result$selected_vars) > 0) {
  selected_snps <- data.frame(
    index = result$selected_vars,
    W_statistic = result$W_stats[result$selected_vars]
  )
  
  selected_file <- file.path(temp_dir, "selected_snps.txt")
  write.table(selected_snps, selected_file, row.names = FALSE, quote = FALSE, sep = "\t")
  cat("   - Selected SNPs saved to:", basename(selected_file), "\n")
}
cat("\n")

# 4. Analysis Summary
cat("4. Analysis Summary\n")
cat("================================================================================\n")
cat("Results:\n")
cat("- Variables selected:", length(result$selected_vars), "\n")
if (length(result$W_stats) > 0) {
  cat("- Total variants tested:", length(result$W_stats), "\n")
  cat("- Selection rate:", round(100 * length(result$selected_vars) / length(result$W_stats), 2), "%\n")
}

cat("\n=== CoxMK Workflow Completed Successfully ===\n")
cat("Files saved to temporary directory:\n")
cat("✓ analysis_results.rds - Complete analysis results\n")
if (length(result$selected_vars) > 0) {
  cat("✓ selected_snps.txt - Selected SNP information\n")
}
cat("\nNote: Files are saved in temporary directory for CRAN compliance.\n")
cat("In practice, specify output_dir parameter to save to desired location.\n")