SNPs marker是全基因组范围应用广泛的分子标记，本文介绍生态基因组学中利用GATK4软件进行SNPs calling的流程（人的研究中可能略有不同）。以下所有分析过程以GX_01这个样本为例子。如果有多个样本，使用for循环时需要注意引号内变量的使用。导出多个脚本也可。
1、质量控制
数据质量检测：使用fastqc软件
fastqc -o output_dir GX_01.1.fq.gz GX_01.2.fq.gz ... seqfileN.fq.gz
# 后面可以接多个fq.gz文件
数据过滤：使用fastp软件
fastp -I GX_01.1.fq.gz -o GX_01.out.1.fq.gz \
-I GX_01.2.fq.gz -O GX_01.out.2.fq.gz
# .1.fq.gz和.2.fq.gz分别为双端测序的两个文件
2、使用bwa软件进行reads mapping
# bwa建立参考基因组的引索
bwa index -a bwtsw reference.fa
#bwa mem比对算法
bwa mem -t 30 -M -P -R '@RG\tID:GX_01\tSM:GX_01\tLB:GX_01\tPL:Illumina' reference.fa GX_01.out.1.fq.gz GX_01.out.2.fq.gz >GX_01.0.sam
# -t为线程数，-R指定引号内的flag信息。ID一定要每个样本对号入座。
3、将sam转换为bam，并进行整理排序
samtools view -bS GX_01.0.sam -o GX_01.1.bam
#将sam转换为bam
samtools sort GX_01.1.bam -o GX_01.sort.2.bam
#对bam进行整理排序
4、标记Duplicates(PCA重复)
java -jar picard.jar MarkDuplicates REMOVE_DUPLICATES= false \
MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000 \
INPUT=GX_01.sort.2.bam OUTPUT=GX_01.marked.3.bam \
METRICS_FILE=GX_01.marked.3.bam.metrics
5、对新的到的bam文件建立引索
samtools index GX_01.marked.3.bam
6、Base (Quality Score) Recalibration(参考https://www.jianshu.com/p/f190f9dfac03)
就是利用机器学习的方式调整原始碱基的质量分数（需要参考基因组有已知SNP的vcf文件，对大多数动物基因组而言不存在，人的研究中用的到，基本可以略过）。
java -Xmx128g -jar gatk-package-4.0.10.1-local.jar BaseRecalibrator \
-R reference.fa -I GX_01.marked.3.bam --known-sites reference.known.vcf \
-O GX_01_recal_data.table
# 输出重校准表格
java -Xmx128g -jar gatk-package-4.0.10.1-local.jar ApplyBQSR \
-R reference.fa -I GX_01.marked.3.bam -bqsr GX_01_recal_data.table \
-O GX_01_recal_reads.bam
# 根据重校准表格导出校准后bam文件
#检测上述生成的bam文件是否可用。
java -Xmx128g -jar gatk-package-4.0.10.1-local.jar ValidateSamFile \
-I GX_01_recal_reads.bam
如果显示no errors found,则可以用HaplotypeCaller call SNP/Indel.
7、GATK4 call variants
如果进行了第6步的话（也就是说如果有参考SNPs库的话）
java -Xmx128G -jar gatk-package-4.0.10.1-local.jar HaplotypeCaller \
-R reference.fa -I GX_01_recal_reads.bam \
--dbsnp reference.known.vcf --native-pair-hmm-threads 96 \
-stand-call-conf 30 -O GX_01.g.vcf.gz -ERC GVCF
# 生成gvcf文件（gvcf格式包括所有的变异类型，包括snps和indel，需要进一步过滤）也就是说现在每个样本一个gvcf文件
如果没有进行第6步（一般动植物不会有SNPs库）
java -Xmx128G -jar gatk-package-4.0.10.1-local.jar HaplotypeCaller \
-R reference.fa -I GX_01.marked.3.bam -O GX_01.g.vcf.gz \
-ERC GVCF -ploidy 2 -stand_call_conf 30.0
# 生成gvcf文件（gvcf格式包括所有的变异类型，包括snps和indel，需要进一步过滤）也就是说现在每个样本一个gvcf文件
————
2021.4.3更新：有朋友私信我说在HaplotypeCaller这一步要加一个“--max-mnp-distance 0”选项，否则后面CombineGVCF会报错。我好像没有遇到这种情况，遇到的朋友可以试一下这个解决方法。感谢“大海龟啦啦啦”的建议。
————
8、将多个样本的gvcf文件整合成一个文件（假设除了GX_01，我们还跑了GB_01，GK_01这两个样本，现在一共三个样本进行整合）
如果进行了第6步的话（也就是说如果有参考SNPs库的话）
java -Xmx128G -jar gatk-package-4.0.10.1-local.jar CombineGVCFs \
-R reference.fa --dbsnp reference.known.vcf --variant GX_01.g.vcf.gz \
--variant GB_01.g.vcf.gz --variant GK_01.g.vcf.gz -O cohort.g.vcf.gz
如果没有进行第6步
java -Xmx128G -jar gatk-package-4.0.10.1-local.jar CombineGVCFs \
-R reference.fa --variant GX_01.g.vcf.gz --variant GB_01.g.vcf.gz \
--variant GK_01.g.vcf.gz -O cohort.g.vcf.gz
# 或者也可以不用一个一个指定gvcf文件，可以用-V gvcf.list
9、Genotyping(提取基因型)
java -Xmx128G -jar gatk-package-4.0.10.1-local.jar GenotypeGVCFs \
-R reference.fa -V cohort.g.vcf.gz -O genotype.vcf.gz
10、提取SNPs（因为上一个文件包含SNPs和indels）
java -Xmx96g -jar gatk-package-4.0.10.1-local.jar SelectVariants \
-R reference.fa -O SNPs.vcf --variant genotype.vcf.gz \
--select-type-to-include SNP
# 同理，也可以提取indels
11、对SNPs进行机械过滤
按照官方建议的来就行
gatk VariantFiltration \
-V SNPs.vcf \
-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 SNPs_filtered.vcf
12、还可以用python脚本、VCFtools进一步筛选
vcftools \
--vcf SNPs_filtered.vcf \
--maf 0.05 \ # 最小等位基因频率为 0.05，去除稀有等位基因 --max-missing-count 0 \ # 最大丢失个体数为零
--recode --recode-INFO-all \
--hwe 0.001 # 去除不满足哈温平衡 (p<0.001)的位点
--out final.vcf.gz
注意，过滤后的位点只是打着过滤标签，还没有被移除，可以用gvcftools处理。
