基于R计算各疾病队列中每种致病性胚系变异的频率

前置说明 #

以下代码默认你数据中：

• 变异唯一标识列名为VariantID（即你举例的1:17588689这类位点ID，若实际为多列组合标识，可先拼接为唯一ID再使用）
• 疾病信息表中疾病类型列名为DiseaseType
• 所有未匹配到疾病信息、或疾病类型不在指定4类中的样本统一归入others组，最终5类占比和严格为1

完整可运行代码 #

你已加载的3个依赖包足够，不需要额外安装新包：

library(stringr)
library(tidyr)
library(dplyr)

# 1. 读取数据（保留原始宽表variantData_raw，不要提前拆长表，方便最后合并结果）
diseaseData <- read.delim(".../disease_cohort.txt", header = T, sep = "\t") 
variantData_raw <- read.delim(".../variant_list.txt", header = T, sep = "\t")

# 2. 处理变异表：拆分杂合样本ID为长表
variantData_long <- variantData_raw %>% 
  mutate(HetSamples = strsplit(as.character(HetSamples), ",")) %>% 
  unnest(HetSamples) %>%
  # 过滤HetSamples为空的无效记录
  filter(HetSamples != "")

# 3. 统一新旧队列ID匹配规则（拆分ID操作本质是对齐两边ID格式，可根据实际ID规则调整）
# 旧队列ID拆分对齐
variantOld <- variantData_long %>% 
  filter(!str_detect(HetSamples, 'U'))
diseaseOld <- diseaseData %>% 
  filter(!str_detect(ClinicalSeqID, 'U'))
variantOld[c('Col1', 'Col2', 'Col3')] <- str_split_fixed(variantOld$HetSamples, '-', 3)
diseaseOld[c('Col1', 'Col2', 'Col3')] <- str_split_fixed(diseaseOld$ClinicalSeqID, '-', 3)
# 生成匹配用ID，若直接用原始HetSamples/ClinicalSeqID就能匹配，可替换match_id生成逻辑
variantOld <- variantOld %>% mutate(match_id = HetSamples)
diseaseOld <- diseaseOld %>% mutate(match_id = ClinicalSeqID)

# 新队列ID不需要拆分，直接生成匹配ID
variantNew <- variantData_long %>% 
  filter(str_detect(HetSamples, 'U'))
diseaseNew <- diseaseData %>% 
  filter(str_detect(ClinicalSeqID, 'U'))
variantNew <- variantNew %>% mutate(match_id = HetSamples)
diseaseNew <- diseaseNew %>% mutate(match_id = ClinicalSeqID)

# 4. 合并新旧数据集，关联样本对应的疾病类型
variantAll <- bind_rows(variantOld, variantNew)
diseaseAll <- bind_rows(diseaseOld, diseaseNew) %>% select(match_id, DiseaseType)

variantWithDisease <- variantAll %>%
  left_join(diseaseAll, by = "match_id") %>%
  # 去重：避免同一个样本同个变异被重复计数
  distinct(VariantID, match_id, .keep_all = T) %>%
  # 疾病类型归类，不在指定4类/未匹配到的统一归为others
  mutate(DiseaseType = case_when(
    is.na(DiseaseType) ~ "others",
    !DiseaseType %in% c("glioma", "meningioma", "schwannoma", "pituitary adenoma") ~ "others",
    TRUE ~ DiseaseType
  ))

# 5. 按变异位点统计各疾病队列占比
targetDiseases <- c("meningioma", "glioma", "schwannoma", "pituitary adenoma", "others")
variantFreq <- variantWithDisease %>%
  group_by(VariantID, DiseaseType) %>%
  summarise(sampleCount = n(), .groups = "drop") %>%
  # 补全每个位点所有疾病类型的计数，无样本的填0
  complete(VariantID, DiseaseType = targetDiseases, fill = list(sampleCount = 0)) %>%
  group_by(VariantID) %>%
  mutate(
    totalSample = sum(sampleCount),
    # 计算占比，保留6位小数消除浮点误差
    freq = round(sampleCount / totalSample, 6)
  ) %>%
  # 转为宽表，每列对应一个疾病队列的占比
  select(VariantID, DiseaseType, freq) %>%
  pivot_wider(names_from = DiseaseType, values_from = freq, names_prefix = "freq_")

# 6. 将占比结果追加到原始变异宽表最后
finalData <- variantData_raw %>%
  left_join(variantFreq, by = "VariantID")

# 可选：校验占比和是否为1
checkSum <- finalData %>%
  mutate(total_freq = rowSums(select(., starts_with("freq_")), na.rm = T)) %>%
  filter(abs(total_freq - 1) > 0.001)
if(nrow(checkSum) == 0){
  message("所有变异位点队列占比之和为1，统计无误")
}else{
  message("存在占比和不为1的位点，请检查ID匹配逻辑")
}

注意事项 #

• 代码中VariantID、DiseaseType列名请替换为你自己数据中对应的实际列名，如果变异唯一标识是多列组合（如染色体+位置+参考等位基因+变异等位基因），请先拼接为唯一ID列再分组统计，避免统计错误。
• 若你的ID匹配规则不是直接用原始样本ID，可修改match_id的生成逻辑，保证变异表中的样本ID和疾病表中的临床ID能一一对应即可。
• 最后输出的finalData就是追加了5类疾病队列占比列的完整数据集，占比列分别为freq_meningioma、freq_glioma、freq_schwannoma、freq_pituitary adenoma、freq_others。

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