肿瘤克隆进化分析¶

# VAF → CCF 转换公式
# CCF = VAF × (purity × CN_total + (1-purity) × 2) / (purity × CN_mutant)
#
# purity: 肿瘤纯度
# CN_total: 突变位点的总拷贝数
# CN_mutant: 携带突变的等位基因拷贝数（通常为1）
#
# 例：纯度=0.8, 二倍体区域(CN=2), 杂合突变
# CCF = VAF × (0.8×2 + 0.2×2) / (0.8×1) = VAF × 2/0.8 = 2.5 × VAF
# 如果 VAF=0.4, 则 CCF=1.0（克隆突变）
# 如果 VAF=0.2, 则 CCF=0.5（亚克隆，50%肿瘤细胞携带）

# === 输入数据准备 ===

# 1. 体细胞突变检测（Mutect2）
gatk Mutect2 \
    -R reference.fasta \
    -I tumor.bam \
    -I normal.bam \
    -normal normal_sample_name \
    -O somatic.vcf.gz \
    --f1r2-tar-gz f1r2.tar.gz

# 过滤
gatk FilterMutectCalls \
    -R reference.fasta \
    -V somatic.vcf.gz \
    -O somatic.filtered.vcf.gz

# 2. 拷贝数分析（FACETS / Sequenza / ASCAT）
# FACETS示例
Rscript run_facets.R tumor.bam normal.bam

# 3. 肿瘤纯度估计（从FACETS/ASCAT获得）
# purity通常从拷贝数分析中同时估计

# 4. 准备PyClone输入
# 从VCF提取必要信息
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t[%AD]\n' \
    somatic.filtered.vcf.gz > mutations_for_pyclone.txt

# 准备PyClone-VI输入文件
import pandas as pd
import numpy as np

# 读取突变信息
mutations = pd.read_csv("mutations_for_pyclone.txt", sep="\t",
                         header=None, names=['chrom', 'pos', 'ref', 'alt', 'AD'])

# 解析AD字段（ref_count,alt_count）
mutations[['ref_counts', 'alt_counts']] = mutations['AD'].str.split(',', expand=True).astype(int)
mutations['total_counts'] = mutations['ref_counts'] + mutations['alt_counts']
mutations['vaf'] = mutations['alt_counts'] / mutations['total_counts']

# 读取拷贝数信息（从FACETS/Sequenza）
cn_segments = pd.read_csv("copy_number_segments.txt", sep="\t")

# 为每个突变匹配CN信息
def get_cn(row, cn_df):
    seg = cn_df[(cn_df['chrom'] == row['chrom']) &
                 (cn_df['start'] <= row['pos']) &
                 (cn_df['end'] >= row['pos'])]
    if len(seg) > 0:
        return seg.iloc[0]['major_cn'], seg.iloc[0]['minor_cn']
    return 1, 1

mutations[['major_cn', 'minor_cn']] = mutations.apply(
    lambda row: pd.Series(get_cn(row, cn_segments)), axis=1)
mutations['total_cn'] = mutations['major_cn'] + mutations['minor_cn']

# 过滤：只保留有足够深度的位点
mutations_filtered = mutations[mutations['total_counts'] >= 20].copy()

# 生成PyClone-VI输入
pyclone_input = pd.DataFrame({
    'mutation_id': mutations_filtered['chrom'] + ':' + mutations_filtered['pos'].astype(str),
    'sample_id': 'tumor_sample',
    'ref_counts': mutations_filtered['ref_counts'],
    'alt_counts': mutations_filtered['alt_counts'],
    'major_cn': mutations_filtered['major_cn'],
    'minor_cn': mutations_filtered['minor_cn'],
    'normal_cn': 2,
    'tumour_content': 0.75  # 肿瘤纯度
})
pyclone_input.to_csv("pyclone_input.tsv", sep="\t", index=False)

# === PyClone-VI 运行 ===
pip install pyclone-vi

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

# 参数说明：
# -c 40: 最大簇数（算法会自动选择实际簇数）
# -d beta-binomial: 数据分布模型（推荐）
# -r 10: 重启次数（取最优ELBO）

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

# PyClone-VI 结果解析与可视化
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# 读取结果
results = pd.read_csv("pyclone_clusters.tsv", sep="\t")

# 查看克隆结构
cluster_summary = results.groupby('cluster_id').agg(
    n_mutations=('mutation_id', 'count'),
    mean_ccf=('cellular_prevalence', 'mean'),
    std_ccf=('cellular_prevalence', 'std')
).sort_values('mean_ccf', ascending=False)

print("Clone structure:")
print(cluster_summary)

# 可视化：CCF分布
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# 左：CCF密度图
for cluster in results['cluster_id'].unique():
    subset = results[results['cluster_id'] == cluster]
    axes[0].hist(subset['cellular_prevalence'], alpha=0.5,
                 label=f"Cluster {cluster} (n={len(subset)})", bins=20)
axes[0].set_xlabel("Cellular Prevalence (CCF)")
axes[0].set_ylabel("Number of mutations")
axes[0].legend()

# 右：簇大小 vs 平均CCF
axes[1].scatter(cluster_summary['mean_ccf'], cluster_summary['n_mutations'],
                s=100, alpha=0.7)
for idx, row in cluster_summary.iterrows():
    axes[1].annotate(f"C{idx}", (row['mean_ccf'], row['n_mutations']))
axes[1].set_xlabel("Mean CCF")
axes[1].set_ylabel("Number of mutations")
plt.tight_layout()
plt.savefig("clone_structure.pdf")

# === SciClone 分析 ===
library(sciClone)

# 准备输入
# 1. VAF表（chr, pos, ref_reads, var_reads, vaf）
vafs <- read.table("sample_vafs.txt", header = TRUE)
# 列：chr, pos, ref_count, var_count, vaf

# 2. 拷贝数segments
cn_segs <- read.table("copy_number_segs.txt", header = TRUE)
# 列：chr, start, end, segment_mean (log2ratio)

# 3. 需排除的区域（LOH区域）
loh_regions <- read.table("loh_regions.txt", header = TRUE)

# 运行SciClone
sc <- sciClone(
  vafs = vafs,
  copyNumberCalls = cn_segs,
  sampleNames = "Tumor_Sample",
  regionsToExclude = loh_regions,
  minimumDepth = 20,
  maximumClusters = 10,
  clusterMethod = "bmm"  # beta mixture model
)

# 输出结果
writeClusterTable(sc, "sciclone_clusters.txt")

# 可视化
sc.plot1d(sc, "sciclone_1d.pdf")
sc.plot2d(sc, "sciclone_2d.pdf")  # 需要多区域采样

# 多区域采样分析（2个区域）
sc_multi <- sciClone(
  vafs = list(vafs_region1, vafs_region2),
  copyNumberCalls = list(cn_region1, cn_region2),
  sampleNames = c("Region1", "Region2"),
  minimumDepth = 20
)
sc.plot2d(sc_multi, "multi_region_2d.pdf")

# === PhyloWGS 运行 ===

# 准备SSM（Simple Somatic Mutations）文件
# 格式：id, gene, a(var_reads), d(total_reads), mu_r, mu_v
python parser/create_phylowgs_inputs.py \
    --vcf-type mutect \
    --tumor-sample TUMOR \
    --cnvs cnv_data.txt \
    --output-variants ssm_data.txt \
    --output-params params.json \
    somatic.filtered.vcf.gz

# 运行PhyloWGS
python multievolve.py \
    --num-chains 4 \
    --ssms ssm_data.txt \
    --cnvs cnv_data.txt \
    --params params.json \
    --output-dir phylowgs_output

# 解析结果
python write_results.py \
    phylowgs_output/trees.zip \
    phylowgs_output/results \
    --include-ssm-names

# PhyloWGS结果可视化
import json
import numpy as np

# 读取最优树
with open("phylowgs_output/results/best_tree.json") as f:
    tree = json.load(f)

# 打印克隆结构
for node_id, node_data in tree['populations'].items():
    print(f"Clone {node_id}:")
    print(f"  CCF: {node_data['cellular_prevalence']:.3f}")
    print(f"  Num SSMs: {len(node_data['ssms'])}")
    print(f"  Num CNVs: {len(node_data['cnvs'])}")
    print(f"  Children: {node_data.get('children', [])}")
    print()

# 使用ETE3绘制进化树
from ete3 import Tree, TreeStyle

def build_newick(tree_data, node_id="0"):
    children = tree_data['populations'][node_id].get('children', [])
    if not children:
        return f"Clone_{node_id}"
    child_strs = [build_newick(tree_data, str(c)) for c in children]
    return f"({','.join(child_strs)})Clone_{node_id}"

newick = build_newick(tree) + ";"
t = Tree(newick)
t.render("evolution_tree.pdf")

# === 多区域/纵向分析 ===

# 使用CITUP（Clonality Inference in Tumors Using Phylogenetics）
# 或 CliP（Clonality in Pairwise samples）

# 方法1：PyClone多样本联合分析
# 准备多样本输入（同一个病人的多个样本）
# pyclone_input_multisample.tsv 包含多个sample_id

# PyClone-VI 自动处理多样本
pyclone-vi fit \
    -i multisample_input.tsv \
    -o multisample_results.h5 \
    -c 40 \
    -d beta-binomial \
    -r 10

# 方法2：使用fishplot可视化克隆动态
library(fishplot)

# 定义时间点
timepoints <- c(0, 30, 100, 200)  # 天数

# 定义各克隆在各时间点的细胞比例
# 行=克隆, 列=时间点
frac.table <- matrix(
  c(100, 50, 10, 0,     # Clone A (对治疗敏感)
     30, 80, 90, 95,    # Clone B (耐药克隆)
     0,  10, 20, 5),    # Clone C (新出现)
  ncol = 4, byrow = TRUE
)

# 定义父子关系
parents <- c(0, 1, 1)  # Clone B和C都从Clone A衍生

# 创建fishplot
fish <- createFishObject(frac.table, parents, timepoints = timepoints,
                          clone.labels = c("Clone A", "Clone B", "Clone C"))
fish <- setCol(fish, c("#E41A1C", "#377EB8", "#4DAF4A"))
fishPlot(fish, shape = "spline", title.btm = "Clonal Evolution Under Treatment",
         vlines = c(30), vlab = c("Treatment Start"))

# === 克隆进化结果的临床解读 ===

# 1. 驱动突变与克隆归属
# 确定每个驱动突变属于哪个克隆
driver_genes <- c("TP53", "KRAS", "PIK3CA", "PTEN")
driver_mutations <- results[results$gene %in% driver_genes, ]
driver_clone_assignment <- merge(driver_mutations, cluster_results,
                                  by = "mutation_id")

# 2. 克隆的时间顺序推断
# 克隆型(CCF≈1)突变：早期/创始突变
# 亚克隆(CCF<1)突变：晚期/进展相关突变
early_mutations <- results[results$cellular_prevalence > 0.9, ]
late_mutations <- results[results$cellular_prevalence < 0.5, ]

# 3. 治疗相关克隆动态
# 比较治疗前后的CCF变化
pre_treatment <- read.csv("pre_treatment_ccf.csv")
post_treatment <- read.csv("post_treatment_ccf.csv")
ccf_change <- merge(pre_treatment, post_treatment, by = "cluster_id",
                     suffixes = c("_pre", "_post"))
ccf_change$delta_ccf <- ccf_change$ccf_post - ccf_change$ccf_pre
# delta_ccf > 0: 耐药克隆扩张
# delta_ccf < 0: 对治疗敏感

# 完整克隆进化分析流程
# 1. 变异检测
gatk Mutect2 -I tumor.bam -I normal.bam -normal NORMAL -O somatic.vcf.gz
# 2. 拷贝数
Rscript run_facets.R tumor.bam normal.bam
# 3. 准备输入
python prepare_pyclone_input.py somatic.vcf.gz cn_segments.txt 0.75
# 4. 克隆推断
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
# 5. 进化树
python phylowgs/multievolve.py --ssms ssm.txt --cnvs cnv.txt