350_结构变异与遗传病¶

# ========== Manta：短读测序SV检测 ==========

# Manta是GATK生态外最常用的SV检测工具
# 安装：conda install -c bioconda manta

# 步骤1：配置分析（需要BAM索引）
configManta.py \
    --bam proband.bam \              # 先证者BAM文件
    --bam father.bam \               # 父亲BAM（三人组模式）
    --bam mother.bam \               # 母亲BAM
    --referenceFasta hg38.fa \       # 参考基因组
    --exome \                        # WES/Exome模式（调整参数）
    --runDir manta_results/          # 输出目录

# 步骤2：运行分析
./manta_results/runWorkflow.py \
    -m local \                       # 本地运行模式
    -j 16                            # 并行线程数

echo "Manta SV检测完成"

# 步骤3：查看结果
ls manta_results/results/variants/   # 列出结果文件
# 主要结果文件：
# candidateSV.vcf.gz    - 所有候选SV（未过滤）
# diploidSV.vcf.gz      - 经过滤的最终SV
# candidateSmallIndels.vcf.gz - 用于传给GATK的小Indel

# ========== CNVnator：读数深度法CNV检测 ==========

# 安装：conda install -c bioconda cnvnator

# 步骤1：提取读数历史（以100bp窗口为例）
cnvnator \
    -root sample.root \              # 输出ROOT文件
    -genome hg38 \                   # 参考基因组版本
    -tree sample.bam                 # 输入BAM文件

# 步骤2：计算覆盖深度直方图
cnvnator \
    -root sample.root \
    -his 100 \                       # 100bp bin大小
    -d /path/to/chromosomes/         # 各染色体FASTA目录

# 步骤3：统计分析
cnvnator \
    -root sample.root \
    -stat 100                        # 计算统计量

# 步骤4：分割检测（关键步骤）
cnvnator \
    -root sample.root \
    -partition 100                   # 分割算法

# 步骤5：CNV检出
cnvnator \
    -root sample.root \
    -call 100 \                      # 检出CNV
    > sample_cnv_calls.txt           # 保存结果

# 查看结果格式（CNV类型、位置、大小、q-value）
head -20 sample_cnv_calls.txt       # 查看前20个CNV

# ========== 解析和过滤SV VCF文件 ==========
import pysam    # VCF/BAM文件处理
import pandas as pd   # 数据处理

def parse_sv_vcf(sv_vcf_file, min_size=500, max_size=10_000_000):
    """
    解析SV VCF文件，提取结构变异信息
    sv_vcf_file: Manta/GATK SV输出的VCF文件
    min_size: 最小SV大小（bp）
    max_size: 最大SV大小（bp）
    """
    variants = []  # 存储变异信息

    # 打开VCF文件
    vcf = pysam.VariantFile(sv_vcf_file)

    for rec in vcf:
        # 获取SV类型
        sv_type = rec.info.get('SVTYPE', 'UNKNOWN')  # DEL/DUP/INV/BND

        # 获取SV大小
        sv_len = abs(rec.info.get('SVLEN', 0))
        if isinstance(sv_len, tuple):
            sv_len = abs(sv_len[0])  # SVLEN有时是元组

        # 过滤大小不合适的SV
        if sv_len < min_size or sv_len > max_size:
            continue  # 跳过太小或太大的SV

        # 过滤未通过质控的SV
        if list(rec.filter)[0] != 'PASS':
            continue  # 只保留PASS状态

        # 获取等位基因频率（通过INFO字段）
        pr = rec.samples[0].get('PR', (0, 0))  # paired reads支持
        sr = rec.samples[0].get('SR', (0, 0))  # split reads支持

        variants.append({
            'chrom': rec.chrom,           # 染色体
            'start': rec.start,           # 起始位置
            'end': rec.stop,              # 终止位置
            'sv_type': sv_type,           # SV类型
            'sv_len': sv_len,             # SV大小
            'PR_ref': pr[0] if pr else 0, # 参考等位基因读数
            'PR_alt': pr[1] if pr else 0, # 变异等位基因读数
        })

    return pd.DataFrame(variants)

# 使用示例
# df_sv = parse_sv_vcf("diploidSV.vcf.gz", min_size=100)
# print(f"过滤后SV数量: {len(df_sv)}")
# print(df_sv.groupby('sv_type').size())  # 按类型统计

# ========== CMA数据分析：PennCNV后处理 ==========
library(dplyr)    # 数据处理
library(GenomicRanges)   # 基因组区间操作

# 读取PennCNV输出的CNV list
cnv_calls <- read.table("sample.cnv", header=FALSE, 
    col.names=c("region","numSNP","length","state","file","startSNP","endSNP","conf"))

# 解析染色体和位置（格式：chr1:123-456）
cnv_calls <- cnv_calls %>%
    mutate(
        chrom = sub(":(.*)", "", region),   # 提取染色体号
        start = as.numeric(sub(".*:(\\d+)-.*", "\\1", region)),  # 提取起始
        end = as.numeric(sub(".*-(\\d+)$", "\\1", region)),       # 提取终止
        cn_type = ifelse(state %in% c(1,2), "DEL", "DUP")         # CNV类型
    ) %>%
    filter(length >= 100000)  # 过滤<100kb的CNV（减少假阳性）

cat("保留CNV数量:", nrow(cnv_calls), "\n")  # 输出保留数量

# 与数据库比对（DGV正常人群CNV数据库）
# 这里使用GRanges进行区间重叠分析
cnv_gr <- GRanges(
    seqnames = cnv_calls$chrom,
    ranges = IRanges(start=cnv_calls$start, end=cnv_calls$end)
)
cat("构建GRanges对象完成，共", length(cnv_gr), "个CNV\n")