进行群体遗传学分析的前提是我们得得到变异信息，之前已经介绍了如何进行
CNV的识别，那么本文来介绍下如何识别基因组中分布更多、更常见的变异SNP（INDEL顺带）
二代测序 reads 的比对
二代测序 reads 的比对分两步：
用 bwa 将 reads 比对到参考基因组
用 picard 去除 PCR 造成的重复 reads
最终得到比对后的 CRAM 文件，这是一种BAM的压缩格式，在 samtools 给的基准测试中，CRAM 大小约为 BAM 的一半。
bwa 软件实现了三种比对算法 BWA-backtrack, BWA-SW 和 BWA-MEM。
- 第一种算法适用于长度在 100bp 以下的 reads，后两种算法适用于70bp至数M的长 reads。
- BWA-MEM 是最新的算法，70bp以上的 reads 用
BWA-MEM就好，70b p以下的用BWA-backtrack。
01. read alignment, sort, and remove PCR duplication
这一步也就是我们常用的数据处理了，包括比对、排序、标记PCR重复
bwa mem -t 4 -M -R '@RG\tID:$sample\tLB:$group\tPL:ILLUMINA\tSM:$sample' $reference.fa $fastq1 $fastq2 | samtools view -b -S -o $sample.bam
java -jar picard.jar SortSam I=$sample.bam O=$sample.sort.bam SORT_ORDER=coordinate VALIDATION_STRINGENCY=LENIENT
java -jar picard.jar MarkDuplicates I=$sample.sort.bam O=$sample.dup.bam ASSUME_SORT_ORDER=coordinate METRICS_FILE=$sample.dup.txt VALIDATION_STRINGENCY=LENIENT REMOVE_DUPLICATES=true
java -jar picard.jar SetNmMdAndUqTags I=$sample.dup.bam O=$sample.dup.sort.fix.bam CREATE_INDEX=true R=$reference.fa VALIDATION_STRINGENCY=LENIENT
# REMOVE_DUPLICATES=true 去除PCR重复
# SetNmMdAndUqTags 设置SAM/BAM文件中的 NM(比对质量)、MD(碱基差异)和 UQ(唯一比对性) 标签,后续过滤可能用得上
# NM(比对质量) > 20 的用于后续 SNP calling
# bwa 中 -R 参数添加头文件信息，
若比对时未加这一参数则可用：
java -jar picard.jar AddOrReplaceReadGroups I=input.bam O=output.bam RGID=$sample RGLB=$sample RGPL=illumina RGPU=$samplePU RGSM=$sample
GATK进行SNP calling
image.png
02.variants calling using GATK4+
GATK4版本之后就不用对INDEL区域进行重比对了，方便
~/gatk-4.1.4.0/gatk HaplotypeCaller -R $reference.fa -L $chr -ERC GVCF -I $bam -o $chr.gvcf.gz
~/gatk-4.1.4.0/gatk CombineGVCFs -R $reference.fa --variant $chr.list -o $chr.gvcf.gz
~/gatk-4.1.4.0/gatk GenotypeGVCFs -R $reference.fa --variant $chr.gvcf.gz --includeNonVariantSites -O $chr.vcf.gz
# -L 分染色体call，占用空间少, CombineGVCFs时指定$chr.list即可
# --includeNonVariantSite 保留非变异位点
-
GATK的HaplotypeCaller模块识别每个个体中基因组变异SNP和INDEL -
GATK的CombineGVCFs模块对分染色体文件及个体文件进行合并，最后得到群体层面的大gvcf -
GATK的GenotypeGVCFs模块对群体层面的大gvcf进行基因分型，得到基因型
02.variants calling using GATK4+（hard filter）
- 对SNP和INDEL分别进行硬过滤，因为他们过滤标准不一样
-
GATK的SelectVariants模块分别提取出SNP和INDEL - 过滤完成后，用
MergeVcfs模块对过滤后vcf文件进行合并。
SNP
~/gatk-4.1.4.0/gatk SelectVariants -R $reference.fa -V $Pop.vcf.gz --select-type-to-include SNP -O Pop.SNP.vcf.gz
~/gatk-4.1.4.0/gatk VariantFiltration -R $reference.fa -V Pop.SNP.vcf.gz --filter-expression "QD < 2.0 || FS > 60.0 || MQ < 40.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "my_snp_filter" -O Pop.HDflt.SNP.vcf.gz
INDEL
~/gatk-4.1.4.0/gatk SelectVariants -R $reference.fa -V $Pop.vcf.gz --select-type-to-include INDEL -O Pop.INDEL.vcf.gz
~/gatk-4.1.4.0/gatk VariantFiltration -R $reference.fa -V Pop.INDEL.vcf.gz --filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -20" --filter-name "my_indel_filter" -O Pop.HDflt.INDEL.vcf.gz
merge and get PASS site
~/gatk-4.1.4.0/gatk MergeVcfs -I Pop.HDflt.SNP.vcf.gz -I Pop.HDflt.INDEL.vcf.gz -O Pop.HDflt.SNP.INDEL.genotype.vcf
~/gatk-4.1.4.0/gatk SelectVariants -R $reference.fa -V Pop.HDflt.SNP.INDEL.genotype.vcf -O Pop.HDflt.SNP.INDEL.genotype.pass.vcf -select "vc.isNotFiltered()"
call完之后可根据质量值QUAL >= 30去除低质量位点，去除多等位基因等等过滤
vcftools --gzvcf $vcf --minQ 30 --min-alleles 2 --max-alleles 2 --max-missing 0.9 --non-ref-ac 2 --remove-indels --recode --recode-INFO-all --out $filter
--minQ 30 保留QUAL >= 30
--min-alleles 2 保留双等位基因
--max-alleles 2 保留双等位基因
--non-ref-ac 2 保留非参考等位基因多余2个个体的位点
--remove-indels 去除INDEL
--remove-snps 去除SNP
--max-missing 0.9 缺失率 < 10% 的位点留下
整个硬过滤及质量控制流程HardFilterSnpIndel.sh
#! /bin/bash
## 精准而优雅
#################### Hard filter for SNP ######################
# 如果想要位点多点，则最后--minQ可不设置
# vcftools --non-ref-ac 2 保留非参考等位基因多余2个个体的位点
# 看自己情况设还是不设
###############################################################
# Set up the file name, software
GATK="/home/software/gatk-4.1.4.0/gatk" # change as you want
vcftools="/home/sll/miniconda3/bin/bcftools" # change as you want
reference="/home/sll/genome" # reference genome fasta file
if [[ $# -ne 2 ]]; then
echo "Usage: bash $0 <vcf.gz file> <outprefix>"
echo "输出文件前缀推荐和输入文件一样吧"
echo "对SNP和INDEL进行硬过滤及SNP去除多等位基因，输出过滤后的SNP及INDEL文件"
exit 1
fi
vcf=$1
out=$2
## get SNP and INDEL
$GATK SelectVariants -R $reference --select-type SNP -V $vcf -O ${out}.snps.vcf.gz
$GATK SelectVariants -R $reference --select-type INDEL -V $vcf -O ${out}.indel.vcf.gz
## filter SNP and INDEL
$GATK VariantFiltration -R $reference -V ${out}.snps.vcf.gz --filter-expression "QD < 2.0 || FS > 60.0 || MQ < 40.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "SNP_FILTER" -O ${out}.snps.HDflt.vcf.gz
$GATK VariantFiltration -R $reference -V ${out}.snps.vcf.gz --filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -20" --filter-name "INDEL_FILTER" -O ${out}.indel.HDflt.vcf.gz
## get pass for SNP and INDEL
$GATK MergeVcfs -I ${out}.snps.HDflt.vcf.gz -I ${out}.indel.HDflt.vcf.gz -O ${out}.SNP.INDEL.HDflt.vcf.gz
$GATK SelectVariants -R $reference -V ${out}.SNP.INDEL.HDflt.vcf.gz -O ${out}.SNP.INDEL.HDflt.pass.vcf.gz -select "vc.isNotFiltered()"
## get SNP and indel
$vcftools --gzvcf ${out}.SNP.INDEL.HDflt.pass.vcf.gz --minQ 30 --min-alleles 2 --max-alleles 2 --remove-indels --recode --recode-INFO-all --out ${out}.SNP.HDflt.pass
$vcftools --gzvcf ${out}.SNP.INDEL.HDflt.pass.vcf.gz --minQ 30 --min-alleles 2 --max-alleles 2 --remove-snps --recode --recode-INFO-all --out ${out}.INDEL.HDflt.pass
GATK最常用的是硬过滤，因为简单直接，但还有机器学习的方法VQSR那个更加适用于自己的数据
ANGSD进行SNP calling
这个就用来 call 出 SNP 之后与GATK结果取共同鉴定出的SNP，让SNP calling结果更准确,然后他输出的是基因型的似然值(GL)，而不是确定的基因型，但下游的大多分析都是基于genotype的，因此，还是将交集部分利用GATK进行genotype calling
angsd -bam $bam.list -only_proper_pairs 1 -uniqueOnly 1 -remove_bads 1 \
-minQ 20 -minMapQ 30 -C 50 -ref $reference.fa -r $chr \
-out $chr -doQsDist 1 -doDepth 1 -doCounts 1 -maxDepth $max
angsd -bam $bam.list -only_proper_pairs 1 -uniqueOnly 1 -skipTriallelic 1 -remove_bads 1 \
-minQ 20 -minMapQ 30 -C 50 -ref $reference.fa -r $chr -out $chr \
-doMaf 1 -doMajorMinor 1 -GL 1 -setMinDepth $mindepth -setMaxDepth $maxdepth \
-doCounts 1 -dosnpstat 1 -SNP_pval 1
# 第一步是call SNP，及计算多个参数
# 第二步是多个过滤参数及估计genotype likelihood（GL），保证SNP的准确性
# -SNP_pval 1 SNP的p_value，可用于过滤
# -doMajorMinor 1 最终文件无 ref alt 而是 major minor 形式
# -doMaf 1 1: Frequency (fixed major and minor)
# -GL 1 1: SAMtools
2: GATK
3: SOAPsnp
4: SYK
# -doGeno 进行genotype calling，
1：print out major minor
2: print the called genotype as -1,0,1,2 (count of minor)
4: print the called genotype as AA, AC, AG, ...
8: print all 3 posts (major,major),(major,minor),(minor,minor)
