176_肿瘤HRD评分计算¶

一句话概述¶

HRD（Homologous Recombination Deficiency）评分通过三个基因组不稳定性指标（LOH、TAI、LST）的综合评估来量化肿瘤同源重组修复缺陷程度，是指导PARP抑制剂和铂类化疗用药的重要生物标志物。

核心知识点表格¶

知识点	说明
HRD定义	同源重组修复通路缺陷，导致基因组不稳定
HRD评分组成	LOH + TAI + LST 三个评分之和
LOH	Loss of Heterozygosity，杂合性丢失区段数
TAI	Telomeric Allelic Imbalance，端粒等位基因不平衡
LST	Large-scale State Transitions，大规模状态转换数
临床阈值	HRD≥42通常判定为HRD阳性（Myriad标准）
关联基因	BRCA1, BRCA2, RAD51, PALB2等HR通路基因
临床应用	PARP抑制剂（奥拉帕利等）用药伴随诊断
计算工具	scarHRD, sigtools, HRDetect

步骤详解¶

第一步：理解HRD的生物学背景¶

白话解释：细胞DNA损伤后需要修复。同源重组修复（HR）是最精确的修复方式，由BRCA1/2等基因控制。当这些基因出了问题，细胞只能用不太精确的方式修复，留下大量"疤痕"。HRD评分就是数这些疤痕的数量。

技术细节：HR通路缺陷导致DNA双链断裂使用非同源末端连接（NHEJ）或替代修复途径，产生三种特征性基因组改变：(1) LOH - 大片段杂合性丢失；(2) TAI - 延伸到端粒的等位基因不平衡；(3) LST - 相邻大片段拷贝数状态之间的转换。这三个指标分别从不同角度反映基因组不稳定性。

# HRD概念图解
DNA双链断裂 → HR修复正常？
    ├─ 是 → 精确修复 → 基因组稳定 → HRD评分低
    └─ 否（BRCA突变等）→ 错误修复 → 基因组疤痕 → HRD评分高
                                        ├─ LOH（杂合性丢失）
                                        ├─ TAI（端粒不平衡）
                                        └─ LST（状态转换）

第二步：准备输入数据¶

白话解释：HRD评分需要肿瘤样本的拷贝数变异（CNV）信息和B等位基因频率（BAF）数据。通常来自全基因组测序（WGS）、全外显子组测序（WES）或SNP芯片。

技术细节：需要逐段的拷贝数（total copy number）和等位基因特异性拷贝数（allele-specific copy number）。常用的上游工具包括ASCAT、Sequenza、FACETS、PureCN等。

# 方法1：从WGS/WES数据用Sequenza获取输入
# 安装Sequenza
# R: install.packages("sequenza")
# Python: pip install sequenza-utils

# Step1: 处理BAM文件
sequenza-utils bam2seqz \
    --normal normal.bam \
    --tumor tumor.bam \
    --fasta ref.fa \
    -gc gc50.wig.gz \
    --output sample.seqz.gz

# Step2: 分箱
sequenza-utils seqz_binning \
    --seqz sample.seqz.gz \
    -w 50 \
    --output sample.small.seqz.gz

# Step3: R中运行Sequenza
library(sequenza)

# 加载数据
seqz <- sequenza.extract(file = "sample.small.seqz.gz")

# 拟合参数（纯度和倍性）
cp <- sequenza.fit(seqz)

# 输出结果
sequenza.results(
  sequenza.extract = seqz,
  cp.table = cp,
  sample.id = "tumor_sample",
  out.dir = "sequenza_output"
)

第三步：使用scarHRD计算HRD评分¶

白话解释：scarHRD是最常用的HRD评分计算R包。你只需要把拷贝数分段数据输入进去，它就能帮你算出LOH、TAI、LST三个分数。

技术细节：scarHRD要求输入包含6列的分段文件：SampleID, Chromosome, Start, End, total_cn, A_cn（次等位基因拷贝数）。它会根据这些分段信息计算三个组分得分。

# 安装scarHRD
# install.packages("devtools")
# devtools::install_github("sztup/scarHRD")
library(scarHRD)

# 准备输入数据（来自Sequenza或其他工具）
# 格式：SampleID, Chromosome, Start, End, total_cn, A_cn
seg <- read.table("sequenza_output/tumor_sample_segments.txt", header = TRUE)

# 转换为scarHRD格式
scar_input <- data.frame(
  SampleID = "tumor_sample",
  Chromosome = seg$chromosome,
  Start_position = seg$start.pos,
  End_position = seg$end.pos,
  total_cn = seg$CNt,
  A_cn = seg$A        # minor allele copy number
)

# 计算HRD评分
hrd_result <- scar_score(
  scar_input,
  reference = "grch38",    # 参考基因组版本
  seqz = FALSE             # 是否直接从seqz文件读取
)

print(hrd_result)
# 输出示例：
#       HRD-LOH  Telomeric AI  LST  HRD
# [1,]    12        15          18   45

# 判定HRD状态
hrd_total <- hrd_result$HRD_sum  # LOH + TAI + LST
hrd_status <- ifelse(hrd_total >= 42, "HRD-positive", "HRD-negative")
cat("HRD Score:", hrd_total, "Status:", hrd_status, "\n")

第四步：使用Sequenza直接计算¶

# Sequenza内置的HRD评分
# 方法：直接从seqz文件计算
library(scarHRD)

hrd_from_seqz <- scar_score(
  "sample.small.seqz.gz",
  reference = "grch38",
  seqz = TRUE
)
print(hrd_from_seqz)

第五步：批量样本处理¶

# 批量处理多个样本
library(parallel)

sample_files <- list.files("segments/", pattern = "*.seg", full.names = TRUE)

process_hrd <- function(seg_file) {
  seg <- read.table(seg_file, header = TRUE, sep = "\t")
  sample_id <- gsub("_segments.seg", "", basename(seg_file))

  scar_input <- data.frame(
    SampleID = sample_id,
    Chromosome = seg$chrom,
    Start_position = seg$loc.start,
    End_position = seg$loc.end,
    total_cn = seg$tcn,
    A_cn = seg$lcn
  )

  tryCatch({
    result <- scar_score(scar_input, reference = "grch38", seqz = FALSE)
    data.frame(
      Sample = sample_id,
      LOH = result[1, 1],
      TAI = result[1, 2],
      LST = result[1, 3],
      HRD = result[1, 4]
    )
  }, error = function(e) {
    data.frame(Sample = sample_id, LOH = NA, TAI = NA, LST = NA, HRD = NA)
  })
}

# 并行计算
results <- mclapply(sample_files, process_hrd, mc.cores = 8)
hrd_table <- do.call(rbind, results)

# 添加HRD状态
hrd_table$Status <- ifelse(hrd_table$HRD >= 42, "HRD+", "HRD-")
write.csv(hrd_table, "hrd_scores_all_samples.csv", row.names = FALSE)

第六步：结果可视化¶

library(ggplot2)
library(tidyr)
library(dplyr)

# 1. HRD评分分布图
ggplot(hrd_table, aes(x = HRD, fill = Status)) +
  geom_histogram(bins = 30, alpha = 0.7) +
  geom_vline(xintercept = 42, linetype = "dashed", color = "red", size = 1) +
  annotate("text", x = 44, y = Inf, label = "Threshold = 42",
           hjust = 0, vjust = 1.5, color = "red") +
  theme_minimal() +
  labs(x = "HRD Score", y = "Count", title = "HRD Score Distribution")

# 2. 三组分堆叠图
hrd_long <- hrd_table %>%
  pivot_longer(cols = c(LOH, TAI, LST), names_to = "Component", values_to = "Score") %>%
  arrange(desc(HRD))

hrd_long$Sample <- factor(hrd_long$Sample, levels = unique(hrd_long$Sample))

ggplot(hrd_long, aes(x = Sample, y = Score, fill = Component)) +
  geom_bar(stat = "identity") +
  geom_hline(yintercept = 42, linetype = "dashed", color = "red") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 6)) +
  labs(y = "HRD Score", title = "HRD Components by Sample")

# 3. BRCA状态与HRD关系
# 假设有BRCA突变信息
ggplot(hrd_table, aes(x = BRCA_status, y = HRD, fill = BRCA_status)) +
  geom_boxplot(alpha = 0.7) +
  geom_jitter(width = 0.2, alpha = 0.5) +
  geom_hline(yintercept = 42, linetype = "dashed", color = "red") +
  theme_minimal() +
  labs(title = "HRD Score by BRCA Status")

实战命令速查¶

# 从BAM到HRD评分完整流程
# 1. 运行FACETS（替代Sequenza）
Rscript run_facets.R --tumor tumor.bam --normal normal.bam --out facets_output

# 2. 运行scarHRD
Rscript -e '
library(scarHRD)
result <- scar_score("segments.txt", reference="grch38", seqz=FALSE)
write.csv(result, "hrd_result.csv")
'

# 3. 使用CHORD（另一种HRD工具，基于突变特征）
Rscript -e '
library(CHORD)
contexts <- extractSigsChord(vcf_snv, vcf_indel, vcf_sv, ref_genome="BSgenome.Hsapiens.UCSC.hg38")
chord_result <- chordPredict(contexts)
'

# 4. 使用HRDetect（全面的HRD检测方法）
# 需要SNV、Indel、SV、CNV多种数据

面试常问点¶

Q1: HRD评分的三个组分分别反映什么？ A: LOH（杂合性丢失）反映大片段的一个等位基因丢失；TAI（端粒等位基因不平衡）反映延伸到端粒但不跨越着丝粒的等位基因不平衡，与断裂修复缺陷相关；LST（大规模状态转换）反映相邻大片段（≥10Mb）之间的拷贝数状态转换，与BRCA1缺陷高度相关。

Q2: HRD阈值42是怎么来的？ A: 阈值42由Myriad Genetics在myChoice CDx检测中建立，基于大规模临床试验数据，该阈值能最佳区分对PARP抑制剂有响应和无响应的患者。不同研究和平台可能采用略有不同的阈值。

Q3: HRD评分与BRCA突变是什么关系？ A: BRCA1/2突变是HRD的主要但非唯一原因。约50%的HRD阳性肿瘤有BRCA突变，其余可能由RAD51、PALB2等其他HR通路基因缺陷或BRCA启动子甲基化引起。HRD评分能捕获所有原因导致的HR缺陷。

Q4: WES和WGS数据算出的HRD评分有差异吗？ A: 有差异。WGS能更精确地检测大片段CNV和LOH，HRD评分更准确。WES由于覆盖间断，可能低估LST和TAI。实际中WES数据也可使用，但需注意平台间的可比性。

Q5: scarHRD和商业检测（myChoice）结果一致吗？ A: 文献报告两者相关性较好（r>0.9），但不完全一致。商业检测使用SNP芯片数据和专有算法，scarHRD基于NGS数据。临床决策建议以FDA批准的商业检测为准。

Q6: HRD评分在什么癌种中最有价值？ A: 卵巢癌中应用最成熟（FDA批准奥拉帕利伴随诊断），乳腺癌（特别是三阴性）也有明确临床价值。前列腺癌、胰腺癌等也有研究支持。

易错点¶

参考基因组不匹配：scarHRD中reference参数必须与数据的参考基因组一致（hg19/grch37 vs hg38/grch38）
次等位基因拷贝数错误：A_cn应该是minor allele的拷贝数，不是major allele
倍性估计不准：上游工具的倍性和纯度估计会直接影响HRD评分，应仔细检查
忽略肿瘤纯度：低纯度样本的CNV检测不准确，导致HRD评分偏低
性染色体处理：HRD评分只计算常染色体，确保输入不包含X/Y染色体

补充知识¶

HRDetect方法¶

HRDetect使用加权逻辑回归模型整合6种突变特征来预测HRD： - SBS3（单碱基替换特征3，与HR缺陷相关） - SBS8 - Rearrangement signatures 3和5 - HRD-LOH指数 - 小缺失中的微同源性

# HRDetect需要更多类型的变异数据
# 优势是在BRCA突变之外也能检测HRD
# 被认为是目前最全面的HRD检测方法

各工具适用场景¶

工具	输入数据	方法	优势
scarHRD	CNV segments	基因组疤痕	简单快速
CHORD	SNV+Indel+SV	突变特征	不依赖CNV
HRDetect	多种变异	综合模型	最全面
myChoice CDx	SNP芯片	商业专有	FDA批准