此处与云养江停的文章有略微差距

一.SNP calling

参照文章gvcf文件与vcf文件 - 简书 (jianshu.com)
可知GVCF文件与VCF文件的差异

这里可采用两种方法
方法1：
是直接执行一次HaplotypeCaller，适用于单样本或固定样本数量的情况。而如果此时增加一次样本，则需重新运行整个HaplotypeCaller，为了解决这个问题，样本数量不固定的采用方法2
方法2：
是每个样本先各自生成gVCF，然后再进行群体joint-genotype

（具体详见文章WGS（重测序）分析详解与脚本-麦子-的博客-CSDN博客_重测序分析）

因为我有多个样本，此处采用方法2

1.为每个样本生成GVCF文件

gatk HaplotypeCaller  -R reference.fasta --emit-ref-confidence GVCF -L chr19 -I 1.sorted.markdup.bam -O 1.g.vcf

然而，基因组上各个不同的染色体之间其实是可以理解为相互独立的（结构性变异除外），也就是说，为了提高效率我们可以按照染色体一条条来独立执行这个步骤，最后再把结果合并起来就好了，这样的话就能够节省很多的时间。增加了一个 -L 参数，通过这个参数我们可以指定特定的染色体（或者基因组区域)

例如L 1指定的是 1 号染色体，有些地方会写成chr1，具体看human.fasta中如何命名，与其保持一致即可。。

2.合并GVCF文件

首先，对GVCF文件进行压缩
压缩命令如下，需要进行批量压缩，下面提供了批量压缩的脚本

bgzip -c -f -@ 10 sample.vcf > sample.vcf.gz

-c, --stdout            write on standard output, keep original files unchanged
-f, --force             overwrite files without asking
-@, --threads INT       number of compression threads to use [1]

以下提供了批量压缩的脚本

madir zff_txt
cd zff_test
vi vcf_zip.sh
ls *.vcf >1.txt
cat 1.txt可查看该文件，此步非必须
cat 1.txt | while read line
do
        arr=$line
        bgzip -c -f -@10 ${arr[0]} > ${arr[0]}.gz
done

madir是新建一个文件夹
cd 打开这个文件
vi 对这个脚本命名，后面必须跟后缀.sh
*表示匹配任意字符,“ls *.vcf” 表示展示所有后缀带 ".vcf" 的文件， “>1.txt”表示将展示的文件名输出为txt文件
arr没有含义，相当于函数中X的作用
$表示引用
do和done之间的命令在开头要加tab键
0表示txt文件中的第一列，若是1，则表示第二列，2则表示第三列

其次，合并

ls*.g.vcf.gz > all_gvcf 
gatk CombineGVCFs -R ref.fa -V all_gvcf.list  -O merged.g.vcf.gz
(炫9明太小明郎君)

3.生成基因型文件，这一步生成VCF文件了

gatk  GenotypeGVCFs -R reference.fasta -V 1.g.vcf -O test.vcf

二、挑选SNP及过滤

在GATK HaplotypeCaller之后，⾸选的质控⽅案是GATK VQSR， 原理是利⽤⾃⾝数据和已知变异位点集的overlap进行质控。因为我研究的物质是鱼，没有已知变异集，因此只能进行硬过滤或自己构建变异集。这里采用硬过滤。

1.选择出SNP

gatk SelectVariants -select-type SNP -V sheep.raw.vcf.gz -O sheep.raw.snp.vcf.gz

2.过滤

SNP位点过滤前需要问自己一个问题，我的数据需要过滤吗？

一般要看后期是否做关联分析（GWAS）；如果只是单纯研究群体结构建议不过滤，因为过滤掉低频位点可能会改变某些样本之间的关系；如果需要和表型联系其来做关联分析，那么建议过滤，因为在后期分析中低频位点是不在考虑范围内的，需要保持前后一致

作者：EwanH
链接：https://www.jianshu.com/p/2677dc3b1383
来源：简书
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

先进行粗过滤

gatk VariantFiltration \
-R ref.fa \
-V sheep.raw.snp.vcf.gz \
--filter-expression "QUAL < 30.0 || QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" \
--filter-name LowQua \
-O sheep.snp.filter.vcf.gz

根据过滤结果选SNP

gatk SelectVariants 、
-R ref.fa \
-V sheep.snp.filter.vcf.gz \
--exclude-filtered --restrict-alleles-to BIALLELIC \
-O  sheep.snp.pass.vcf.gz

（炫9明太小明郎君9）

进一步过滤

这里我的过滤条件为：MAF>0.05,J基因缺失率<0.15，只保留bi-allelic SNP

bcftools view -i 'F_MISSING < 15' -q 0.05:minor -m2 -M2 test.vcf.gz -Oz -o test.flt.vcf.gz

使用 bcftools 进行基因型过滤（genotypes QC） - 橙子牛奶糖 - 博客园 (cnblogs.com)
及SNP过滤 - 简书 (jianshu.io)

参数

-i/e: --include/--exclude <expr> select/exclude sites for which the expression is true
-m/M: --min-alleles/--max-alleles <int>; -m2 -M2即保留二肽SNP
-O: --output-type <b|u|z|v> b: compressed BCF, u: uncompressed BCF, z: compressed VCF, v: uncompressed VCF
0.05:minor 指的是次等位基因频率大于0.05

最⼩等位基因频率（Minor Allele frequencies）
MAF是最小等位基因频率,通常是指在给定人群中的不常见的等位基因发生频率，例如TT，TC，CC三个基因型，在人群中C的频率=0.36，T的频率=0.64，则等位基因C就为最小等位基因频率，MAF=0.36。
Hapmap计划将MAF>0.05的SNPs作为首要研究目标，MAF广泛应用于复杂疾病的全基因组关联研究。在关联研究中，较小的MAF将会使统计效能降低，从而造成假阴性的结果。为了研究人群中罕见突变与疾病的关联，通常通过加大样本量的方法来弥补MAF降低所带来的统计效能的损失。在研究中，这一突变位点往往属于同一基因或同一通路上一组罕见突变中的一个。

此外，已知MAF还可以估算样本量或检验效能，并可以确定基因型的频率。

(百度百科MAF_百度百科 (baidu.com)）

次要等位基因的count数⼤于3: 相当于MAF=0.01
(作者：0玉萧元甲资料集)

缺失⽐例（missing rates）
假如缺失⽐例为0.05，总共100个个体的情况下，则该SNP在100*0.05=5个个体中丢失。
bi-allelic位点

为什么⼀般只保留bi-allelic SNP，要去除multi-allelic SNP？
bi-allelic位点是指基因组的某个位置上有两个allele，其中参考基因组在该位点上的碱基算作⼀个allele，样本在该位置上的变异算作⼀个allele。所以bi-allelic 位点即该位点只有⼀种变异。例如下图所⽰的位点7只有⼀种变异，样本1-3的deletion。

image.png

⽽下图所⽰的位点7则是⼀个multi-allelic位点，有参考基因组的G和样本2的C和样本3的T两种SNP。

image.png

三、计算Fst

fst的介绍见Fst的含义、计算与应用 - 简书 (jianshu.com)

vcf 转为 ped/map

plink格式的ped和map文件及转化为012的方法 (360doc.com)

vcftools --vcf snp.vcf --plink --out snp

计算fst

vcftools --vcf test.vcf --weir-fst-pop 1_population.txt --weir-fst-pop 2_population.txt --out p_1_2—single

此处为单点计算fst的方法
test.vcf是SNP calling 过滤后生成的vcf 文件；
p_1_2_3 生成结果的prefix
1_population.txt是一个文件包含同一个群体中所有个体，一般每行一个个体。个体名字要和vcf的名字对应。2_population.txt 包含了群体二中所有个体。

作者：资料降龙纪花大全
链接：https://wenku.baidu.com/view/ce3d5a23b868a98271fe910ef12d2af90242a8ba.html
来源：百度文库
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

过滤的更多命令详见下文章
群体变异数据vcf文件过滤概念及使用方法 - 简书 (jianshu.com)

作者：0玉萧元甲资料集0
链接：https://wenku.baidu.com/view/0277ac175427a5e9856a561252d380eb62942332.html

作者：炫9明太小明郎君9
链接：https://wenku.baidu.com/view/3d8ebd366f175f0e7cd184254b35eefdc8d3157b.html