以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
- （未完待续）
SNP Calling的基本步骤
- 比对到参考基因组；
- 找到差异。
可能遇到的问题：
a.比对位置是否唯一？
*Reads很短
*基因组很复杂
b.测序错误？
*Reads数量大
-
使用软件：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，流程上参考了：
#!/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
