347_罕见病基因panel设计¶

# ========== Panel设计：提取目标区域BED文件 ==========

# 使用PanelApp数据库API获取panel基因列表（英国NHS推荐）
python3 << 'EOF'
import requests   # HTTP请求
import json       # JSON解析

# PanelApp API：获取特定panel的基因列表
# 示例：遗传性癫痫panel (panel ID: 402)
panel_id = 402
url = f"https://panelapp.genomicsengland.co.uk/api/v1/panels/{panel_id}/genes/"

response = requests.get(url, timeout=30)  # 发送请求
data = response.json()

# 提取基因名称（仅绿色=Definitive级别）
genes = []
for gene in data['results']:
    if gene['confidence_level'] == '3':  # 3=绿色=高证据
        genes.append(gene['gene_data']['gene_symbol'])  # 提取基因名

print(f"Panel基因数量: {len(genes)}")
print("基因列表:", genes[:10], "...")  # 显示前10个

# 保存基因列表
with open("panel_genes.txt", "w") as f:
    f.write("\n".join(genes))  # 每行一个基因名
EOF

# ========== 从基因列表生成目标区域BED ==========

# 使用UCSC Table Browser或BioMart获取基因外显子坐标
# 这里用biomaRt R包自动化获取

Rscript << 'EOF'
library(biomaRt)    # 基因组数据库接口
library(dplyr)      # 数据处理

# 读取Panel基因列表
genes <- readLines("panel_genes.txt")  # 读取基因名

# 连接Ensembl数据库（GRCh38）
mart <- useMart("ensembl", 
                dataset = "hsapiens_gene_ensembl",
                host = "https://grch38.ensembl.org")  # GRCh38版本

# 获取外显子坐标
exons <- getBM(
    attributes = c("chromosome_name", "exon_chrom_start", 
                   "exon_chrom_end", "hgnc_symbol"),  # 需要的字段
    filters = "hgnc_symbol",           # 按基因名过滤
    values = genes,                    # 目标基因列表
    mart = mart
)

# 添加±50bp捕获剪接位点
exons$start_padded <- exons$exon_chrom_start - 50  # 左侧延伸50bp
exons$end_padded <- exons$exon_chrom_end + 50      # 右侧延伸50bp

# 过滤标准染色体
exons <- exons[exons$chromosome_name %in% c(1:22, "X", "Y"), ]

# 格式化为BED格式（0-based坐标）
bed <- data.frame(
    chrom = paste0("chr", exons$chromosome_name),
    start = exons$start_padded - 1,   # BED格式0-based
    end = exons$end_padded,
    name = exons$hgnc_symbol
)

# 排序和去重
bed <- bed %>% 
    arrange(chrom, start) %>%   # 按染色体位置排序
    distinct()                   # 去除重复区域

# 保存BED文件
write.table(bed, "panel_targets.bed", 
            sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)

cat("BED文件生成完成，共", nrow(bed), "个区间\n")
EOF

# ========== Panel测序数据质控 ==========

# 计算Panel覆盖统计
mosdepth \
    --by panel_targets.bed \     # Panel目标区域
    --no-abbrev \                # 不简写报告
    --quantize 0:10:20:50: \    # 按深度分层统计（0x,10x,20x,50x）
    --threads 4 \               # 并行线程
    sample_coverage \           # 输出前缀
    sample.bam                  # 输入BAM

# 查看覆盖汇总
cat sample_coverage.mosdepth.summary.txt  # 均值、中位数覆盖深度

# 计算各深度覆盖率
python3 << 'EOF'
import pandas as pd   # 数据处理

# 读取mosdepth分位数文件
df = pd.read_csv("sample_coverage.quantized.bed.gz", 
                 sep='\t', names=['chr','start','end','depth_range'])

# 计算Panel各深度覆盖比例
total_bases = df['end'].sub(df['start']).sum()  # 总碱基数

for threshold in [10, 20, 50, 100]:
    # 筛选超过阈值深度的区域
    covered = df[df['depth_range'].apply(
        lambda x: int(x.split(':')[0]) >= threshold if ':' in x else False
    )]
    covered_bases = covered['end'].sub(covered['start']).sum()
    pct = covered_bases / total_bases * 100
    print(f">={threshold}x 覆盖率: {pct:.1f}%")
EOF

# ========== Panel CNV检测（ExomeDepth）==========
# ExomeDepth利用多样本对比检测CNV

# 这里展示R语言ExomeDepth基本流程
import subprocess   # 调用R脚本

r_script = """
library(ExomeDepth)   # CNV检测包

# 读取目标区域BED
exons <- read.table("panel_targets.bed", header=FALSE,
                    col.names=c("chromosome","start","end","name"))

# 从BAM文件统计各区域读数
# 需要先用samtools对每个样本生成计数
bam_files <- c("sample1.bam", "sample2.bam", "control1.bam")
counts <- getBamCounts(bed.frame = exons,
                       bam.files = bam_files,
                       min.mapq = 20)  # 最低比对质量

# 选择参考样本（正常对照）
my_choice <- select.reference.set(
    test.counts = as.vector(counts[, "sample1"]),     # 待分析样本
    reference.counts = counts[, c("control1", "control2")],  # 对照
    bin.length = exons$end - exons$start             # 区域长度
)

# 运行ExomeDepth CNV检测
mytest <- new("ExomeDepth",
    test = as.vector(counts[, "sample1"]),           # 测试样本计数
    reference = as.vector(my_choice$reference.count), # 参考计数
    formula = "cbind(test, reference) ~ 1")

mytest <- CallCNVs(x = mytest, transition.probability = 10^-4,
                   chromosome = exons$chromosome,
                   start = exons$start,
                   end = exons$end, name = exons$name)

# 输出CNV结果
write.csv(mytest@CNV.calls, "cnv_results.csv", row.names=FALSE)
cat("CNV检测完成，发现", nrow(mytest@CNV.calls), "个CNV\n")
"""

with open("run_exomedepth.R", "w") as f:
    f.write(r_script)  # 保存R脚本

subprocess.run(["Rscript", "run_exomedepth.R"])  # 运行R脚本