一、胚系突变 SNP/INDEL（Germline SNPs + Indels）的介绍
胚系突变又叫生殖细胞突变，是来源于精子或卵子这些生殖细胞的突变，因此通常身上所有细胞都带有突变。胚系突变可以遗传，一般用于分析遗传病。当胚系突变发生在BRCA1，BRCA2等基因上的时候，就可能导致肿瘤的风险大大升高。这一步介绍的是GATK call 胚系突变的SNP和INDEL。
单个基因组是孤立的，通常并不能得出什么变异是致病性变异。一般要求家系或者群体数据，可以发现稀有变异，de novo变异，以及种族背景信息。这一步测试用到的是N/T单个样品分别进行分析。
GATK官网介绍：Germline short variant discovery (SNPs + Indels)
官方提供了六个流程：
Prod* germline short variant per-sample calling
Prod* germline short variant joint genotyping
$5 Genome Analysis Pipeline
Generic germline short variant per-sample calling
Generic germline short variant joint genotyping
Intel germline short variant per-sample calling
选取第四个流程（Generic germline short variant per-sample calling）进行学习
git地址：https://github.com/gatk-workflows/gatk4-germline-snps-indels
git clone https://github.com/gatk-workflows/gatk4-germline-snps-indels
下载后如下：
README.md
LICENSE
generic.google-papi.options.json
haplotypecaller-gvcf-gatk4.wdl（下面用到的）
haplotypecaller-gvcf-gatk4-nio.wdl
haplotypecaller-gvcf-gatk4.hg38.wgs.inputs.json
joint-discovery-gatk4-local.hg38.wgs.inputs.json
joint-discovery-gatk4-fc.wdl
joint-discovery-gatk4-local.wdl
joint-discovery-gatk4.hg38.wgs.inputs.json
joint-discovery-gatk4.wdl（下面用到的）
二、WDL中各个task的介绍
学习以上两个wdl，把这两个文件中的命令拆分出来。
主要是看task中的command。里面共包含了16个task（haplotypecaller-gvcf-gatk4.wdl包含了3个task/joint-discovery-gatk4.wdl包含了13个task），一一来看一下功能及命令：
haplotypecaller-gvcf-gatk4.wdl（3个）:
1. CramToBamTask
这一步可以省略，输入的bam直接采用数据预处理后的bam即可。如果输入是CRAM格式，才需要用到这一步。
2. HaplotypeCaller
gatk --java-options "-Xmx6G -XX:GCTimeLimit=50 -XX:GCHeapFreeLimit=10" \
HaplotypeCaller \
-R ref_fasta \
-I input_bam \
-L bed \
-O output_filename.vcf \
-contamination 0 \
--dbsnp dbsnp_138.hg19.vcf (-ERC GVCF)
--dbsnp 这个参数加上，结果中第三列ID会显示对应的RS号 GATK3如何选择选择known sites数据库。如果没有这个参数，第三列ID为.
-ERC GVCF是输出GVCF格式，默认是输出VCF格式。GVCF格式对于群体call突变比较适用，这里只是一个样品，所以用不到-ERC GVCF这个参数。
GVCF输出格式的一些介绍见官网：
What is a GVCF and how is it different from a 'regular' VCF?
HaplotypeCaller Reference Confidence Model (GVCF mode)
3. MergeGVCFs
功能：Merge GVCFs generated per-interval for the same sample
joint-discovery-gatk4.wdl（13个）：
1. GetNumberOfSamples
单个样品 用不到 忽略
2. ImportGVCFs
单个样品 用不到 忽略
3. GenotypeGVCFs
功能：合并每个样本的突变信息到单一vcf文件来方便进行下一步的过滤分析。
参考：http://www.php-master.com/post/326294.html
单个样品 用不到 忽略
4. HardFilterAndMakeSitesOnlyVcf
这一步包含了两步，其中一步是：
gatk --java-options "-Xmx3g -Xms3g" VariantFiltration \
--filter-expression "ExcessHet >54.69" --filter-name ExcessHet \
-O samplename_N.HaplotypeCaller.Filter.vcf \
-V samplename_T.HaplotypeCaller.vcf
这一步的作用是给VCF第七列FILTER增加PASS或者Filter，如果没有这一步，FILTER列会是.
运行完这一步，不满足条件的记录也会出现在最终生成的vcf文件中， 只不过对应的Filter字段的信息不是PASS。
--filter-expression是筛选条件 过滤掉ExcessHet > ${excess_het_threshold}的行。
但是经过测试，发现用了这一步以后，结果反而不准确，所以忽略这一步。（IGV核对）
参考：
https://www.jianshu.com/p/6f3198b7a070
https://www.plob.org/article/7023.html
5. IndelsVariantRecalibrator
建模，用于之后的VQSR
6. SNPsVariantRecalibratorCreateModel
建模，用于之后的VQSR
7. SNPsVariantRecalibrator
建模，用于之后的VQSR
8. GatherTranches
9. ApplyRecalibration
功能：Train on high-confidence known sites to determine the probability that other sites are true or false
样本数量较少的WES/非全基因组测序无法使用VQSR进行质控（https://blog.csdn.net/hanli1992/article/details/83652352），所以这一步在我的测试中省略。
官网对于VQSR介绍：Variant Quality Score Recalibration (VQSR)
10. GatherVcfs
11. CollectVariantCallingMetrics
12. GatherMetrics
13. DynamicallyCombineIntervals
三、执行命令
这一步对于单样品来讲比较简单，只需要运行HaplotypeCaller一步，再用SelectVariants筛选出SNP和INDEL即可。
实际运行命令（以normal为例）：
time gatk \
--java-options "-Xmx6G -XX:GCTimeLimit=50 -XX:GCHeapFreeLimit=10" \
HaplotypeCaller \
-R ref_hg19/ucsc.hg19.fasta -I samplename_N.recal.bam \
-L Covered.bed -O samplename_N.HaplotypeCaller.vcf \
-contamination 0 --dbsnp dbsnp_138.hg19.vcf
time gatk SelectVariants -V samplename_N.HaplotypeCaller.vcf -O samplename_N.Snp.Filter.vcf -select-type SNP
time gatk SelectVariants -V samplename_N.HaplotypeCaller.vcf -O samplename_N.Indel.Filter.vcf -select-type INDEL
normal所用时间：
13m0.890s
1m8.137s
1m8.058s
得到：
samplename_N.HaplotypeCaller.vcf
samplename_N.HaplotypeCaller.vcf.idx
samplename_N.Snp.Filter.vcf
samplename_N.Snp.Filter.vcf.idx
samplename_N.Indel.Filter.vcf
samplename_N.Indel.Filter.vcf.idx
对应的tumour所用时间：
7m53.815s
1m8.085s
1m10.476s
这一步的学习基本完成。
下一个step是对体细胞突变 SNV/INDEL（Somatic SNVs + Indels）的学习
网友评论