Skip to content

2020.11.17
SNPs calling流程（GATK4）

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处理。