对全基因重测序数据进行变异检测获取 SNP 、indel、SV 等变异信息，是对群体进行遗传学分析前的必须步骤。如果检测个体较多或测序深度（测序量）较大，那么该过程还需要计算机具有较大的 RAM 容量以及高性能 CPU。

进行变异检测需要用到不同的软件，并且实际执行的步骤在不同的操作者手中也各有差异。下面介绍一种我自己执行后可用的SNP calling 步骤教程，仅供参考。

软件和环境准备：

conda create -n snp_calling python=3.9 
conda activate snp_calling
conda install -c bioconda bwa hisat2 samtools picard gatk4

这里我使用 conda 管理软件，同时提前创建好分析环境。

1.比对

首先对参考基因组进行索引，这里使用 bwa 进行。

#使用BWA对参考基因组进行索引，这里使用100个线程进行索引
bwa index -t 100 ref_genome.fna

我使用的计算服务器具有 384 个 CPU 线程，这里使用 100 个线程计算。实际数字需要根据你使用的计算机或者服务器的实际 CPU 核心和线程数决定，不要超过线程数。

索引完成后，我们使用下面的命令比对

bwa mem -t 10 ref_genome.fna R1.clean.fq.gz R2.clean.fq.gz 01.sam

该命令表示使用 10 个线程进行比对，R1.clean.fq.gz、R2.clean.fq.gz分别代表已质控的不同方向的双端测序下机数据，01.sam 是结果输出文件。相关路径可以自行再进行设置。

如果个体数量比较多，可以使用 wait 命令多线程同时对所有个体进行比对。

2.转换和排序SAM/BAM文件

使用 samtools 处理时，标准的命令如下：

samtools view -@ 4 -bS 01.sam > 01.bam
samtools sort -@ 4 01.bam -o 01_sorted.bam
samtools index 01_sorted.bam

同样的，面对个体数量较多的情况，我们可以将整个过程多线程且自动化，可以使用下面的脚本：

#Step 1: Convert .sam to .bam
#这里识别输入文件夹d的所有的 sam 文件，所以请提前把上一步比对sh生成的 sam 文件放在同一个目录下
#注意这一步输出的目录也要设置号
for file in "$input_dir"/*.sam; do
  samtools view -@ 12 -bS "$file" > "$output_dir_que/$(basename ${file%.sam}.bam)" &
done
wait
#Step 2: Sort .bam files
for file in "$output_dir_que"/*.bam; do
  samtools sort -@ 12 "$file" -o "$output_dir_modify/$(basename ${file%.bam}_sorted.bam)" &
done
wait
#Step 3: Index sorted .bam files
for file in "$output_dir_modify"/*_sorted.bam; do
  samtools index "$file" &
done
wait
echo "运行结束"

上面的脚本，三个步骤，不同步骤的输出目录可以在脚本里修改好后运行，也可以通过bash 命令将变量送入。

另外线程数一样需要保证不要大于 CPU 最大线程数。如果你使用脚本，还需要注意，单个命令的线程数乘以你的个体数量得到的总线程需要数，不能大于你的 CPU 线程数。

3.标记 PCR 重复

这里我们是用到的是Picard，需要注意的是，Java openjdk的版本需要大于17。如果出现了Java版本问题，可以使用conda构建合适的环境：

conda install -c conda-forge openjdk=17

标准的标记重复和索引的命令：

picard MarkDuplicates I=01_sorted.bam O=01_marked.bam M=01_metrics.txt
samtools index 01_marked.bam

使用脚本自动化执行的代码如下：

#!/bin/bash
#设置picard工具的路径，这里根据实际情况修改
PICARD_JAR="/你的picard.jar的路径/picard.jar"
#创建输出目录
mkdir -p modify/marked_bams
#定义一个函数进行标记重复
mark_duplicates() {
    local bam_file=$1
    #这里教程中命名是 01.bam或者01_sorted.bam，这样的格式。意思是如果你有很多个体，后面的就是 02、03 这样子。如果你的文件名不是这样的格式，需要修改下面的代码
    local base_name=(basename bam_file "_sorted.bam")
    #这里自动识别目录中以_sorted.bam结尾的文件，然后去掉_sorted.bam，得到文件名
    local output_bam="modify/marked_bams/${base_name}_marked.bam"
    #这里定义输出文件名
    local metrics_file="modify/marked_bams/${base_name}_metrics.txt"
    #这里定义输出metrics文件名
    java -Xmx32g -jar $PICARD_JAR MarkDuplicates         I=$bam_file         O=$output_bam         M=$metrics_file         CREATE_INDEX=true 2>&1 | tee modify/marked_bams/${base_name}_markdup.log
    #这里调用picard工具进行标记重复
    if [ ! -f $output_bam ]; then
        echo "Failed to generate output_bam, check modify/marked_bams/{base_name}_markdup.log for details." >> error_log.txt
    fi
    #如果输出文件不存在，则输出错误信息
}
#并行处理5个任务，也可以根据需要进行修改。可以通过调整并行数量，然后尝试运行，如果2-3分钟未报错中止，则可正常进行
max_jobs=5
current_jobs=0
for bam_file in modify/*_sorted.bam; do
    mark_duplicates $bam_file &
    current_jobs=$((current_jobs + 1))
    if [ current_jobs -ge max_jobs ]; then
        wait
        current_jobs=0
    fi
done
#等待所有后台任务完成
wait
echo "所有BAM文件的重复标记已完成"

4.添加样本信息到 BAM 文件

在实际处理中，对标记后的bam文件进行处理时，GATK 4会报错，未发现样本信息。所以需要将样本信息标记到bam文件中。使用Pichard进行标记

picard AddOrReplaceReadGroups  I=modify/marked_bams/01_marked.bam  O=modify/marked_bams/01_marked_with_RG.bam RGID=1 RGLB=lib1 RGPL=illumina RGPU=unit1 RGSM=sample_name 
//这里的sample_name等信息需要根据你的实际情况针对性修改

下面是一个进行处理的脚本，如果你需要添加的信息不一样，做针对性修改即可。

#!/bin/bash

input_dir="modify/marked_bams"
output_dir="modify/marked_bams_with_RG"

mkdir -p $output_dir

for bam_file in $input_dir/*_marked.bam; do
    sample_name=$(basename $bam_file _marked.bam)
    output_file="$output_dir/${sample_name}_marked_with_RG.bam"
  
    picard AddOrReplaceReadGroups \
        I=$bam_file \
        O=$output_file \
        RGID=1 \
        RGLB=lib1 \
        RGPL=illumina \
        RGPU=unit1 \
        RGSM=$sample_name
  
    samtools index $output_file
done

5.GATK变异检测

使用GATK进行SNP和INDEL检测。下面是其中一个样本的实际分析过程中的命令，使用15个CPU线程。如果你需要获取 GVCF 文件，请跳过本部分代码，朝后翻看。

gatk HaplotypeCaller 
   -R ../../genome_data/ref_genome.fna 
   -I modify/marked_bams_with_RG/01_marked_with_RG.bam 
   -O output_vcf/01.vcf 
   --native-pair-hmm-threads 15

下面是多线程操作脚本

#!/bin/bash
input_dir="modify/marked_bams_with_RG"
output_dir="output_vcf"
reference_genome="../../genome_data/ref_genome.fna"
mkdir -p $output_dir
for bam_file in $input_dir/*_marked_with_RG.bam; do

    sample_name=$(basename $bam_file _marked_with_RG.bam)

    output_vcf="$output_dir/${sample_name}.vcf"
gatk HaplotypeCaller \

       -R $reference_genome \

       -I $bam_file \

       -O $output_vcf \

       --native-pair-hmm-threads 15 &
done
wait

耗时参考信息：我自己的本次变异检测工作，每个个体的SNP检测使用15个CPU线程，最终完成时间为12000个CPU分钟，共200个CPU小时，即13.3小时。2.2GB基因组大小，测序深度35X测序双端序列文件大小各为7GB。

上面的命令没有获取 GVCF 文件，下面是获取 GVCF 的版本。

gatk HaplotypeCaller 
   -R ../../genome_data/ref_genome.fna 
   -I modify/marked_bams_with_RG/01_marked_with_RG.bam 
   -O output_vcf/01.vcf 
   --native-pair-hmm-threads 15 \ -ERC GVCF -ploidy 2 -stand-call-conf 30.0

批处理脚本：

#!/bin/bash
input_dir="modify/marked_bams_with_RG"
output_dir="output_vcf"
reference_genome="../../genome_data/ref_genome.fna"
mkdir -p $output_dir
for bam_file in $input_dir/*_marked_with_RG.bam; do

    sample_name=$(basename $bam_file _marked_with_RG.bam)
output_vcf="$output_dir/${sample_name}.g.vcf.gz"
gatk HaplotypeCaller \

       -R $reference_genome \

       -I $bam_file \

       -O $output_vcf \

       -ERC GVCF -ploidy 2 -stand-call-conf 30.0 \

       --native-pair-hmm-threads 15 &
done
wait

6.个体变异信息VCF文件合并

gatk CombineGVCFs -R ../../genome_data/ref_genome.fna 
--variant output_vcf/01.g.vcf.gz 
--variant output_vcf/02.g.vcf.gz
--variant ...
-O ./final_combine/combined.g.vcf.gz

如果你的文件数量很多，你可以分批次合并加快速度。例如如果你有 100 个个体，你可以一次合并 10 个，最后一起合并 10 个。

7.Genotyping & 提取 SNP

使用下面命令对变异位点进行基因分型。

gatk GenotypeGVCFs -R ../../genome_data/ref_genome.vcf -V ./final_combine/combined.g.vcf.gz -O ./final_genotype/genotype.vcf.gz

使用下面的命令提取 SNP

gatk SelectVariants -R ../../genome_data/ref_genome.vcf -O SNP.vcf.gz --variant ./final_genotype/genotype.vcf.gz --select-type-to-include SNP

8.SNP 机械过滤

机械过滤根据自己的需求设置相关参数和过滤标准，下面是一个参考

gatk VariantFiltration \     
-V SNP.vcf.gz \
-filter "QD < 2.0" --filter-name "QD2" \
-filter "QUAL < 30.0" --filter-name "QUAL30" \    
-filter "SOR > 3.0" --filter-name "SOR3" \    
-filter "FS > 60.0" --filter-name "FS60" \   
-filter "MQ < 40.0" --filter-name "MQ40" \    
-filter "MQRankSum < -12.5" --filter-name "MQRankSum-12.5" \   
-filter "ReadPosRankSum < -8.0" --filter-name "ReadPosRankSum-8" \
-O SNP_filtered.vcf.gz

9.后续操作

获取 SNP.vcf文件后，可以在使用 vcftools 进行分析前的过滤等操作。

菜单

重测序数据SNP Calling 教程

分享

重测序数据SNP Calling 教程

1.比对

2.转换和排序SAM/BAM文件

3.标记 PCR 重复

4.添加样本信息到 BAM 文件

5.GATK变异检测

6.个体变异信息VCF文件合并

7.Genotyping & 提取 SNP

8.SNP 机械过滤

9.后续操作

庆祝新中国成立 75 周年！

Excalidraw: 轻量强大的自由画布 (路线图、演示图绘制推荐)

重测序数据SNP Calling 教程

基于腾讯云EO实现对被封443端口IP的 443 端口访问

国产良心！GLM-4-Flash 免费模型接入 Zotero

测序技术发展对动物遗传育种的一些影响

Article: Multiple populations' genetic structure research(WGRS base)

重测序群体绘制 NJ Tree 进化树

Zotero 翻译插件的安装并基于 AI 模型调用翻译

若干国内访问 ChatGPT 的方法