166_DNA损伤修复通路分析¶

# === 使用scarHRD计算HRD评分 ===
# install.packages("devtools")
# devtools::install_github("sztup/scarHRD")
library(scarHRD)

# 输入: Sequenza输出的拷贝数segment文件
# 格式: SampleID, Chromosome, Start_position, End_position, 
#        total_cn, A_cn, B_cn, ploidy
hrd_result <- scar_score(
  seg_file = "sequenza_segments.txt",
  reference = "grch38",
  seqz = TRUE  # 输入来自Sequenza
)

# 输出三个基因组瘢痕指标:
# HRD-LOH: 超过15Mb但不覆盖整条染色体的LOH区域数
# LST: 相邻segment间>10Mb的大规模状态转换数
# TAI: 延伸到端粒的等位不平衡区域数
print(hrd_result)
# HRD-LOH  LST  TAI  HRD_sum
# 12       15   10   37

# HRD评分 = LOH + LST + TAI
# 临床阈值: HRD >= 42 (Myriad myChoice CDx标准)
if (hrd_result$HRD_sum >= 42) {
  cat("HRD-positive: 可能对PARPi敏感\n")
} else {
  cat("HRD-negative\n")
}

# === 从ASCAT结果计算 ===
# ASCAT输出: allele-specific拷贝数
# 也可使用scarHRD处理ASCAT segments
hrd_ascat <- scar_score(
  seg_file = "ascat_segments.txt",
  reference = "grch38",
  seqz = FALSE  # 非Sequenza格式
)

library(maftools)

# === 读取MAF文件 ===
maf <- read.maf("tumor_mutations.maf")

# === DDR基因列表 ===
hr_genes <- c("BRCA1", "BRCA2", "PALB2", "RAD51B", "RAD51C", "RAD51D",
              "ATM", "ATR", "CHEK1", "CHEK2", "BARD1", "BRIP1", "FANCA",
              "FANCC", "FANCD2", "FANCE", "FANCF", "FANCG", "FANCL", "NBN")
mmr_genes <- c("MLH1", "MSH2", "MSH6", "PMS2", "EPCAM")
nhej_genes <- c("TP53BP1", "LIG4", "XRCC4", "XRCC5", "XRCC6")
ber_genes <- c("PARP1", "PARP2", "OGG1", "MUTYH", "XRCC1", "APEX1")

all_ddr_genes <- c(hr_genes, mmr_genes, nhej_genes, ber_genes)

# === 提取DDR基因突变 ===
ddr_mutations <- subsetMaf(maf, genes = all_ddr_genes)

# 瀑布图
oncoplot(maf, genes = all_ddr_genes[1:20])

# === BRCA突变分类 ===
brca_muts <- subsetMaf(maf, genes = c("BRCA1", "BRCA2"))
# 区分致病性突变类型
# Frameshift, Nonsense, Splice_Site → 功能丧失(LOF)
lof_types <- c("Frame_Shift_Del", "Frame_Shift_Ins", "Nonsense_Mutation", 
               "Splice_Site")
brca_lof <- brca_muts@data[Variant_Classification %in% lof_types]

# === 胚系vs体细胞BRCA突变 ===
# 胚系突变需要normal样本的VAF确认
# VAF ~ 50% in normal → 胚系
# VAF ~ 0% in normal → 体细胞

# === DDR通路缺陷评估 ===
# 综合判断: 突变 + 表达沉默 + HRD评分
ddr_status <- data.frame(
  Pathway = c("HR", "MMR", "NHEJ", "BER"),
  Gene_Mutation = c(
    any(ddr_mutations@data$Hugo_Symbol %in% hr_genes),
    any(ddr_mutations@data$Hugo_Symbol %in% mmr_genes),
    any(ddr_mutations@data$Hugo_Symbol %in% nhej_genes),
    any(ddr_mutations@data$Hugo_Symbol %in% ber_genes)
  ),
  HRD_positive = hrd_result$HRD_sum >= 42,
  Signature = NA  # 需要突变签名分析
)

# === MSIsensor2: 从BAM检测MSI ===
# 安装: conda install -c bioconda msisensor2

# 生成微卫星位点列表
msisensor2 scan -d reference.fa -o microsatellites.list

# 检测MSI状态
msisensor2 msi \
  -d microsatellites.list \
  -n normal.bam \
  -t tumor.bam \
  -o msi_output \
  -b 4

# 输出: msi_output
# Total_Number_of_Sites  Number_of_Somatic_Sites  %
# MSI score ≥ 3.5% → MSI-H; < 3.5% → MSS

# === 替代方法: MANTIS ===
# conda install -c bioconda mantis
mantis --bedfile microsatellites.bed \
  --genome reference.fa \
  -n normal.bam \
  -t tumor.bam \
  -o mantis_output.txt

# MANTIS score > 0.4 → MSI-H

# === 基于表达的MMR缺陷检测 ===
# MLH1甲基化导致的表达沉默是散发性MSI的主要原因

# 检查MMR基因表达
mmr_expression <- expression_matrix[mmr_genes, ]

# MLH1低表达 + MSI-H → 可能是MLH1启动子甲基化
# MLH1正常表达 + MSI-H → 可能是其他MMR基因突变(Lynch综合征)

library(GSVA)
library(msigdbr)

# === 基于基因集的DDR通路活性评估 ===
# 获取DDR相关基因集
ddr_genesets <- list(
  HR_repair = c("BRCA1", "BRCA2", "RAD51", "PALB2", "ATM", "ATR", "CHEK1", "CHEK2"),
  MMR = c("MLH1", "MSH2", "MSH6", "PMS2", "EXO1", "RFC1"),
  NHEJ = c("XRCC5", "XRCC6", "LIG4", "XRCC4", "DCLRE1C", "PRKDC"),
  BER = c("OGG1", "APEX1", "XRCC1", "PARP1", "PARP2", "LIG3"),
  NER = c("XPA", "XPB", "XPC", "XPD", "ERCC1", "DDB1", "DDB2"),
  FA_pathway = c("FANCA", "FANCB", "FANCC", "FANCD2", "FANCE", "FANCF", "FANCG")
)

# GSVA评估通路活性
gsva_result <- gsva(
  as.matrix(expression_tpm),
  ddr_genesets,
  method = "gsva",
  verbose = TRUE
)

# 热图展示各样本的DDR通路活性
library(pheatmap)
pheatmap(gsva_result,
         scale = "row",
         annotation_col = clinical_annotation,
         color = colorRampPalette(c("blue", "white", "red"))(100),
         main = "DDR Pathway Activity Scores")

# === 与PARP抑制剂响应关联 ===
# HR通路低活性 → 可能PARPi敏感
hr_low <- gsva_result["HR_repair", ] < median(gsva_result["HR_repair", ])
cat(sprintf("HR-low样本比例: %.1f%%\n", sum(hr_low)/length(hr_low)*100))

# === 完整DDR分析流程 ===
# 1. 变异检测
gatk Mutect2 -R ref.fa -I tumor.bam -I normal.bam -O somatic.vcf

# 2. 拷贝数分析(HRD需要)
sequenza-utils bam2seqz -n normal.bam --tumor tumor.bam --fasta ref.fa -o seqz.gz

# 3. HRD评分
Rscript -e 'library(scarHRD); scar_score("segments.txt", reference="grch38", seqz=TRUE)'

# 4. MSI检测
msisensor2 scan -d ref.fa -o ms.list
msisensor2 msi -d ms.list -n normal.bam -t tumor.bam -o msi_result

# 5. DDR基因突变提取
Rscript -e 'library(maftools); maf=read.maf("somatic.maf"); subsetMaf(maf, genes=ddr_genes)'