SNPs marker是全基因組范圍應用廣泛的分子標記，本文介紹生態(tài)基因組學中利用GATK4軟件進行SNPs calling的流程（人的研究中可能略有不同）。以下所有分析過程以GX_01這個樣本為例子衰絮。如果有多個樣本，使用for循環(huán)時需要注意引號內變量的使用扼雏。導出多個腳本也可南捂。
1攒砖、質量控制
數(shù)據(jù)質量檢測：使用fastqc軟件
fastqc -o output_dir GX_01.1.fq.gz GX_01.2.fq.gz ... seqfileN.fq.gz
# 后面可以接多個fq.gz文件
數(shù)據(jù)過濾：使用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為線程數(shù)日裙，-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(參考http://www.reibang.com/p/f190f9dfac03)
就是利用機器學習的方式調整原始堿基的質量分數(shù)（需要參考基因組有已知SNP的vcf文件财著，對大多數(shù)動物基因組而言不存在联四，人的研究中用的到，基本可以略過）撑教。
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
# 根據(jù)重校準表格導出校準后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收苏，需要進一步過濾）也就是說現(xiàn)在每個樣本一個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愤兵，需要進一步過濾）也就是說現(xiàn)在每個樣本一個gvcf文件
————
2021.4.3更新：有朋友私信我說在HaplotypeCaller這一步要加一個“--max-mnp-distance 0”選項鹿霸，否則后面CombineGVCF會報錯。我好像沒有遇到這種情況秆乳，遇到的朋友可以試一下這個解決方法懦鼠。感謝“大海龜啦啦啦”的建議。
————
8屹堰、將多個樣本的gvcf文件整合成一個文件（假設除了GX_01肛冶，我們還跑了GB_01，GK_01這兩個樣本双藕，現(xiàn)在一共三個樣本進行整合）
如果進行了第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 \ # 最大丟失個體數(shù)為零
--recode --recode-INFO-all \
--hwe 0.001 # 去除不滿足哈溫平衡 (p<0.001)的位點
--out final.vcf.gz
注意叶堆，過濾后的位點只是打著過濾標簽，還沒有被移除斥杜，可以用gvcftools處理虱颗。
