346_新生儿基因筛查数据分析¶

# ========== 靶向panel数据分析流程 ==========

# 步骤1：质控（FastQC + Trimmomatic）
fastqc sample_R1.fastq.gz sample_R2.fastq.gz \
    -o qc_reports/           # 质控报告目录

trimmomatic PE \
    sample_R1.fastq.gz sample_R2.fastq.gz \         # 输入双端reads
    sample_R1_clean.fastq.gz sample_R1_unpaired.fastq.gz \  # 输出
    sample_R2_clean.fastq.gz sample_R2_unpaired.fastq.gz \
    LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36  # 质控参数

# 步骤2：比对参考基因组
bwa mem \
    -t 8 \                     # 8个线程
    -R "@RG\tID:SAMPLE\tSM:SAMPLE\tPL:ILLUMINA" \  # 读组信息
    hg38.fa \                  # 参考基因组
    sample_R1_clean.fastq.gz sample_R2_clean.fastq.gz \
    | samtools sort -o sample.sorted.bam  # 直接排序输出

samtools index sample.sorted.bam  # 建立索引

# 步骤3：捕获效率统计
picard CollectHsMetrics \
    I=sample.sorted.bam \                # 输入BAM
    O=hs_metrics.txt \                   # 输出指标文件
    BAIT_INTERVALS=panel_baits.bed \     # 探针区域BED
    TARGET_INTERVALS=panel_targets.bed   # 目标区域BED

# 查看平均覆盖深度
grep "MEAN_TARGET_COVERAGE" hs_metrics.txt  # 目标区域平均覆盖深度

# ========== GATK变异检测（针对panel数据）==========

# 标记PCR重复
gatk MarkDuplicates \
    -I sample.sorted.bam \         # 输入BAM
    -O sample.dedup.bam \          # 去重后BAM
    -M dedup_metrics.txt           # 去重统计

gatk BuildBamIndex -I sample.dedup.bam  # 建立BAM索引

# BQSR碱基质量分数校正（可选，但推荐）
gatk BaseRecalibrator \
    -I sample.dedup.bam \
    -R hg38.fa \
    --known-sites dbsnp_146.hg38.vcf.gz \  # 已知变异数据库
    -O recal_data.table

gatk ApplyBQSR \
    -I sample.dedup.bam \
    -R hg38.fa \
    --bqsr-recal-file recal_data.table \
    -O sample.recal.bam  # 校正后BAM

# 变异检测
gatk HaplotypeCaller \
    -R hg38.fa \
    -I sample.recal.bam \
    -L panel_targets.bed \         # 仅分析目标区域
    -O sample.g.vcf.gz \           # GVCF格式输出
    -ERC GVCF                      # 输出GVCF用于联合分型

# ========== SMA SMN1拷贝数分析 ==========
# 使用SMNCopyNumberCaller工具分析SMN1/SMN2拷贝数

import subprocess   # 调用外部程序
import pandas as pd # 数据处理

def run_smn_analysis(bam_file, output_dir):
    """
    运行SMNCopyNumberCaller分析SMN1/SMN2拷贝数
    """
    # SMNCopyNumberCaller（推荐工具，专为SMA分析设计）
    cmd = [
        "smn_caller.py",
        "--bam", bam_file,                 # 输入BAM文件
        "--genome", "38",                   # 参考基因组版本
        "--prefix", "sample",              # 输出文件前缀
        "--outDir", output_dir,            # 输出目录
        "--depth", "10"                    # 最低覆盖深度
    ]

    result = subprocess.run(cmd, capture_output=True, text=True)

    if result.returncode != 0:
        print(f"分析失败: {result.stderr}")  # 输出错误信息
        return None

    # 读取拷贝数结果
    result_file = f"{output_dir}/sample_smn.tsv"
    df = pd.read_csv(result_file, sep='\t')

    # 解读结果
    smn1_copy = df['SMN1_CN'].values[0]   # SMN1拷贝数
    smn2_copy = df['SMN2_CN'].values[0]   # SMN2拷贝数

    print(f"SMN1拷贝数: {smn1_copy}")
    print(f"SMN2拷贝数: {smn2_copy}")

    # SMA诊断判断
    if smn1_copy == 0:
        print("警告：SMN1拷贝数为0，提示SMA高风险！")
        print(f"SMN2拷贝数={smn2_copy}，疾病严重程度预测参考")
    elif smn1_copy == 1:
        print("SMN1单拷贝，可能为SMA携带者")
    else:
        print("SMN1拷贝数正常")

    return df

# 使用示例
# result = run_smn_analysis("sample.bam", "smn_results/")

# ========== 批量样本NBS结果汇总 ==========
import pandas as pd    # 数据处理
import os              # 文件操作

def summarize_nbs_results(results_dir, disease_genes_file):
    """
    汇总多个样本的NBS分析结果
    """
    # 读取疾病-基因对应表
    disease_genes = pd.read_csv(disease_genes_file)  # 疾病对应基因表

    all_results = []  # 存储所有样本结果

    # 遍历所有样本结果文件
    for f in os.listdir(results_dir):
        if f.endswith("_variants.tsv"):  # 变异文件后缀
            sample_id = f.replace("_variants.tsv", "")  # 提取样本ID
            df = pd.read_csv(os.path.join(results_dir, f), sep='\t')

            # 筛选致病/可能致病变异
            pathogenic = df[df['ClinVar_sig'].str.contains(
                'Pathogenic|Likely_pathogenic', na=False
            )]

            if len(pathogenic) > 0:
                # 有P/LP变异，标记为需要随访
                for _, row in pathogenic.iterrows():
                    all_results.append({
                        'sample_id': sample_id,
                        'gene': row['Gene'],
                        'variant': f"{row['CHROM']}:{row['POS']} {row['REF']}>{row['ALT']}",
                        'classification': row['ClinVar_sig'],
                        'action': '需要遗传咨询'  # 标记处理意见
                    })

    # 输出汇总表
    summary = pd.DataFrame(all_results)
    summary.to_csv("nbs_positive_results.csv", index=False)
    print(f"共{len(summary)}条阳性结果，涉及{summary['sample_id'].nunique()}个样本")
    return summary