627_基因组结构变异检测¶

# === 安装三大SV检测工具 ===
conda install -c bioconda delly      # Delly安装
conda install -c bioconda manta      # Manta安装
conda install -c bioconda svaba      # SvABA安装

# === 辅助工具 ===
conda install -c bioconda survivor   # SV合并和比较
conda install -c bioconda samtools   # BAM处理
conda install -c bioconda bcftools   # VCF处理

# 验证安装
delly --version      # 查看Delly版本
configManta.py --version  # 查看Manta版本
svaba --version      # 查看SvABA版本

# === Delly SV检测 ===
# 输入：排序好的BAM文件 + 参考基因组

# 第1步：检测SV
delly call \
    -g reference.fa \              # 参考基因组
    -o delly_sv.bcf \              # 输出BCF格式
    -x delly_exclude.tsv \         # 排除区域（如着丝粒、端粒）
    tumor.bam \                    # 肿瘤样本BAM
    normal.bam                     # 配对正常样本BAM

# 第2步：过滤
delly filter \
    -f somatic \                   # 体细胞过滤模式
    -o delly_filtered.bcf \        # 过滤后的输出
    -s samples.tsv \               # 样本信息（样本名\ttumor/control）
    delly_sv.bcf                   # 输入BCF

# samples.tsv格式：
# tumor_sample   tumor
# normal_sample  control

# 第3步：BCF转VCF
bcftools view delly_filtered.bcf > delly_filtered.vcf  # 转换格式

# === 下载Delly排除区域文件 ===
# 这些区域的SV检测不可靠（重复序列、着丝粒等）
wget https://github.com/dellytools/delly/blob/main/excludeTemplates/human.hg38.excl.tsv

# 查看排除区域
head human.hg38.excl.tsv  # 格式：chr\tstart\tend

# === 种系SV检测（群体遗传学/遗传病）===
# 第1步：每个样本独立检测
for bam in samples/*.bam; do
    name=$(basename "$bam" .bam)
    delly call \
        -g reference.fa \
        -o ${name}_delly.bcf \
        -x human.hg38.excl.tsv \
        ${bam}
done

# 第2步：合并所有样本
delly merge \
    -o sites.bcf \                 # 合并的SV位点
    *_delly.bcf                    # 所有样本的BCF

# 第3步：重新基因分型
for bam in samples/*.bam; do
    name=$(basename "$bam" .bam)
    delly call \
        -g reference.fa \
        -v sites.bcf \             # 使用合并的位点
        -o ${name}_genotyped.bcf \
        ${bam}
done

# 第4步：合并基因分型结果
bcftools merge -m id -O b -o merged_genotyped.bcf *_genotyped.bcf

# 第5步：过滤
delly filter \
    -f germline \                  # 种系过滤模式
    -o germline_filtered.bcf \
    merged_genotyped.bcf

# === Manta：配置+运行 两步走 ===

# 第1步：配置分析
configManta.py \
    --tumorBam tumor.bam \         # 肿瘤BAM
    --normalBam normal.bam \       # 正常BAM
    --referenceFasta reference.fa \  # 参考基因组
    --runDir manta_output           # 输出目录

# 第2步：运行分析
python2 manta_output/runWorkflow.py \
    -j 16 \                        # 16线程
    -m local                       # 本地运行模式

# === 输出文件说明 ===
# manta_output/results/variants/
#   candidateSmallIndels.vcf.gz    — 小indel候选
#   candidateSV.vcf.gz             — SV候选（未过滤）
#   diploidSV.vcf.gz               — 种系SV（二倍体）
#   somaticSV.vcf.gz               — 体细胞SV（关注这个）

# 解压查看结果
zcat manta_output/results/variants/somaticSV.vcf.gz | grep -v "^#" | wc -l
# 查看体细胞SV数量

# === 仅肿瘤样本（无配对正常）===
configManta.py \
    --tumorBam tumor.bam \         # 只有肿瘤BAM
    --referenceFasta reference.fa \
    --runDir manta_tumor_only

python2 manta_tumor_only/runWorkflow.py -j 16 -m local

# 注意：没有配对正常样本时，假阳性会更多
# 需要额外过滤步骤（如与gnomAD-SV数据库比对去除常见多态性）

# === 种系SV检测（遗传病分析）===
configManta.py \
    --bam sample.bam \             # 单个样本（不加--tumorBam）
    --referenceFasta reference.fa \
    --runDir manta_germline

python2 manta_germline/runWorkflow.py -j 16 -m local

# 种系结果在 diploidSV.vcf.gz 中

# === SvABA：局部组装+比对 ===
svaba run \
    -t tumor.bam \                 # 肿瘤BAM
    -n normal.bam \                # 正常BAM
    -p 16 \                        # 16线程
    -G reference.fa \              # 参考基因组
    -a svaba_output \              # 输出前缀
    -D dbsnp.vcf.gz                # dbSNP数据库（过滤已知变异）

# === 输出文件 ===
# svaba_output.svaba.somatic.sv.vcf       — 体细胞SV
# svaba_output.svaba.somatic.indel.vcf    — 体细胞indel
# svaba_output.svaba.germline.sv.vcf      — 种系SV
# svaba_output.svaba.germline.indel.vcf   — 种系indel

# SvABA的优势：对小SV（50-300bp）检测灵敏度最高
# 因为它使用局部de novo组装重建断点序列

# === 用SURVIVOR合并多个工具的结果 ===

# 第1步：准备VCF文件列表
ls delly_filtered.vcf manta_somatic.vcf svaba_somatic.vcf > vcf_list.txt

# 第2步：合并
SURVIVOR merge \
    vcf_list.txt \                 # VCF文件列表
    1000 \                         # 断点距离阈值（bp）
    2 \                            # 至少2个工具支持
    1 \                            # 考虑SV类型
    1 \                            # 考虑SV方向（strand）
    0 \                            # 不考虑距离估计
    50 \                           # 最小SV长度（bp）
    merged_sv.vcf                  # 输出合并VCF

# 统计合并结果
echo "=== 合并后SV统计 ==="
grep -v "^#" merged_sv.vcf | awk '{print $5}' | sort | uniq -c | sort -rn
# 输出各类型SV数量

# 第3步：统计每个工具的贡献
SURVIVOR stats merged_sv.vcf -1 -1 -1 stats_output

# === 用AnnotSV注释SV ===
conda install -c bioconda annotsv  # 安装

AnnotSV \
    -SVinputFile merged_sv.vcf \   # 输入SV VCF
    -outputFile annotated_sv \     # 输出前缀
    -genomeBuild GRCh38 \          # 基因组版本
    -SVminSize 50                  # 最小SV大小

# AnnotSV注释信息包括：
# 受影响的基因、基因功能区域、已知SV数据库比对
# OMIM疾病关联、DGV频率、临床分级

# === R语言SV统计可视化 ===
library(ggplot2)  # 加载ggplot2
library(dplyr)    # 数据处理

# 读取VCF
library(VariantAnnotation)  # VCF处理
vcf <- readVcf("merged_sv.vcf", "hg38")  # 读取VCF

# 提取SV信息
sv_info <- data.frame(
    type = info(vcf)$SVTYPE,               # SV类型
    length = abs(info(vcf)$SVLEN),          # SV长度
    chr = as.character(seqnames(rowRanges(vcf))),  # 染色体
    support = info(vcf)$SUPP                # 工具支持数
)

# 柱状图：各类型SV数量
ggplot(sv_info, aes(x = type, fill = type)) +
    geom_bar() +
    labs(x = "SV Type", y = "Count",
         title = "Structural Variant Types") +
    theme_minimal() +
    scale_fill_brewer(palette = "Set2")

ggsave("sv_type_distribution.pdf", width = 8, height = 5)

# SV长度分布
ggplot(sv_info[sv_info$length > 0 & sv_info$length < 1e6, ],
       aes(x = log10(length), fill = type)) +
    geom_histogram(bins = 50, alpha = 0.7) +
    facet_wrap(~type, scales = "free_y") +
    labs(x = "log10(SV Length)", y = "Count",
         title = "SV Size Distribution") +
    theme_minimal()

ggsave("sv_size_distribution.pdf", width = 12, height = 8)