肿瘤转移克隆分析¶

# === PyClone-VI 克隆结构推断 ===
import subprocess                      # 命令执行

# 1. 准备输入文件（TSV格式）
# 列：mutation_id, sample_id, ref_counts, alt_counts,
#     major_cn, minor_cn, normal_cn, tumour_content
input_columns = """mutation_id\tsample_id\tref_counts\talt_counts\tmajor_cn\tminor_cn\tnormal_cn\ttumour_content
mut_1\tprimary\t100\t50\t2\t1\t2\t0.7
mut_1\tmetastasis\t120\t45\t2\t1\t2\t0.6
mut_2\tprimary\t90\t30\t2\t0\t2\t0.7
mut_2\tmetastasis\t110\t5\t2\t0\t2\t0.6
"""

# 2. 运行 PyClone-VI
cmd = (
    "pyclone-vi fit "                  # 拟合模型
    "-i input.tsv "                    # 输入文件
    "-o pyclone_results.h5 "           # 输出文件
    "-c 40 "                           # 最大簇数
    "-d beta-binomial "                # 分布模型
    "-r 10"                            # 重启次数
)
# subprocess.run(cmd, shell=True)

# 3. 提取结果
cmd2 = (
    "pyclone-vi write-results-file "   # 写出结果
    "-i pyclone_results.h5 "           # 输入
    "-o pyclone_results.tsv"           # 输出TSV
)
# subprocess.run(cmd2, shell=True)

# === ClonEvol 克隆演化可视化 ===
library(clonevol)                      # 克隆演化分析包

# 1. 准备输入数据（PyClone结果整理）
clone_data <- data.frame(
  cluster = c(1, 1, 2, 2, 3, 3),      # 克隆簇编号
  sample = c("primary", "met",         # 样本来源
             "primary", "met",
             "primary", "met"),
  ccf = c(0.95, 0.90,                  # 克隆细胞比例
          0.45, 0.70,                   # 克隆2在转移灶中扩增！
          0.30, 0.05)                   # 克隆3在转移灶中缩小
)

# 2. 推断克隆树
tree <- infer.clonal.models(
  variants = clone_data,               # 克隆数据
  cluster.col.name = "cluster",        # 簇列名
  ccf.col.name = "ccf",              # CCF列名
  sample.names = c("primary", "met"), # 样本名
  cancer.initiation.model = "monoclonal",  # 单克隆起源
  subclonal.test = "bootstrap",        # 自助法检验
  num.boots = 1000                     # 自助次数
)

# 3. 可视化克隆演化（鱼图）
plot.clonal.models(
  tree,                                # 克隆树
  matched = NULL,                      # 匹配参数
  variants = clone_data,               # 数据
  box.plot = TRUE,                     # 箱线图模式
  out.format = "pdf",                  # 输出格式
  out.dir = "clonevol_output/"         # 输出目录
)

# 鱼图（fishplot）展示克隆随时间/空间的演变

# 原发灶 vs 转移灶突变比较分析
import pandas as pd                    # 数据处理
import matplotlib.pyplot as plt        # 绑图
from matplotlib_venn import venn2      # 韦恩图

# 读取两个样本的突变
primary = pd.read_csv("primary_mutations.vcf.gz",
                       sep="\t", comment="#",
                       usecols=[0,1,3,4],
                       names=["chr","pos","ref","alt"])
met = pd.read_csv("metastasis_mutations.vcf.gz",
                   sep="\t", comment="#",
                   usecols=[0,1,3,4],
                   names=["chr","pos","ref","alt"])

# 创建突变ID
primary["mut_id"] = primary.apply(
    lambda r: f"{r['chr']}:{r['pos']}:{r['ref']}>{r['alt']}", axis=1)
met["mut_id"] = met.apply(
    lambda r: f"{r['chr']}:{r['pos']}:{r['ref']}>{r['alt']}", axis=1)

# 集合比较
shared = set(primary["mut_id"]) & set(met["mut_id"])     # 共有突变（trunk）
primary_only = set(primary["mut_id"]) - set(met["mut_id"])  # 原发独有
met_only = set(met["mut_id"]) - set(primary["mut_id"])      # 转移独有

print(f"共有突变（trunk）: {len(shared)}")
print(f"原发灶独有: {len(primary_only)}")
print(f"转移灶独有: {len(met_only)}")

# 韦恩图可视化
plt.figure(figsize=(8, 6))
venn2([set(primary["mut_id"]), set(met["mut_id"])],
      set_labels=("原发灶", "转移灶"))  # 两组标签
plt.title("原发灶 vs 转移灶突变对比")
plt.savefig("primary_vs_met_venn.png", dpi=150)