肿瘤纯度与倍性估计¶

# === 使用 Sequenza 估计纯度和倍性 ===

# 1. 从 BAM 提取数据（需要配对 Tumor/Normal）
sequenza-utils bam2seqz \
  -n normal.bam \                      # 正常样本BAM
  -t tumor.bam \                       # 肿瘤样本BAM
  --fasta hg38.fa \                    # 参考基因组
  -gc gc50Base.wig.gz \                # GC含量文件
  -o sample.seqz.gz                    # 输出seqz文件

# 2. 分箱（减少数据量）
sequenza-utils seqz_binning \
  --seqz sample.seqz.gz \             # 输入seqz
  -w 50 \                              # 窗口50bp
  -o sample_bin.seqz.gz               # 输出分箱结果

# === 使用 FACETS（R语言） ===
# 3. 从 BAM 生成 pileup（使用 snp-pileup）
snp-pileup \
  -g \                                 # 仅保留杂合SNP
  -q15 -Q20 \                         # 最小比对和碱基质量
  dbsnp_common.vcf.gz \               # 已知SNP数据库
  pileup.csv.gz \                      # 输出pileup
  normal.bam \                         # 正常BAM
  tumor.bam                            # 肿瘤BAM

# === Sequenza R 分析（接上面的步骤2） ===
library(sequenza)                      # 加载sequenza包

# 读取分箱后的seqz数据
data <- sequenza.extract(              # 提取数据
  "sample_bin.seqz.gz"                 # 输入分箱文件
)

# 拟合模型（搜索最优纯度和倍性组合）
fit <- sequenza.fit(data)              # 模型拟合

# 查看结果
print(fit$purity)                      # 最优纯度估计
print(fit$ploidy)                      # 最优倍性估计

# 生成完整报告
sequenza.results(                      # 生成结果
  sequenza.extract = data,             # 数据
  cp.table = fit,                      # 拟合结果
  sample.id = "tumor_sample",          # 样本ID
  out.dir = "sequenza_results/"        # 输出目录
)

# === FACETS R 分析 ===
library(facets)                        # 加载facets包

# 读取 pileup 数据
rcmat <- readSnpMatrix("pileup.csv.gz")  # 读取SNP矩阵

# 预处理（GC校正+分段）
xx <- preProcSample(rcmat,             # 输入数据
                    gbuild = "hg38")   # 基因组版本

# 运行 FACETS
oo <- procSample(xx,                   # 处理样本
                 cval = 300)           # 分段灵敏度（越大越不灵敏）

# 拟合纯度和倍性
fit2 <- emcncf(oo)                     # EM算法拟合

cat("纯度:", fit2$purity, "\n")        # 打印纯度
cat("倍性:", fit2$ploidy, "\n")        # 打印倍性

# 可视化
plotSample(x = oo, emfit = fit2)       # 绘制全基因组CNV图