149_肿瘤突变签名分析¶

# 96种SBS突变类型构成
substitutions <- c("C>A", "C>G", "C>T", "T>A", "T>C", "T>G")
flanking_5 <- c("A", "C", "G", "T")
flanking_3 <- c("A", "C", "G", "T")

# 生成96种三核苷酸上下文
contexts <- expand.grid(
  sub = substitutions,
  f5 = flanking_5,
  f3 = flanking_3
)
contexts$label <- paste0(contexts$f5, "[", contexts$sub, "]", contexts$f3)
cat("Total mutation types:", nrow(contexts), "\n")  # 96

library(MutationalPatterns)
library(BSgenome.Hsapiens.UCSC.hg38)

# 1. 加载VCF文件
vcf_files <- list.files("vcf_dir/", pattern = "*.vcf$", full.names = TRUE)
sample_names <- gsub(".vcf", "", basename(vcf_files))

# 读取VCF (GRanges格式)
grl <- read_vcfs_as_granges(
  vcf_files = vcf_files,
  sample_names = sample_names,
  genome = "hg38"
)

# 2. 获取突变类型和三核苷酸上下文
ref_genome <- "BSgenome.Hsapiens.UCSC.hg38"
muts <- get_mut_type(grl, type = "snv")

# 3. 构建96类突变计数矩阵
mut_matrix <- mut_matrix(vcf_list = grl, ref_genome = ref_genome)
dim(mut_matrix)  # 96 × N(样本数)

# 4. 可视化突变谱
plot_96_profile(mut_matrix[, 1:4])  # 显示前4个样本的突变谱

# 5. 突变谱热图
plot_contribution_heatmap(mut_matrix, cluster_samples = TRUE)

# 安装: pip install SigProfilerExtractor SigProfilerMatrixGenerator

# === Step 1: 生成突变矩阵 ===
from SigProfilerMatrixGenerator.scripts import SigProfilerMatrixGeneratorFunc as matGen

# 从VCF生成突变矩阵
matrices = matGen.SigProfilerMatrixGeneratorFunc(
    project="my_project",
    genome="GRCh38",
    vcfpath="path/to/vcf_folder/",  # VCF文件所在目录
    plot=True,
    exome=False  # True用于外显子组数据
)

# === Step 2: 提取突变签名 ===
from SigProfilerExtractor import sigpro as sig

sig.sigProfilerExtractor(
    input_type="matrix",
    output="sig_output/",
    input_data="path/to/my_project/output/SBS/my_project.SBS96.all",
    reference_genome="GRCh38",
    minimum_signatures=1,
    maximum_signatures=10,    # 搜索1-10个签名
    nmf_replicates=100,       # NMF重复次数
    cpu=-1,                   # 使用所有CPU
    gpu=False
)

# === Step 3: 查看结果 ===
# 结果目录结构:
# sig_output/
#   SBS96/
#     Suggested_Solution/
#       SBS96_De-Novo_Solution/
#         Signatures/          - 提取的签名
#         Activities/          - 各样本的签名贡献
#       COSMIC_SBS96_Decomposed_Solution/
#         Signatures/          - 分解到COSMIC签名
#         Activities/          - COSMIC签名贡献

library(MutationalPatterns)

# === De novo签名提取(NMF) ===
# 估计最优签名数
estimate <- estimate_rank(
  mut_matrix = mut_matrix,
  rank_range = 2:8,
  nrun = 50
)

# 可视化选择最优K
plot(estimate)  # 看cophenetic/silhouette指标

# 提取签名 (假设选择K=4)
nmf_res <- extract_signatures(
  mut_matrix = mut_matrix,
  rank = 4,
  nrun = 100
)

# 签名矩阵和贡献矩阵
signatures <- nmf_res$signatures
contributions <- nmf_res$contribution

# 可视化提取的签名
plot_96_profile(signatures)

# === 与COSMIC签名比对 ===
# 加载COSMIC签名
cosmic_sbs <- get_known_signatures(muttype = "snv", source = "COSMIC_v3.3")

# 计算余弦相似度
sim_matrix <- cos_sim_matrix(signatures, cosmic_sbs)
plot_cosine_heatmap(sim_matrix, cluster_rows = TRUE)

# === 签名归因(Fitting): 将COSMIC签名拟合到样本 ===
# 严格拟合
fit_res <- fit_to_signatures(
  mut_matrix = mut_matrix,
  signatures = cosmic_sbs
)

# 查看各签名贡献
contributions_cosmic <- fit_res$contribution
plot_contribution(contributions_cosmic, coord_flip = TRUE, mode = "absolute")
plot_contribution(contributions_cosmic, coord_flip = TRUE, mode = "relative")

# === 严格归因(减少过拟合) ===
# 使用strict拟合, 移除贡献极小的签名
strict_res <- fit_to_signatures_strict(
  mut_matrix = mut_matrix,
  signatures = cosmic_sbs,
  max_delta = 0.004  # 余弦相似度下降阈值
)
plot_contribution(strict_res$fit_res$contribution, mode = "relative")

# === DBS签名 ===
# 获取DBS突变矩阵
dbs_contexts <- get_dbs_context(grl)
dbs_matrix <- count_dbs_contexts(dbs_contexts)

# 可视化DBS谱
plot_dbs_contexts(dbs_matrix[, 1:4])

# 拟合COSMIC DBS签名
cosmic_dbs <- get_known_signatures(muttype = "dbs", source = "COSMIC_v3.3")
dbs_fit <- fit_to_signatures(dbs_matrix, cosmic_dbs)
plot_contribution(dbs_fit$contribution, mode = "relative")

# === ID签名 ===
# 获取Indel突变矩阵
indel_contexts <- get_indel_context(grl, ref_genome = ref_genome)
indel_matrix <- count_indel_contexts(indel_contexts)

# 可视化ID谱
plot_indel_contexts(indel_matrix[, 1:4])

# 拟合COSMIC ID签名
cosmic_id <- get_known_signatures(muttype = "indel", source = "COSMIC_v3.3")
id_fit <- fit_to_signatures(indel_matrix, cosmic_id)
plot_contribution(id_fit$contribution, mode = "relative")

# === HRD检测 ===
# SBS3 是 HRD 的标志性签名 (BRCA1/2突变相关)
# 结合ID6 (NHEJ修复留下的微同源介导删除)

# 计算SBS3贡献比例
sbs3_fraction <- contributions_cosmic["SBS3", ] / colSums(contributions_cosmic)

# HRD评分综合
# 通常结合3个基因组瘢痕指标：
# LOH (Loss of Heterozygosity)
# TAI (Telomeric Allelic Imbalance)
# LST (Large-Scale State Transitions)

# 使用scarHRD包
# BiocManager::install("sztup/scarHRD")
library(scarHRD)
# 输入: 从Sequenza或ASCAT得到的等位基因特异拷贝数
hrd_result <- scar_score(
  seg_file = "sample.sequenza.segments.txt",
  reference = "grch38",
  seqz = TRUE
)
# 输出: HRD-LOH, LST, TAI, 以及总HRD score

# === MMR缺陷检测 ===
# MMR缺陷相关签名: SBS6, SBS14, SBS15, SBS20, SBS21, SBS26, SBS44
# ID签名: ID1, ID2 (微卫星不稳定相关)
mmr_sigs <- c("SBS6", "SBS14", "SBS15", "SBS20", "SBS21", "SBS26", "SBS44")
mmr_contribution <- colSums(contributions_cosmic[mmr_sigs, , drop = FALSE])
mmr_fraction <- mmr_contribution / colSums(contributions_cosmic)

# 临床判定标准 (简化版)
classify_mmr <- function(mmr_frac, tmb) {
  # TMB: tumor mutation burden (mutations/Mb)
  if (tmb > 10 & mmr_frac > 0.15) return("dMMR/MSI-H")
  else return("pMMR/MSS")
}

# 签名贡献的时间变化
# 比较配对样本(如原发vs转移)
library(ggplot2)

# 假设有配对数据
primary_contrib <- contributions_cosmic[, grep("Primary", colnames(contributions_cosmic))]
metastasis_contrib <- contributions_cosmic[, grep("Metastasis", colnames(contributions_cosmic))]

# 计算各签名的变化
delta_sigs <- rowMeans(metastasis_contrib) - rowMeans(primary_contrib)
delta_df <- data.frame(
  Signature = names(delta_sigs),
  Delta = delta_sigs
)

ggplot(delta_df, aes(x = reorder(Signature, Delta), y = Delta)) +
  geom_bar(stat = "identity", fill = ifelse(delta_df$Delta > 0, "red", "blue")) +
  coord_flip() +
  labs(x = "COSMIC Signature", y = "Change (Metastasis - Primary)",
       title = "Mutational Signature Changes") +
  theme_minimal()

# === SigProfiler 安装 ===
pip install SigProfilerExtractor SigProfilerMatrixGenerator SigProfilerAssignment
# 安装参考基因组
python -c "from SigProfilerMatrixGenerator import install as genInstall; genInstall.install('GRCh38')"

# === 从VCF生成突变矩阵 ===
python -c "
from SigProfilerMatrixGenerator.scripts import SigProfilerMatrixGeneratorFunc as matGen
matGen.SigProfilerMatrixGeneratorFunc('project', 'GRCh38', 'vcf_dir/', plot=True)
"

# === 签名提取 ===
python -c "
from SigProfilerExtractor import sigpro as sig
sig.sigProfilerExtractor('matrix', 'output/', 'project/output/SBS/project.SBS96.all',
  reference_genome='GRCh38', minimum_signatures=1, maximum_signatures=15, nmf_replicates=100)
"

# === 签名归因(Assignment) ===
python -c "
from SigProfilerAssignment import Analyzer as Analyze
Analyze.cosmic_fit(samples='project/output/SBS/project.SBS96.all',
  output='assignment_output/', genome_build='GRCh38')
"

# === R环境: MutationalPatterns ===
Rscript -e '
BiocManager::install("MutationalPatterns")
BiocManager::install("BSgenome.Hsapiens.UCSC.hg38")
'