长读测序结构变异检测方法¶

# ===== 长读测序结构变异检测流程 =====

# 1. 比对（已完成的话跳过）
minimap2 -a -x map-hifi -t 16 \
    hg38.fa hifi_reads.fastq.gz | \
    samtools sort -@ 8 -o aligned.sorted.bam
samtools index aligned.sorted.bam

# ===== 方法1：pbsv（PacBio 官方） =====
# 第一步：发现 SV 信号（signatures）
pbsv discover --tandem-repeats hg38_trf.bed \
    aligned.sorted.bam signatures.svsig.gz

# 第二步：根据信号调用 SV
pbsv call --ccs hg38.fa signatures.svsig.gz \
    pbsv_output.vcf

# ===== 方法2：Sniffles2 =====
sniffles --input aligned.sorted.bam \
    --vcf sniffles_output.vcf \
    --reference hg38.fa \
    --threads 16 \
    --minsvlen 50 \           # 最小 SV 长度 50bp
    --minsupport auto         # 自动确定最小支持 reads

# ===== 方法3：cuteSV =====
mkdir -p cutesv_tmp
cuteSV aligned.sorted.bam hg38.fa \
    cutesv_output.vcf cutesv_tmp/ \
    --threads 16 \
    --min_size 50 \            # 最小 SV 50bp
    --min_support 3 \          # 至少 3 条 reads 支持
    --sample sample_name \
    --max_cluster_bias_INS 1000 \
    --diff_ratio_merging_INS 0.9 \
    --max_cluster_bias_DEL 1000 \
    --diff_ratio_merging_DEL 0.5

# ===== 合并多个 caller 的结果 =====
# 使用 SURVIVOR 合并
# 参数：最大距离 1000bp，至少 2 个 caller 支持，最小 50bp
ls pbsv_output.vcf sniffles_output.vcf cutesv_output.vcf > caller_list.txt
SURVIVOR merge caller_list.txt \
    1000 \   # 断点距离 ≤1000bp 视为同一个 SV
    2 \      # 至少 2 个 caller 检测到
    1 \      # 考虑 SV 类型
    1 \      # 考虑链方向
    0 \      # 不考虑大小差异比例
    50 \     # 最小 SV 大小 50bp
    merged_sv.vcf

# 统计合并后的 SV 数量和类型
bcftools query -f '%INFO/SVTYPE\n' merged_sv.vcf | \
    sort | uniq -c | sort -rn
# 输出示例：
# 15234 DEL    （缺失最多）
# 12456 INS    （插入次多）
# 1523 DUP     （重复）
# 456 INV      （倒位）
# 123 BND      （易位/断裂）

# ===== SV 注释 =====
# 使用 AnnotSV 注释 SV 的基因影响
AnnotSV -SVinputFile merged_sv.vcf \
    -outputFile annotated_sv \
    -genomeBuild GRCh38 \
    -SVminSize 50

# ===== Python: SV 统计与可视化 =====
import pandas as pd
import matplotlib.pyplot as plt

# 读取 VCF 文件中的 SV 信息
sv_data = []
with open("merged_sv.vcf") as f:
    for line in f:
        if line.startswith("#"):
            continue
        fields = line.strip().split("\t")
        info = dict(item.split("=") for item in fields[7].split(";")
                    if "=" in item)
        sv_data.append({
            "chr": fields[0],
            "pos": int(fields[1]),
            "type": info.get("SVTYPE", "UNK"),
            "length": abs(int(info.get("SVLEN", "0").split(",")[0]))
        })

df = pd.DataFrame(sv_data)

# SV 长度分布（按类型）
fig, ax = plt.subplots(figsize=(10, 6))
for sv_type in ["DEL", "INS", "DUP", "INV"]:
    subset = df[df["type"] == sv_type]
    ax.hist(subset["length"].clip(upper=10000),  # 截断到 10kb
            bins=50, alpha=0.5, label=sv_type)
ax.set_xlabel("SV 长度 (bp)")
ax.set_ylabel("数量")
ax.set_title("结构变异长度分布")
ax.legend()
plt.savefig("sv_length_distribution.png", dpi=150)