3. SNP Calling

以GATK流程为例，简单来说，SNP Calling主要包括以下几步：

1. 给参考基因组建立索引：samtools faidx、bwa index, gatk CreateSequenceDictionary
2. 比对到参考基因组，生成SAM文件：bwa mem
3. 排序并生成BAM文件，并对BAM文件进行PCR重复标记（这一步生成MarkDup.bam）：samtools sort、gatk MarkDuplicates
4. 对MarkDup.bam建立索引：samtools index
5. SNP Calling（产生gVCF文件）：gatk HaplotypeCaller
6. 合并多个样本gVCFs，生成VCF文件：gatk GenomicsDBImport、gatk GenotypeGVCFs

7. （未完待续）

SNP Calling的基本步骤

1. 比对到参考基因组；
2. 找到差异。

可能遇到的问题：
a.比对位置是否唯一？
*Reads很短
*基因组很复杂
b.测序错误？
*Reads数量大

3. 使用软件：GATK (Genome Analysis Toolkit)

   
   GATK

• GATK主要包括以下几种流程：

  
  GATK主要流程

• 当进行BSA-seq分析时，主要使用Germline SNPs & Indel流程：

  
  GATK-Germline SNPs & Indel pipeline
  

  3.1 GATK语法结构：

gatk ToolName --argument-name value
#argument names follow a "kebab" convention，即开头两个破折号，中间用单个破折号隔开

#添加Java参数，需要将参数写在双引号内，"-Xmx"指定了内存分配的数量：
gatk --java-options "-Xmx4G" [program arguments]

#一个栗子：
gatk HaplotypeCaller -R reference.fasta -I sample1.bam -O variants.vcf -ERC GVCF

#若使其更可读，可以用反斜杠隔开：
gatk HaplotypeCaller \
    -R reference.fasta \
    -I sample1.bam \
    -O variants.g.vcf \
    -ERC GVCF

一个实操：

拿SARS-CoV作为参考序列，新冠病毒测序数据作为练手：
1.创建环境

conda create -n SNPCalling gatk4 bwa samtools vcftools freebayes bcftools snpeff
conda activate SNPCalling

2.序列下载与预处理：

mkdir 2019-CoV
cd 2019-CoV/
mkdir ref
#参考基因组选用SARS序列，命名为sars_seq.fasta
mkdir raw_data
cd raw_data/
prefetch SRR11085733 -O ./
prefetch SRR11092056 -O ./

fastq-dump --split-files ./SRR11085733.sra
fastq-dump --split-files ./SRR11092056.sra

3.SNP calling
首先用bwa进行序列比对，产生bam文件，再进行SNP calling。目前，能够进行SNP calling的主流软件包括：GATK、bcftools，和freebayes，分别用这三款软件进行SNP calling，流程上参考了：

• https://www.jianshu.com/p/a84ff2b2ab59
• https://cloud.tencent.com/developer/article/1445600
  在此感谢！

#!/usr/bin/bash
##定义变量
COV=/home/Leon_Yang/data/2019-CoV
REF=$COV/ref/sars_seq.fasta
SRR1=SRR11085733
SRR2=SRR11092056
BAM1=$COV/align/$SRR1.bam
BAM2=$COV/align/$SRR1.bam
R1_1=$COV/raw_data/${SRR1}_1.fastq
R1_2=$COV/raw_data/${SRR1}_2.fastq
R2_1=$COV/raw_data/${SRR2}_1.fastq
R2_2=$COV/raw_data/${SRR2}_2.fastq
TAG1="@RG\tID:$SRR1\tSM:$SRR1\tLB:$SRR1"
TAG2="@RG\tID:$SRR2\tSM:$SRR2\tLB:$SRR2"
VARI=$COV/variant

##bwa比对，samtools排序并构建索引，为了之后更快调用比对文件
mkdir -p $COV/align && cd $COV/align
bwa index $REF
samtools faidx $REF
bwa mem -R $TAG1 $REF $R1_1 $R1_2 | samtools sort > $BAM1
bwa mem -R $TAG2 $REF $R2_1 $R2_2 | samtools sort > $BAM2
samtools index $BAM1
samtools index $BAM2

mkdir -p $VARI

##SNP calling
##第一种方法：GATK4（版本4.1.4.1）
#注意：使用GATK之前，需要先建立参考基因组索引文件.dict和.fai
#.dict中包含了基因组中contigs的名字，也就是一个字典；
#.fai也就是fasta index file，索引文件，可以快速找出参考基因组的碱基，由samtools faidx构建
#构建.dict文件（原来要使用picard的CreateSequenceDictionary模块，但是现在gatk整合了此模块，可以直接使用）
gatk --java-options "-Xmx8G -Djava.io.tmpdir=./tmp" CreateSequenceDictionary \
    -R $REF \
    -O $COV/ref/sars_seq.dict

#标记PCR重复reads
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./tmp" MarkDuplicates \
    -I $BAM1 -O ${SRR1}_MarkDup.bam \
    -M ${SRR1}.metrics
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./tmp" MarkDuplicates \
    -I $BAM2 -O ${SRR2}_MarkDup.bam \
    -M ${SRR2}.metrics

#用samtools对标记后的文件创建索引
samtools index ${SRR1}_MarkDup.bam
samtools index ${SRR2}_MarkDup.bam

#gatk开始：必选 -I -O -R，代表输入、输出、参考
#接下来可以按照字母顺序依次写出来，这样比较清晰
#-bamout：将一整套经过gatk程序重新组装的单倍体基因型（haplotypes）输出到文件
#-stand-call-conf :低于这个数字的变异位点被忽略，可以设成标准30（默认是10）
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./tmp" HaplotypeCaller \
    -R $REF \
    -I ${SRR1}_MarkDup.bam -O ${VARI}/${SRR1}_gatk.gvcf \
    --emit-ref-confidence GVCF \
    -stand-call-conf 30
gatk --java-options "-Xmx20G -Djava.io.tmpdir=./tmp" HaplotypeCaller \
    -R $REF \
    -I ${SRR2}_MarkDup.bam -O ${VARI}/${SRR2}_gatk.gvcf \
    --emit-ref-confidence GVCF \
    -stand-call-conf 30

#通常，GATK的HaplotypeCaller流程跑完后，得到的是单个样本的gVCF文件，即全基因组层面的变异信息，根据需要进行样本间gVCF文件的合并，并转化为VCF文件：
#合并多样本gVCF文件
gatk GenomicsDBImport \
    -L 1 \
    -V ${VARI}/${SRR1}_gatk.gvcf \
    -V ${VARI}/${SRR2}_gatk.gvcf \
    --genomicsdb-workspace-path $VARI/vcfdb
#gVCF转为VCF
gatk GenotypeGVCFs \
    -R $REF \
    -V gendb://$VARI/vcfdb \
    -O ${VARI}/HaplotypeCaller_gatk.vcf

##第二种方法：bcftools
#首先用samtools mpileup分析参考序列上每个碱基位点的比对结果，并生成VCF文件；
#再使用bcftools对VCF文件进行SNP calling 
samtools mpileup -uvf $REF $BAM1 $BAM2 | bcftools call -vm -Oz > ${VARI}/bcftools.vcf.gz
#参数说明：samtools: -v, 进行基因分型，结果输出到VCF格式，默认输出bgzip压缩格式文件，加-u后，不进行bgzip压缩，用于管道输入下一个命令
#参数说明：samtools: -f, 输入参考基因组序列fasta文件，fasta文件必须经过faidx产生fai索引文件
#参数说明：bcftools: -v, 仅输出变异位点信息，-m, 适合多种allele的算法，适合多样本，与-c只能二选一；
#参数说明：bcftools: -Oz, 设置输出文件格式，z为压缩VCF格式

##第三种方法：freebayes
freebayes -f $REF $BAM1 $BAM2 > ${VARI}/freebayes.vcf