学习笔记-如何对SNP进行Hardy-Weinberg检验

这个分析针对那种列名为各个SNP的名称，行名为各个样本的基因型的数据集开展Hardy-Weinberg检验。

也就是大概长成下面这个样子的数据（只是一个示意图）：

	rs561698	rs561699	rs561700	rs561701	rs561802
sample1	0	1	2	1	0
sample2	0	0	0	0	2

主要是利用卡方检验来对该数据是否符合Hardy-Weinberg平衡，也就是我们的数据集本身是一个Case的组（在代码里我们称其成为x组），而我们要去构造一个符合Hardy-Weinberg平衡的Control组（在代码里我们称其为y组），使用卡方检验探究两组数据之间是否存在差异（p>0.05就可以简单粗暴地说两组没有显著差异，也就是我们的数据集符合Hardy-Weinberg平衡，就可以开展后续进一步的分析了！）

主要的难点还是在于处理数据集，在算Hardy-Weinberg时并没有什么难度。就是把p和q算出来，根据p和q构造出Expected的数据就行。

以下是完整的代码。

library(readr)
library(dplyr)

# 读取数据
df <- read_csv("你的csv文件路径和名称")

# 初始化存储卡方检验结果的向量
result_columns <- c()

# 遍历数据集的每一列
for (col_name in names(df)) {
  # 计算每列各个数值的计数
  value_counts <- df %>%
    count(!!sym(col_name)) %>%
    rename(Value = !!sym(col_name), Count = n)
  
  # 处理缺少的值
  value_1_count <- value_counts %>%
    filter(Value == 0) %>%
    pull(Count)
  value_2_count <- value_counts %>%
    filter(Value == 1) %>%
    pull(Count)
  value_3_count <- value_counts %>%
    filter(Value == 2) %>%
    pull(Count)
  
  # 如果缺少某些值，设置计数为 0
  value_1_count <- ifelse(length(value_1_count) == 0, 0, value_1_count)
  value_2_count <- ifelse(length(value_2_count) == 0, 0, value_2_count)
  value_3_count <- ifelse(length(value_3_count) == 0, 0, value_3_count)
  
  # 重新计算总数
  total_count <- value_1_count + value_2_count + value_3_count
  
  # 如果总数为0，则跳过该列
  if (total_count == 0) next
  
  # 计算 q 和 p
  q <- (2 * value_1_count + value_2_count) / (2 * total_count)
  p <- (2 * value_3_count + value_2_count) / (2 * total_count)
  
  # 计算 ExpectedNumber 列
  # 注意: 只计算存在的值
  expected_counts <- numeric(3)
  
  if (value_1_count > 0) {
    expected_counts[1] <- (q^2) * total_count
  }
  if (value_2_count > 0) {
    expected_counts[2] <- 2 * p * q * total_count
  }
  if (value_3_count > 0) {
    expected_counts[3] <- (p^2) * total_count
  }
  
  # 对应的实际计数
  observed_counts <- c(value_1_count, value_2_count, value_3_count)
  
  # 确保 `expected_counts` 中对应的值存在
  valid_indices <- which(expected_counts > 0)
  observed_counts <- observed_counts[valid_indices]
  expected_counts <- expected_counts[valid_indices]
  
  # 计算概率分布
  if (length(expected_counts) > 0) {
    expected_prob <- expected_counts / sum(expected_counts)
    
    # 进行 Hardy-Weinberg 卡方检验
    test_result <- chisq.test(x = observed_counts, p = expected_prob)
    
    # 提取 p-value
    p_value <- test_result$p.value
    
    # 如果 p_value <= 0.05，将列名添加到结果向量
    if (p_value <= 0.05) {
      result_columns <- c(result_columns, col_name)
    }
  }
}

# 将结果写入文件
result_df <- tibble(ColumnName = result_columns)
write_csv(result_df, "你的结果输出路径和名称")