Apr 15, 2026
Variants Calling

让AI整理了一下用于Variants Calling的代码，以后会是Vibe Coding的天下了可以确定自己在编程是确实没有天赋，但似乎也能发现技术和能力的分布不均

#!/bin/bash
# ============================================================
# Variants Calling Pipeline for Population Genomics
# ============================================================
# 功能: fastq到VCF的变异检测
# 作者: the_sun
# 日期: 2026-04-15
# 
# 流程步骤:
#   Step 1 - 质量过滤 (fastp)
#   Step 2 - 比对与标记重复 (bwa + GATK MarkDuplicates)
#   Step 3 - GVCF生成 (GATK HaplotypeCaller)
#   Step 4 - GenomicsDB构建 (GATK GenomicsDBImport)
#   Step 5 - 联合基因分型与合并 (GATK GenotypeGVCFs + MergeVcfs)
#
# 用法:
#   ./Variants_calling.sh              # 执行全部步骤
#   ./Variants_calling.sh 2 4          # 执行步骤2到4
#   ./Variants_calling.sh -s 3         # 仅执行步骤3
#   ./Variants_calling.sh -h           # 显示帮助信息
# ============================================================

set -e  # 遇错即停
set -u  # 使用未定义变量时报错
set -o pipefail  # 管道命令任一失败即报错

# ==================== 用户配置区 ====================
# 目录配置 (请根据实际情况修改)
INPUT_DIR="/ds3200_1/users_root/sunxiangqian/Analysis_test/01_SNP_Calling/sequence"
OUTPUT_DIR="/ds3200_1/users_root/sunxiangqian/Analysis_test/01_SNP_Calling"
REF="/ds3200_1/users_root/sunxiangqian/Analysis_test/01_SNP_Calling/reference/mira.fasta"

# 并行参数
FASTP_THREADS=4                        # fastp线程数
BWA_THREADS=4                          # bwa线程数
SAMTOOLS_THREADS=4                     # samtools线程数
GATK_THREADS=4                         # GATK线程数
PARALLEL_JOBS=2                        # GNU parallel并行任务数

# 质量过滤参数
MIN_LENGTH=50                          # 最小序列长度
QUALITY_CUTOFF=30                      # 质量 trimming 阈值 (Q值)

# 染色体列表 (请根据参考基因组修改)
CHROMOSOMES=("Pmi01" "Pmi02" "Pmi03" "Pmi04" "Pmi05" "Pmi06" "Pmi07" "Pmi08")

# Java内存设置 (根据服务器内存调整)
GATK_XMS="4g"                          # GATK初始内存
GATK_XMX="60g"                         # GATK最大内存
# ==================== 配置结束 ====================

# ==================== 全局变量 ====================
FILTER_DIR="${OUTPUT_DIR}/01_Filter"
BAM_DIR="${OUTPUT_DIR}/02_BAM"
GVCF_DIR="${OUTPUT_DIR}/03_GVCF"
DB_DIR="${OUTPUT_DIR}/04_DB"

# 颜色输出
RED='\033[0;31m'
GREEN='\033[0;32m'
YELLOW='\033[1;33m'
BLUE='\033[0;34m'
NC='\033[0m' # No Color

# ==================== 工具函数 ====================

log_info() {
    echo -e "${GREEN}[INFO]${NC} $(date '+%Y-%m-%d %H:%M:%S') $*"
}

log_error() {
    echo -e "${RED}[ERROR]${NC} $(date '+%Y-%m-%d %H:%M:%S') $*" >&2
}

log_warn() {
    echo -e "${YELLOW}[WARN]${NC} $(date '+%Y-%m-%d %H:%M:%S') $*"
}

log_step() {
    echo -e "${BLUE}[STEP ${1}]${NC} $(date '+%Y-%m-%d %H:%M:%S') ${2}"
}

check_tools() {
    log_info "检查必要工具..."
    local tools=("fastp" "bwa" "samtools" "gatk" "parallel" "tabix")
    local missing_tools=()
    
    for tool in "${tools[@]}"; do
        if ! command -v "$tool" &> /dev/null; then
            missing_tools+=("$tool")
        fi
    done
    
    if [[ ${#missing_tools[@]} -gt 0 ]]; then
        log_error "以下工具未找到: ${missing_tools[*]}"
        log_error "请确保已安装并添加到PATH环境变量"
        exit 1
    fi
    
    log_info "所有工具检查通过"
}

check_reference() {
    log_info "检查参考基因组及索引..."
    
    if [[ ! -f "$REF" ]]; then
        log_error "参考基因组文件不存在: $REF"
        exit 1
    fi
    
    # 生成 samtools faidx 索引
    if [[ ! -f "${REF}.fai" ]]; then
        log_info "生成 samtools faidx 索引..."
        samtools faidx "$REF"
    fi
    
    # 生成 bwa 索引
    if [[ ! -f "${REF}.bwt" ]]; then
        log_info "生成 bwa 索引..."
        bwa index "$REF"
    fi
    
    # 生成 GATK 序列字典
    local dict="${REF%.*}.dict"
    if [[ ! -f "$dict" ]]; then
        log_info "生成 GATK 序列字典..."
        gatk CreateSequenceDictionary -R "$REF"
    fi
    
    log_info "参考基因组索引检查完成"
}

create_directories() {
    log_info "创建输出目录..."
    mkdir -p "$FILTER_DIR"
    mkdir -p "$BAM_DIR"
    mkdir -p "$GVCF_DIR"
    mkdir -p "$DB_DIR"
    log_info "输出目录创建完成"
}

# ==================== Step 1: 质量过滤 ====================

run_fastp() {
    local R1_FILE=$1
    local BASENAME=$(basename "$R1_FILE" _R1.fastq.gz)
    local R2_FILE="${INPUT_DIR}/${BASENAME}_R2.fastq.gz"
    
    if [[ ! -f "$R2_FILE" ]]; then
        log_error "R2文件不存在: $R2_FILE，跳过样本 $BASENAME"
        return 1
    fi
    
    log_info "处理样本: $BASENAME"
    
    fastp \
        -i "$R1_FILE" \
        -I "$R2_FILE" \
        -o "${FILTER_DIR}/${BASENAME}_filter_R1.fastq.gz" \
        -O "${FILTER_DIR}/${BASENAME}_filter_R2.fastq.gz" \
        -l ${MIN_LENGTH} \
        --thread ${FASTP_THREADS} \
        --detect_adapter_for_pe \
        --cut_front \
        --cut_front_window_size 1 \
        --cut_front_mean_quality ${QUALITY_CUTOFF} \
        --cut_tail \
        --cut_tail_window_size 1 \
        --cut_tail_mean_quality ${QUALITY_CUTOFF} \
        --json "${FILTER_DIR}/${BASENAME}_fastp.json" \
        --html "${FILTER_DIR}/${BASENAME}_fastp.html"
    
    if [[ $? -eq 0 ]]; then
        log_info "样本 $BASENAME 质量过滤完成"
    else
        log_error "样本 $BASENAME 质量过滤失败"
        return 1
    fi
}

step1_filter() {
    log_step 1 "质量过滤 (fastp) 开始..."
    
    local input_files=("${INPUT_DIR}"/*_R1.fastq.gz)
    local total_samples=${#input_files[@]}
    
    if [[ $total_samples -eq 0 ]]; then
        log_error "未找到输入文件: ${INPUT_DIR}/*_R1.fastq.gz"
        exit 1
    fi
    
    log_info "共发现 $total_samples 个样本"
    
    export -f run_fastp
    export INPUT_DIR FILTER_DIR MIN_LENGTH FASTP_THREADS QUALITY_CUTOFF
    export log_info log_error
    
    parallel -j ${PARALLEL_JOBS} run_fastp {} ::: "${INPUT_DIR}"/*_R1.fastq.gz
    
    log_step 1 "质量过滤完成"
}

# ==================== Step 2: 比对与标记重复 ====================

run_mapping() {
    local R1_FILE=$1
    local BASENAME=$(basename "$R1_FILE" _filter_R1.fastq.gz)
    local R2_FILE="${FILTER_DIR}/${BASENAME}_filter_R2.fastq.gz"
    
    if [[ ! -f "$R2_FILE" ]]; then
        log_error "R2文件不存在: $R2_FILE，跳过样本 $BASENAME"
        return 1
    fi
    
    mkdir -p "${BAM_DIR}/${BASENAME}"
    log_info "处理样本: $BASENAME"
    
    # BWA 比对并排序
    log_info "[$BASENAME] 执行 BWA 比对..."
    bwa mem -t ${BWA_THREADS} \
        -R "@RG\tID:${BASENAME}\tPL:ILLUMINA\tSM:${BASENAME}" \
        "$REF" "$R1_FILE" "$R2_FILE" | \
        samtools sort --threads ${SAMTOOLS_THREADS} \
        -o "${BAM_DIR}/${BASENAME}/${BASENAME}.sorted.bam"
    
    # 标记重复序列
    log_info "[$BASENAME] 标记重复序列..."
    gatk MarkDuplicates \
        -I "${BAM_DIR}/${BASENAME}/${BASENAME}.sorted.bam" \
        -O "${BAM_DIR}/${BASENAME}/${BASENAME}.marked_duplicates.bam" \
        -M "${BAM_DIR}/${BASENAME}/${BASENAME}.marked_dup_metrics.txt" \
        --OPTICAL_DUPLICATE_PIXEL_DISTANCE 2500
    
    # 清理中间文件
    rm -f "${BAM_DIR}/${BASENAME}/${BASENAME}.sorted.bam"
    
    log_info "样本 $BASENAME 比对完成"
}

step2_mapping() {
    log_step 2 "比对与标记重复 (bwa + GATK MarkDuplicates) 开始..."
    
    check_reference
    
    local input_files=("${FILTER_DIR}"/*_filter_R1.fastq.gz)
    local total_samples=${#input_files[@]}
    
    if [[ $total_samples -eq 0 ]]; then
        log_error "未找到过滤后的文件: ${FILTER_DIR}/*_filter_R1.fastq.gz"
        exit 1
    fi
    
    log_info "共发现 $total_samples 个样本"
    
    export -f run_mapping
    export FILTER_DIR BAM_DIR REF BWA_THREADS SAMTOOLS_THREADS
    export log_info log_error
    
    parallel -j ${PARALLEL_JOBS} run_mapping {} ::: "${FILTER_DIR}"/*_filter_R1.fastq.gz
    
    log_step 2 "比对完成"
}

# ==================== Step 3: GVCF生成 ====================

run_gvcf() {
    local R1_FILE=$1
    local BASENAME=$(basename "$R1_FILE" _filter_R1.fastq.gz)
    local BAM_FILE="${BAM_DIR}/${BASENAME}/${BASENAME}.marked_duplicates.bam"
    
    if [[ ! -f "$BAM_FILE" ]]; then
        log_error "BAM文件不存在: $BAM_FILE，跳过样本 $BASENAME"
        return 1
    fi
    
    mkdir -p "${GVCF_DIR}/${BASENAME}"
    log_info "处理样本: $BASENAME"
    
    # 建立 BAM 索引
    log_info "[$BASENAME] 建立 BAM 索引..."
    samtools index "$BAM_FILE"
    
    # 质控统计
    samtools flagstat "$BAM_FILE" > "${BAM_DIR}/${BASENAME}/${BASENAME}.flagstat"
    
    # 生成 GVCF
    log_info "[$BASENAME] 生成 GVCF..."
    gatk HaplotypeCaller \
        -R "$REF" \
        --emit-ref-confidence GVCF \
        -I "$BAM_FILE" \
        -O "${GVCF_DIR}/${BASENAME}/${BASENAME}.g.vcf.gz" \
        --sample-name "$BASENAME"
    
    # 建立 GVCF 索引
    tabix -p vcf "${GVCF_DIR}/${BASENAME}/${BASENAME}.g.vcf.gz"
    
    log_info "样本 $BASENAME GVCF 生成完成"
}

step3_gvcf() {
    log_step 3 "GVCF生成 (GATK HaplotypeCaller) 开始..."
    
    local input_files=("${FILTER_DIR}"/*_filter_R1.fastq.gz)
    local total_samples=${#input_files[@]}
    
    if [[ $total_samples -eq 0 ]]; then
        log_error "未找到过滤后的文件: ${FILTER_DIR}/*_filter_R1.fastq.gz"
        exit 1
    fi
    
    log_info "共发现 $total_samples 个样本"
    
    export -f run_gvcf
    export BAM_DIR GVCF_DIR REF
    export log_info log_error
    
    parallel -j ${PARALLEL_JOBS} run_gvcf {} ::: "${FILTER_DIR}"/*_filter_R1.fastq.gz
    
    log_step 3 "GVCF生成完成"
}

# ==================== Step 4: GenomicsDB构建 ====================

run_db_import() {
    local chr=$1
    local tmp_dir="${DB_DIR}/tmp_${chr}"
    mkdir -p "$tmp_dir"
    
    log_info "处理染色体: $chr"
    
    gatk --java-options "-Xms${GATK_XMS} -Xmx${GATK_XMX} -XX:ParallelGCThreads=${GATK_THREADS}" \
        GenomicsDBImport \
        --batch-size 10 \
        --reader-threads ${GATK_THREADS} \
        -L "$chr" \
        --sample-name-map "${DB_DIR}/sample_map.txt" \
        --genomicsdb-workspace-path "${DB_DIR}/DB_${chr}" \
        --tmp-dir "$tmp_dir" \
        --genomicsdb-shared-posixfs-optimizations true \
        --reference "$REF"
    
    # 清理临时目录
    rm -rf "$tmp_dir"
    
    log_info "染色体 $chr GenomicsDB 构建完成"
}

step4_db() {
    log_step 4 "GenomicsDB构建 (GATK GenomicsDBImport) 开始..."
    
    # 生成样本映射文件
    local map_file="${DB_DIR}/sample_map.txt"
    log_info "生成样本映射文件..."
    
    find "$GVCF_DIR" -name "*.g.vcf.gz" | while read -r file; do
        local sample=$(basename "$file" .g.vcf.gz)
        echo -e "${sample}\t${file}"
    done > "$map_file"
    
    local sample_count=$(wc -l < "$map_file")
    log_info "共发现 $sample_count 个样本"
    
    if [[ $sample_count -eq 0 ]]; then
        log_error "未找到 GVCF 文件"
        exit 1
    fi
    
    export -f run_db_import
    export DB_DIR GVCF_DIR REF GATK_THREADS GATK_XMS GATK_XMX
    export log_info log_error
    
    parallel -j ${PARALLEL_JOBS} run_db_import {} ::: "${CHROMOSOMES[@]}"
    
    log_step 4 "GenomicsDB构建完成"
}

# ==================== Step 5: 联合基因分型与合并 ====================

run_genotype() {
    local chr=$1
    local out_dir="${DB_DIR}/${chr}"
    local tmp_dir="${DB_DIR}/tmp_genotype_${chr}"
    
    mkdir -p "$out_dir" "$tmp_dir"
    
    log_info "处理染色体: $chr"
    
    gatk --java-options "-Xms${GATK_XMS} -Xmx${GATK_XMX} -XX:ParallelGCThreads=${GATK_THREADS}" \
        GenotypeGVCFs \
        -R "$REF" \
        -V "gendb://${DB_DIR}/DB_${chr}" \
        -O "${out_dir}/${chr}.vcf.gz" \
        --allow-old-rms-mapping-quality-annotation-data \
        --max-alternate-alleles 4 \
        --genomicsdb-max-alternate-alleles 10 \
        --genomicsdb-shared-posixfs-optimizations \
        --tmp-dir "$tmp_dir"
    
    # 清理临时目录
    rm -rf "$tmp_dir"
    
    log_info "染色体 $chr 基因分型完成"
}

step5_genotype() {
    log_step 5 "联合基因分型 (GATK GenotypeGVCFs) 开始..."
    
    export -f run_genotype
    export DB_DIR REF GATK_THREADS GATK_XMS GATK_XMX
    export log_info log_error
    
    parallel -j ${PARALLEL_JOBS} run_genotype {} ::: "${CHROMOSOMES[@]}"
    
    # 合并所有染色体的 VCF 文件
    log_info "合并所有染色体 VCF 文件..."
    
    local merge_inputs=""
    for chr in "${CHROMOSOMES[@]}"; do
        merge_inputs+=" -I ${DB_DIR}/${chr}/${chr}.vcf.gz"
    done
    
    gatk --java-options "-Xms${GATK_XMS} -Xmx${GATK_XMX}" \
        MergeVcfs \
        $merge_inputs \
        -O "${OUTPUT_DIR}/final_variants.vcf.gz"
    
    # 建立最终 VCF 索引
    tabix -p vcf "${OUTPUT_DIR}/final_variants.vcf.gz"
    
    log_step 5 "联合基因分型完成"
    log_info "最终 VCF 文件: ${OUTPUT_DIR}/final_variants.vcf.gz"
}

# ==================== 主函数 ====================

main() {
    local start_step=${1:-1}
    local end_step=${2:-5}
    
    echo "=========================================="
    echo "  Variants Calling Pipeline"
    echo "=========================================="
    echo "输入目录: $INPUT_DIR"
    echo "输出目录: $OUTPUT_DIR"
    echo "参考基因组: $REF"
    echo "并行任务数: ${PARALLEL_JOBS}"
    echo "执行步骤: ${start_step} - ${end_step}"
    echo "=========================================="
    echo ""
    
    check_tools
    create_directories
    
    if [[ $start_step -le 1 ]] && [[ $end_step -ge 1 ]]; then step1_filter; fi
    if [[ $start_step -le 2 ]] && [[ $end_step -ge 2 ]]; then step2_mapping; fi
    if [[ $start_step -le 3 ]] && [[ $end_step -ge 3 ]]; then step3_gvcf; fi
    if [[ $start_step -le 4 ]] && [[ $end_step -ge 4 ]]; then step4_db; fi
    if [[ $start_step -le 5 ]] && [[ $end_step -ge 5 ]]; then step5_genotype; fi
    
    echo ""
    echo "=========================================="
    log_info "流程全部完成！"
    echo "=========================================="
}

# ==================== 帮助信息 ====================

usage() {
    cat << EOF
用法: $0 [选项] [起始步骤] [结束步骤]

Variants Calling Pipeline - 基于GATK最佳实践的变异检测流程

选项:
    -h, --help      显示此帮助信息
    -a, --all       执行全部步骤 (1-5)
    -s, --step N    仅执行单个步骤 N

步骤说明:
    1   质量过滤 (fastp)
        - 过滤低质量序列
        - 自动检测并切除接头
        - 前后端质量修剪

    2   比对与标记重复 (bwa + GATK MarkDuplicates)
        - BWA比对到参考基因组
        - 标记PCR重复序列

    3   GVCF生成 (GATK HaplotypeCaller)
        - 为每个样本生成GVCF文件
        - 包含所有位点信息

    4   GenomicsDB构建 (GATK GenomicsDBImport)
        - 按染色体导入GVCF到GenomicsDB
        - 便于大规模样本处理

    5   联合基因分型与合并 (GATK GenotypeGVCFs + MergeVcfs)
        - 跨样本联合基因分型
        - 合并所有染色体VCF

示例:
    $0                      # 执行全部步骤 (1-5)
    $0 2 4                  # 执行步骤2到4
    $0 -s 3                 # 仅执行步骤3
    $0 -a                   # 执行全部步骤

配置:
    请修改脚本顶部的"用户配置区"设置以下参数:
    - INPUT_DIR     : 原始fastq文件目录
    - OUTPUT_DIR    : 输出根目录
    - REF           : 参考基因组路径
    - CHROMOSOMES   : 染色体名称列表
    - PARALLEL_JOBS : 并行任务数
    - GATK_XMX      : GATK最大内存

输出文件:
    最终结果: \${OUTPUT_DIR}/final_variants.vcf.gz
    各步骤中间文件位于:
      - 01_Filter/  : 过滤后的fastq文件
      - 02_BAM/      : 比对结果BAM文件
      - 03_GVCF/     : GVCF文件
      - 04_DB/       : GenomicsDB数据库

EOF
}

# ==================== 参数解析与执行 ====================

case "${1:-}" in
    -h|--help)
        usage
        exit 0
        ;;
    -a|--all)
        main 1 5
        ;;
    -s|--step)
        if [[ -z "${2:-}" ]]; then
            log_error "请指定步骤编号"
            usage
            exit 1
        fi
        main "$2" "$2"
        ;;
    "")
        main 1 5
        ;;
    *)
        if [[ "$1" =~ ^[0-9]+$ ]]; then
            main "$1" "${2:-$1}"
        else
            log_error "无效参数: $1"
            usage
            exit 1
        fi
        ;;
esac