循环肿瘤细胞 CTC 分析¶

# === CTC 单细胞 CNV 分析 ===
import numpy as np                     # 数值计算
import pandas as pd                    # 数据处理
import matplotlib.pyplot as plt        # 绑图

def ctc_cnv_analysis(bam_path, bin_size=500000):
    """
    CTC 单细胞 CNV 分析（简化版）
    基于固定窗口的读段计数方法
    """
    import pysam                       # BAM处理

    bam = pysam.AlignmentFile(bam_path, "rb")

    # 按固定窗口统计读段数
    bin_counts = {}                    # 存储每个bin的计数
    for read in bam:                   # 遍历读段
        if read.mapping_quality >= 20 and not read.is_duplicate:
            chrom = read.reference_name  # 染色体
            pos = read.reference_start   # 位置
            bin_id = (chrom, pos // bin_size * bin_size)  # bin标识
            bin_counts[bin_id] = bin_counts.get(bin_id, 0) + 1

    # 转为 DataFrame
    df = pd.DataFrame([
        {"chrom": k[0], "start": k[1], "count": v}
        for k, v in bin_counts.items()
    ])

    # GC 校正（简化版）
    median_count = df["count"].median()  # 中位数
    df["ratio"] = df["count"] / median_count  # 比值
    df["log2"] = np.log2(df["ratio"] + 0.01)  # log2变换

    return df

# 可视化 CTC 全基因组 CNV
def plot_ctc_cnv(cnv_df, title="CTC CNV Profile"):
    """绘制 CTC 单细胞 CNV 图"""
    # 按染色体排序
    chrom_order = [f"chr{i}" for i in range(1, 23)] + ["chrX"]
    cnv_df["chrom_idx"] = cnv_df["chrom"].map(
        {c: i for i, c in enumerate(chrom_order)})
    cnv_df = cnv_df.dropna(subset=["chrom_idx"])
    cnv_df = cnv_df.sort_values(["chrom_idx", "start"])

    # 绘图
    fig, ax = plt.subplots(figsize=(16, 4))
    x_pos = range(len(cnv_df))         # X坐标
    colors = ["steelblue" if i % 2 == 0 else "coral"
              for i in cnv_df["chrom_idx"].astype(int)]

    ax.scatter(x_pos, cnv_df["log2"],  # 散点图
               c=colors, s=1, alpha=0.5)
    ax.axhline(y=0, color="black", linewidth=0.5)     # 基线
    ax.axhline(y=0.58, color="red", linewidth=0.5, linestyle="--")  # 扩增
    ax.axhline(y=-1, color="blue", linewidth=0.5, linestyle="--")   # 缺失
    ax.set_ylabel("log2 ratio")
    ax.set_title(title)
    ax.set_ylim(-3, 3)                 # Y轴范围
    plt.tight_layout()
    plt.savefig("ctc_cnv.png", dpi=150)

# === CTC 单细胞 WGS 分析流程 ===

# 1. 单细胞扩增后的 FASTQ 质控
fastp -i CTC_R1.fq.gz -I CTC_R2.fq.gz \  # 输入
  -o clean_R1.fq.gz -O clean_R2.fq.gz \   # 输出
  --detect_adapter_for_pe                   # 自动检测接头

# 2. 比对到参考基因组
bwa mem -t 8 \                         # 8线程
  hg38.fa \                            # 参考基因组
  clean_R1.fq.gz clean_R2.fq.gz | \   # 输入
  samtools sort -o ctc.bam             # 排序输出

# 3. 去除 PCR 重复（WGA 产生大量重复）
picard MarkDuplicates \
  I=ctc.bam \                          # 输入
  O=ctc_dedup.bam \                    # 输出
  M=dup_metrics.txt \                  # 重复统计
  REMOVE_DUPLICATES=true               # 删除重复

# 4. 检查 WGA 覆盖均匀性
picard CollectWgsMetrics \
  I=ctc_dedup.bam \                    # 输入
  O=wgs_metrics.txt \                  # 输出
  R=hg38.fa                            # 参考基因组

# 5. 使用 Ginkgo 或 SCOPE 做单细胞 CNV
# Ginkgo 是专门为单细胞 CNV 设计的工具
# 在线版：http://qb.cshl.edu/ginkgo

# 6. 使用 CopywriteR 做 CNV 分析
# （R 脚本）
cat > run_cnv.R << 'RSCRIPT'
library(CopywriteR)                    # 单细胞CNV工具

# 运行 CopywriteR
CopywriteR(
  sample.control = data.frame(
    samples = "ctc_dedup.bam",         # CTC BAM
    controls = "bulk_normal.bam"       # 正常对照BAM
  ),
  destination.folder = "cnv_output/",  # 输出目录
  reference.folder = "ref_folder/",    # 参考文件夹
  bp.param = BiocParallel::MulticoreParam(workers = 4)
)
RSCRIPT