SNPs marker是全基因組范圍應(yīng)用廣泛的分子標(biāo)記，本文介紹生態(tài)基因組學(xué)中利用GATK4軟件進(jìn)行SNPs calling的流程（人的研究中可能略有不同）。以下所有分析過(guò)程以GX_01這個(gè)樣本為例子。如果有多個(gè)樣本，使用for循環(huán)時(shí)需要注意引號(hào)內(nèi)變量的使用。導(dǎo)出多個(gè)腳本也可。

1、質(zhì)量控制

數(shù)據(jù)質(zhì)量檢測(cè)：使用fastqc軟件

fastqc -o output_dir GX_01.1.fq.gz GX_01.2.fq.gz ... seqfileN.fq.gz
# 后面可以接多個(gè)fq.gz文件

數(shù)據(jù)過(guò)濾：使用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分別為雙端測(cè)序的兩個(gè)文件

2、使用bwa軟件進(jìn)行reads mapping

# bwa建立參考基因組的引索
bwa index -a bwtsw reference.fa

#bwa mem比對(duì)算法
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為線程數(shù)，-R指定引號(hào)內(nèi)的flag信息。ID一定要每個(gè)樣本對(duì)號(hào)入座。

3、將sam轉(zhuǎn)換為bam，并進(jìn)行整理排序

samtools view -bS GX_01.0.sam -o GX_01.1.bam
#將sam轉(zhuǎn)換為bam

samtools sort GX_01.1.bam -o GX_01.sort.2.bam
#對(duì)bam進(jìn)行整理排序

4、標(biāo)記Duplicates(PCA重復(fù))

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、對(duì)新的到的bam文件建立引索

samtools index GX_01.marked.3.bam

6、Base (Quality Score) Recalibration(參考http://www.itdecent.cn/p/f190f9dfac03)

就是利用機(jī)器學(xué)習(xí)的方式調(diào)整原始?jí)A基的質(zhì)量分?jǐn)?shù)（需要參考基因組有已知SNP的vcf文件，對(duì)大多數(shù)動(dòng)物基因組而言不存在，人的研究中用的到，基本可以略過(guò)）。

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
# 輸出重校準(zhǔn)表格

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
# 根據(jù)重校準(zhǔn)表格導(dǎo)出校準(zhǔn)后bam文件

#檢測(cè)上述生成的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

如果進(jìn)行了第6步的話（也就是說(shuō)如果有參考SNPs庫(kù)的話）

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格式包括所有的變異類(lèi)型，包括snps和indel，需要進(jìn)一步過(guò)濾）也就是說(shuō)現(xiàn)在每個(gè)樣本一個(gè)gvcf文件

如果沒(méi)有進(jìn)行第6步（一般動(dòng)植物不會(huì)有SNPs庫(kù)）

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格式包括所有的變異類(lèi)型，包括snps和indel，需要進(jìn)一步過(guò)濾）也就是說(shuō)現(xiàn)在每個(gè)樣本一個(gè)gvcf文件

————
2021.4.3更新：有朋友私信我說(shuō)在HaplotypeCaller這一步要加一個(gè)“--max-mnp-distance 0”選項(xiàng)，否則后面CombineGVCF會(huì)報(bào)錯(cuò)。我好像沒(méi)有遇到這種情況，遇到的朋友可以試一下這個(gè)解決方法。感謝“大海龜啦啦啦”的建議。
————

8、將多個(gè)樣本的gvcf文件整合成一個(gè)文件（假設(shè)除了GX_01，我們還跑了GB_01，GK_01這兩個(gè)樣本，現(xiàn)在一共三個(gè)樣本進(jìn)行整合）

如果進(jìn)行了第6步的話（也就是說(shuō)如果有參考SNPs庫(kù)的話）

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 

如果沒(méi)有進(jìn)行第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 
# 或者也可以不用一個(gè)一個(gè)指定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（因?yàn)樯弦粋€(gè)文件包含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、對(duì)SNPs進(jìn)行機(jī)械過(guò)濾

按照官方建議的來(lái)就行

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進(jìn)一步篩選

vcftools \
--vcf SNPs_filtered.vcf \
--maf 0.05 \ # 最小等位基因頻率為 0.05，去除稀有等位基因 --max-missing-count 0 \ # 最大丟失個(gè)體數(shù)為零
--recode --recode-INFO-all \
--hwe 0.001 # 去除不滿足哈溫平衡 (p<0.001)的位點(diǎn)
--out final.vcf.gz

注意，過(guò)濾后的位點(diǎn)只是打著過(guò)濾標(biāo)簽，還沒(méi)有被移除，可以用gvcftools處理。