164_肿瘤异质性定量¶

# 基本概念：从VAF理解克隆结构
library(ggplot2)

# 模拟一个肿瘤的VAF分布
set.seed(42)
# 克隆性突变 (CCF=1.0, 在纯度80%肿瘤中VAF≈0.4)
clonal_vaf <- rnorm(100, mean = 0.40, sd = 0.05)
# 亚克隆1 (CCF=0.5, VAF≈0.20)
subclone1_vaf <- rnorm(50, mean = 0.20, sd = 0.04)
# 亚克隆2 (CCF=0.3, VAF≈0.12)
subclone2_vaf <- rnorm(30, mean = 0.12, sd = 0.03)

all_vaf <- c(clonal_vaf, subclone1_vaf, subclone2_vaf)

# VAF密度图
df <- data.frame(VAF = all_vaf)
ggplot(df, aes(x = VAF)) +
  geom_density(fill = "steelblue", alpha = 0.5) +
  geom_vline(xintercept = c(0.40, 0.20, 0.12), linetype = "dashed", color = "red") +
  annotate("text", x = c(0.40, 0.20, 0.12), y = c(3, 2, 1.5),
           label = c("Clonal", "Subclone1", "Subclone2"), color = "red") +
  labs(x = "Variant Allele Frequency (VAF)", y = "Density",
       title = "VAF Distribution Reflecting Clonal Structure") +
  theme_minimal()

# === MATH评分计算 ===
calculate_MATH <- function(vaf_vector) {
  # MATH = 100 × MAD / median
  # MAD = median absolute deviation
  med <- median(vaf_vector)
  mad_val <- median(abs(vaf_vector - med))
  math_score <- 100 * mad_val / med
  return(math_score)
}

# 从VCF/MAF文件读取数据
library(maftools)
maf <- read.maf("tumor.maf")

# 提取VAF
vaf_data <- subsetMaf(maf)
vaf_values <- vaf_data$t_alt_count / (vaf_data$t_alt_count + vaf_data$t_ref_count)

# 计算MATH
math_score <- calculate_MATH(vaf_values)
cat(sprintf("MATH Score: %.2f\n", math_score))

# MATH评分解读:
# < 20: 低异质性(均匀肿瘤)
# 20-40: 中等异质性
# > 40: 高异质性

# === 使用maftools直接计算 ===
math_result <- math(maf, vafCol = "VAF")
print(math_result)

# 多样本比较
maf_list <- list(tumor1 = maf1, tumor2 = maf2, tumor3 = maf3)
math_scores <- sapply(maf_list, function(m) {
  math(m, vafCol = "VAF")$MATH
})

# === 基于克隆结构的多样性指数 ===

# 假设已通过PyClone等工具获得克隆cluster及其细胞比例
clonal_fractions <- c(0.50, 0.30, 0.15, 0.05)  # 4个克隆的细胞比例

# Shannon熵 (信息熵)
shannon_entropy <- function(fractions) {
  fractions <- fractions[fractions > 0]
  -sum(fractions * log2(fractions))
}

# Simpson指数 (1 - sum(pi^2))
simpson_index <- function(fractions) {
  1 - sum(fractions^2)
}

# Richness (克隆数)
richness <- length(clonal_fractions)

# Evenness (均匀度 = H / log2(S))
evenness <- function(fractions) {
  H <- shannon_entropy(fractions)
  S <- length(fractions)
  H / log2(S)
}

cat(sprintf("Shannon Entropy: %.3f bits\n", shannon_entropy(clonal_fractions)))
cat(sprintf("Simpson Index: %.3f\n", simpson_index(clonal_fractions)))
cat(sprintf("Richness: %d clones\n", richness))
cat(sprintf("Evenness: %.3f\n", evenness(clonal_fractions)))

# === 亚克隆比例 ===
# 定义: subclonal mutations / total mutations
# 通常用VAF阈值或PyClone聚类结果

calculate_subclonal_fraction <- function(vaf, purity = 0.8, ploidy = 2) {
  # 在给定纯度和倍性下，克隆性突变的期望VAF
  expected_clonal_vaf <- purity / (purity * ploidy + (1 - purity) * 2)

  # 低于期望VAF*0.7的突变视为亚克隆
  threshold <- expected_clonal_vaf * 0.7
  subclonal <- sum(vaf < threshold) / length(vaf)
  return(subclonal)
}

sub_frac <- calculate_subclonal_fraction(vaf_values, purity = 0.8)
cat(sprintf("亚克隆突变比例: %.1f%%\n", sub_frac * 100))

# === 安装PyClone-VI ===
pip install pyclone-vi

# === 准备输入文件 ===
# TSV格式: mutation_id, sample_id, ref_counts, alt_counts, 
#           major_cn, minor_cn, normal_cn, tumour_content

import pandas as pd

# 准备PyClone-VI输入
mutations = pd.DataFrame({
    'mutation_id': [f"mut_{i}" for i in range(200)],
    'sample_id': ['tumor1'] * 200,
    'ref_counts': ref_counts,     # 参考等位基因读数
    'alt_counts': alt_counts,     # 变异等位基因读数
    'major_cn': major_copy_number,  # 主要拷贝数
    'minor_cn': minor_copy_number,  # 次要拷贝数
    'normal_cn': [2] * 200,       # 正常拷贝数
    'tumour_content': [0.8] * 200  # 肿瘤纯度
})
mutations.to_csv("pyclone_input.tsv", sep='\t', index=False)

# === 运行PyClone-VI ===
pyclone-vi fit \
  -i pyclone_input.tsv \
  -o pyclone_results.h5 \
  -c 40 \
  -d beta-binomial \
  -r 10

# 参数:
# -c 40: 最大cluster数
# -d beta-binomial: 似然模型(推荐)
# -r 10: 随机重启次数

# 提取结果
pyclone-vi write-results-file \
  -i pyclone_results.h5 \
  -o pyclone_clusters.tsv

# === 多样本联合分析 ===
# 多区域测序或配对原发-转移
# 输入中不同sample_id对应不同区域
# PyClone-VI会联合推断克隆结构

# 结果可视化
library(ggplot2)
clusters <- read.delim("pyclone_clusters.tsv")

# 克隆频率分布
ggplot(clusters, aes(x = factor(cluster_id), y = cellular_prevalence)) +
  geom_boxplot(fill = "skyblue") +
  labs(x = "Clone Cluster", y = "Cellular Prevalence (CCF)",
       title = "Clonal Structure (PyClone-VI)") +
  theme_minimal()

# 鱼骨图 (多时间点/区域)
# 使用fishplot或ClonEvol包

# === 多区域测序分析 ===
# 从同一肿瘤的多个区域取样测序，可以直接观察空间异质性

# 读取多区域MAF数据
regions <- c("Region_A", "Region_B", "Region_C", "Region_D")

# 突变分类: ubiquitous(所有区域共有) vs private(某区域特有)
classify_mutations <- function(mutation_matrix) {
  # mutation_matrix: 突变×区域的presence/absence矩阵
  n_regions <- ncol(mutation_matrix)

  present_count <- rowSums(mutation_matrix > 0)

  classification <- case_when(
    present_count == n_regions ~ "Ubiquitous (Clonal)",
    present_count == 1 ~ "Private",
    TRUE ~ "Shared (Subclonal)"
  )
  return(classification)
}

# ITH指数: 1 - (ubiquitous mutations / total mutations)
ith_index <- 1 - sum(classification == "Ubiquitous (Clonal)") / length(classification)
cat(sprintf("ITH Index: %.3f\n", ith_index))

# === Jaccard距离矩阵 ===
# 区域间的突变组成差异
library(vegan)
dist_matrix <- vegdist(t(mutation_matrix > 0), method = "jaccard")
hclust_result <- hclust(dist_matrix, method = "ward.D2")
plot(hclust_result, main = "Tumor Regional Heterogeneity")

# === ITH分析完整流程 ===
# 1. 变异检测
gatk Mutect2 -R ref.fa -I tumor.bam -I normal.bam -O mutations.vcf
# 2. 拷贝数分析(PyClone需要)
sequenza-utils bam2seqz --fasta ref.fa -n normal.bam --tumor tumor.bam -o seqz.gz
# 3. PyClone-VI克隆结构推断
pyclone-vi fit -i input.tsv -o results.h5 -c 40 -d beta-binomial -r 10
pyclone-vi write-results-file -i results.h5 -o clusters.tsv
# 4. MATH评分
Rscript -e 'library(maftools); maf=read.maf("tumor.maf"); math(maf)'