在R中无染色体位置数据时如何开展GWAS分析？

方法一：使用statgenGWAS构建gData对象运行GWAS #

你的代码报错核心原因是statgenGWAS要求输入gData格式对象，且初始数据因cbind被强制转换为字符型，需先修正这两个问题：

1. 修正数据类型，避免字符型转换

# SNPs are in A/B notation, with 0 = AA, 1 = AB, and 2 = BB
library(statgenGWAS)

id <- c("person1", "person2", "person3", "person4", "person5", "person6", "person7", "person8", "person9", "person10")
snp1 <- c(0, 1, 2, 2, 1, 0, 0, 0, 1, 1)
snp2 <- c(2, 2, 2, 1, 1, 1, 0, 0, 0, 1)
snp3 <- c(0, 0, 2, 2, 0, 2, 1, 0, 2, 2)
diagnosis <- c(0, 1, 1, 0, 0, 1, 1, 0, 1, 1)

# 直接创建数据框，确保数值型变量类型正确
data <- data.frame(
  id = id,
  snp1 = snp1,
  snp2 = snp2,
  snp3 = snp3,
  diagnosis = as.factor(diagnosis) # 二元结局转因子，适配逻辑回归模型
)

2. 导入并关联SNP位置数据

# 示例SNP位置数据，替换为你实际的文件读取代码（如read.csv）
snp_info <- data.frame(
  snp = c("snp1", "snp2", "snp3"),
  chr = c(1, 2, 3), # 替换为实际染色体编号
  pos = c(1000, 2000, 3000) # 替换为实际物理位置
)

3. 构建gData对象并运行GWAS

# 分离基因型与表型数据
geno_data <- data[, c("id", "snp1", "snp2", "snp3")]
pheno_data <- data[, c("id", "diagnosis")]

# 创建符合要求的gData对象
gData <- createGData(
  geno = geno_data,
  pheno = pheno_data,
  map = snp_info
)

# 针对二元性状，运行逻辑回归模型的GWAS
gwas1a <- runSingleTraitGwas(
  gData = gData,
  traits = "diagnosis",
  model = "logistic"
)

# 查看分析结果
summary(gwas1a)

方法二：无需statgenGWAS的轻量分析方法 #

如果不想依赖特定包的对象格式，可通过以下两种方式快速完成分析：

方式1：手动循环实现逻辑回归 #

# 确保数据类型正确（同方法一的data创建）
data <- data.frame(
  id = id,
  snp1 = snp1,
  snp2 = snp2,
  snp3 = snp3,
  diagnosis = as.factor(diagnosis)
)

# 提取所有SNP列名
snp_cols <- grep("snp", colnames(data), value = TRUE)

# 循环拟合每个SNP与结局的逻辑回归，提取关键结果
gwas_results <- lapply(snp_cols, function(snp) {
  fit_formula <- as.formula(paste("diagnosis ~", snp))
  glm_fit <- glm(fit_formula, data = data, family = binomial)
  fit_summary <- summary(glm_fit)
  
  data.frame(
    SNP = snp,
    OR = exp(coef(glm_fit)[2]), # 计算优势比
    P_value = fit_summary$coefficients[2, 4],
    stringsAsFactors = FALSE
  )
})

# 合并结果并按p值排序
gwas_results_df <- do.call(rbind, gwas_results)
gwas_results_df <- gwas_results_df[order(gwas_results_df$P_value), ]

# 添加显著性校正判断
gwas_results_df$Bonferroni_Threshold <- 0.05 / nrow(gwas_results_df)
gwas_results_df$Significant_Bonferroni <- gwas_results_df$P_value < gwas_results_df$Bonferroni_Threshold
gwas_results_df$FDR <- p.adjust(gwas_results_df$P_value, method = "fdr")
gwas_results_df$Significant_FDR <- gwas_results_df$FDR < 0.05

# 查看最终结果
print(gwas_results_df)

方式2：使用SNPassoc包简化分析 #

library(SNPassoc)

# 准备数据（同方法一的data创建）
data <- data.frame(
  id = id,
  snp1 = snp1,
  snp2 = snp2,
  snp3 = snp3,
  diagnosis = as.factor(diagnosis)
)

# 初始化SNP分析对象
snp_data <- setupSNP(
  data,
  colSNPs = grep("snp", colnames(data), value = TRUE),
  status = "diagnosis"
)

# 运行二元性状GWAS分析
gwas_fit <- snp.rhs.tests(diagnosis ~ 1, data = snp_data, model = "logistic")

# 查看包含校正p值的结果
summary(gwas_fit)

筛选显著SNPs的常用方法 #

• Bonferroni校正：将显著性阈值设为0.05 / 总SNP数，严格控制家族错误率，适合小数据集
• FDR校正：通过p.adjust(..., method = "fdr")计算错误发现率，平衡性更好，适合中等规模数据集
• 可视化验证：后续可将筛选出的显著SNP与染色体位置数据结合，用qqman包绘制曼哈顿图或QQ图，直观判断关联信号

内容的提问来源于stack exchange，提问作者Scott Hebert