本文中vcf文件的获取主要是使用gatk的HaplotypeCaller模块，该模块的原理大概如下：
- 定义active regions 沿着参考基因组以一定的窗口滑动，统计比对的 mismatches, indels 和 softclips等信息计算基因组每个位置的活跃得分，使用平滑算法进行，此处相当于测定该区域熵值。当熵值达到某个设定的阈值是即确定该区域为active region，用于后续组装。
- 通过active regions的组装确定单倍型 对于每个active regions，程序构建一个类似 De Bruijn 的图来重新组装active regions并识别 数据中可能存在哪些单倍型。然后，程序将每个单倍型与参考重新对齐 使用 Smith-Waterman 算法进行单倍型分析，以识别潜在的变异位点。
- 确定每个read的单倍型的似然值 对于每个active regions，该程序使用PairHMM算法对每个单倍型的每个读数进行成对比对。这就产生了给定读取数据的单倍型可能性矩阵。然后将这些可能性边缘化，以获得给定读取数据的每个潜在变异位点的等位基因的可能性。
- 分配样本基因型 对于每个潜在的变异位点，该程序应用贝叶斯规则，使用给定读取数据的等位基因的可能性来计算给定该样本的读取数据的每个样本的每个基因型的可能性。然后将最有可能的基因型分配给样本。
GATK
下载安装
wget https://github.com/broadinstitute/gatk/releases/download/4.2.1.0/gatk-4.4.0.0.zip #下载安装包
unzip gatk-4.2.1.0 #解压缩安装包
cd gatk-4.2.1.0/
#导入环境变量
echo "PATH=$PATH:~/softwarepath/gatk-4.2.1.0" >>~/.bashrc
source ~/.bashrc
#同样有懒人安装方法
conda install -c bioconda gatk4
#测试是否安装成功
gatk --help #检查是否可正常运行
Usage template for all tools (uses --spark-runner LOCAL when used with a Spark tool)
gatk AnyTool toolArgs
Usage template for Spark tools (will NOT work on non-Spark tools)
gatk SparkTool toolArgs [ -- --spark-runner <LOCAL | SPARK | GCS> sparkArgs ]
Getting help
gatk --list Print the list of available tools
gatk Tool --help Print help on a particular tool
Configuration File Specification
--gatk-config-file PATH/TO/GATK/PROPERTIES/FILE
gatk forwards commands to GATK and adds some sugar for submitting spark jobs
--spark-runner <target> controls how spark tools are run
valid targets are:
LOCAL: run using the in-memory spark runner
SPARK: run using spark-submit on an existing cluster
--spark-master must be specified
--spark-submit-command may be specified to control the Spark submit command
arguments to spark-submit may optionally be specified after --
GCS: run using Google cloud dataproc
commands after the -- will be passed to dataproc
--cluster <your-cluster> must be specified after the --
spark properties and some common spark-submit parameters will be translated
to dataproc equivalents
--dry-run may be specified to output the generated command line without running it
--java-options 'OPTION1[ OPTION2=Y ... ]' optional - pass the given string of options to the
java JVM at runtime.
Java options MUST be passed inside a single string with space-separated values.
使用
因为这边只是介绍该软件在BSA分析中如何使用，所以介绍相对单一，也是为了初学者不会被那么多冗余的信息干扰。普通变异检测同样适用，当样本太多的时候可以试试gatk4的GenomicsDB模块，代替这里的CombineGVCFs + GenotypeGVCFs策略。
BSA分析需要通过该软件获得最终用于分析的vcf文件，所以这一步可以说是入门的钥匙。
首先准备输入文件：混池数据比对结果bam文件和参考基因组，比对过程参考前面的推文
这里文件命名分别为：
！！！一定要区分清楚样本之间的对应关系，不能混了
P1.sort.rmdup.bam 亲本1混池比对结果文件
P2.sort.rmdup.bam 亲本2混池比对结果文件
S1.sort.rmdup.bam 子代1混池比对结果文件（与亲本1表型一致）
S2.sort.rmdup.bam 子代2混池比对结果文件（与亲本2表型一致）
ref.fa 参考基因组
# 建立索引
gatk CreateSequenceDictionary -R ref.fa \
-O ref.dict
# GVCF输出
gatk --java-options "-Xmx4g" HaplotypeCaller \
-R ref.fa \
-I P1.sort.rmdup.bam \
-O P1.g.vcf.gz \
-ERC GVCF
##参数解读
"-Xmx4g" 这个是指定给予的运行内存，可以比4g更大，但是要结合服务器本身的资源决定。
-R 参考基因组
-I bam文件
-O 输出文件
-ERC 输出GVCF的格式
##其他样本同样如此操作，分别得到对应的GVCF
P1.g.vcf.gz P2.g.vcf.gz S1.g.vcf.gz S2.g.vcf.gz
# 合并GVCF
gatk CombineGVCFs \
-R ref.fa \
-V P1.g.vcf.gz \
-V P2.g.vcf.gz \
-V S1.g.vcf.gz \
-V S2.g.vcf.gz \
-O merge.g.vcf.gz
##参数解读
-V 不同样本的gvcf文件
-O 合并后的gvcf文件
# GVCF转VCF
gatk --java-options "-Xmx4g" GenotypeGVCFs \
-R ref.fa \
-V merge.g.vcf.gz \
-O merge.vcf.gz
##参数解读
-V GVCF文件
-O VCF文件
结果文件格式解读
GVCF
##fileformat=VCFv4.2
##ALT=<ID=NON_REF,Description="Represents any possible alternative allele at this location">
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read depth">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=MIN_DP,Number=1,Type=Integer,Description="Minimum DP observed within the GVCF block">
##FORMAT=<ID=PGT,Number=1,Type=String,Description="Physical phasing haplotype information, describing how the alternate alleles are phased in relation to one another">
##FORMAT=<ID=PID,Number=1,Type=String,Description="Physical phasing ID information, where each unique ID within a given sample (but not across samples) connects records within a phasing group">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="The phred-scaled genotype likelihoods rounded to the closest integer">
##FORMAT=<ID=SB,Number=4,Type=Integer,Description="Per-sample component statistics which comprise the Fisher's Exact Test to detect strand bias">
##INFO=<ID=BaseQRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt Vs. Ref base qualities">
##INFO=<ID=ClippingRankSum,Number=1,Type=Float,Description="Z-score From Wilcoxon rank sum test of Alt vs. Ref number of hard clipped bases">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Combined depth across samples">
##INFO=<ID=END,Number=1,Type=Integer,Description="Stop position of the interval">
##INFO=<ID=ExcessHet,Number=1,Type=Float,Description="Phred-scaled p-value for exact test of excess heterozygosity">
##INFO=<ID=InbreedingCoeff,Number=1,Type=Float,Description="Inbreeding coefficient as estimated from the genotype likelihoods per-sample when compared against the Hardy-Weinberg expectation">
##INFO=<ID=MLEAC,Number=A,Type=Integer,Description="Maximum likelihood expectation (MLE) for the allele counts, for each ALT allele, in the same order as listed">
##INFO=<ID=MLEAF,Number=A,Type=Float,Description="Maximum likelihood expectation (MLE) for the allele frequency, for each ALT allele, in the same order as listed">
##INFO=<ID=MQRankSum,Number=1,Type=Float,Description="Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities">
##INFO=<ID=RAW_MQ,Number=1,Type=Float,Description="Raw data for RMS mapping quality">
##INFO=<ID=ReadPosRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias">
##CommandLine.Haplotyper=<ID=Haplotyper,Version="gatk.4.0",Date="2023-10-25T18:22:05Z",CommandLine="gatk --java-options "-Xmx4g" HaplotypeCaller -R ref.fa -I P1.sort.rmdup.bam -O P1.g.vcf.gz -ERC GVCFgatk --java-options "-Xmx4g" HaplotypeCaller -R ref.fa -I P1.sort.rmdup.bam -O P1.g.vcf.gz -ERC GVCF "
##contig=<ID=Chr01,length=91897250,assembly=unknown>
##contig=<ID=Chr02,length=34677250,assembly=unknown>
##contig=<ID=Chr03,length=58976660,assembly=unknown>
##contig=<ID=Chr04,length=23665986,assembly=unknown>
##contig=<ID=Chr05,length=32975341,assembly=unknown>
##contig=<ID=Chr06,length=89148132,assembly=unknown>
##reference=file:///path/to/ref.fa
##GVCFBlock0-1=minGQ=0(inclusive),maxGQ=1(exclusive)
##GVCFBlock1-2=minGQ=1(inclusive),maxGQ=2(exclusive)
##GVCFBlock2-3=minGQ=2(inclusive),maxGQ=3(exclusive)
##GVCFBlock3-4=minGQ=3(inclusive),maxGQ=4(exclusive)
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT CY180
Chr01 1 . A <NON_REF> . . END=5538 GT:DP:GQ:MIN_DP:PL 0/0:0:0:0:0,0,0
Chr01 5539 . A <NON_REF> . . END=5539 GT:DP:GQ:MIN_DP:PL 0/0:1:3:1:0,3,15
Chr01 5540 . C <NON_REF> . . END=5540 GT:DP:GQ:MIN_DP:PL 0/0:1:3:1:0,3,38
Chr01 5541 . C <NON_REF> . . END=5541 GT:DP:GQ:MIN_DP:PL 0/0:1:3:1:0,3,24
Chr01 5542 . C <NON_REF> . . END=5548 GT:DP:GQ:MIN_DP:PL 0/0:1:3:1:0,3,28
......
Chr01 10959 . A G,<NON_REF> 15.43 . DP=2;ExcessHet=3.0103;MLEAC=1,0;MLEAF=0.5,0;RAW_MQ=5200.00 GT:AD:DP:GQ:PL:SB 0/1:0,2,0:2:0:40,0,0,43,6,49:0,0,1,1ChrA01_S1 10959 . A G,<NON_REF> 15.43 . DP=2;ExcessHet=3.0103;MLEAC=1,0;MLEAF=0.5,0;RAW_MQ=5200.00 GT:AD:DP:GQ:PL:SB 0/1:0,2,0:2:0:40,0,0,43,6,49:0,0,1,1
"#"开头的都是注释说明行
##FORMAT=... --> FORMAT那一列对应的参数信息及解释，比如GT是字符串，表示基因型
##INFO=... --> INFO那一列对应的参数信息及解释，比如DP是整数，表示样本组合深度
##CommandLine... --> 执行命令行
##contig=... --> 染色体编号，长度信息
##reference=file:///path/to/ref.fa --> 参考基因组
##GVCFBlock0-1=minGQ=0(inclusive),maxGQ=1(exclusive) --> GVCF才有，VCF没有
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT P1 --> 标题行
以序列编号开头的即为文件主体内容
每列含义与标题行一致
第一列：序列编号
第二列:物理起始位置
第三列：变异位点的rsID号，如果没有的话用"."表示。
第四列：参考基因组在该位点的碱基，REF
第五列：样本在该位点的碱基，ALT
第六列：call出这个位点的质量。这个值等于-10log10(p)，p值是call错alt allele错误的概率。也就是QUAL越大出错概率越小。 第七列：对变异位点进行过滤，如果通过则为PASS，如果没有进行过滤就是"."。第八列 INFO：这一列是额外信息。可能是平台的信息，也可以是DP等的信息：
第八列：INFO列，里面具体参数解释在表头的注释行里。这里使用第一行数据作为例子着重说明一下，"END=5538"表示这一行表示的是从1到5538的变异信息。所以GVCF是包含所有位点的信息文件，不管位点是否存在变异。
第九列：FORMAT列，里面具体参数解释在表头的注释行里
第十列：样本位点信息列，数据以：为分隔符，数据含义与第九列保持一致。比如第九列是GT:DP:GQ:MIN_DP:PL，第十列是0/0:1:3:1:0,3,15，那第十列就表示该样本的GT为0/0，DP为1，GQ为3，MIN_DP为1，PL为0,3,15。
VCF
##fileformat=VCFv4.2
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read depth">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=PGT,Number=1,Type=String,Description="Physical phasing haplotype information, describing how the alternate alleles are phased in relation to one another">
##FORMAT=<ID=PID,Number=1,Type=String,Description="Physical phasing ID information, where each unique ID within a given sample (but not across samples) connects records within a phasing group">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="The phred-scaled genotype likelihoods rounded to the closest integer">
##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele frequency, for each ALT allele, in the same order as listed">
##fileformat=VCFv4.2
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read depth">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=PGT,Number=1,Type=String,Description="Physical phasing haplotype information, describing how the alternate alleles are phased in relation to one another">
##FORMAT=<ID=PID,Number=1,Type=String,Description="Physical phasing ID information, where each unique ID within a given sample (but not across samples) connects records within a phasing group">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="The phred-scaled genotype likelihoods rounded to the closest integer">
##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele frequency, for each ALT allele, in the same order as listed">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
##INFO=<ID=BaseQRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt Vs. Ref base qualities">
##INFO=<ID=ClippingRankSum,Number=1,Type=Float,Description="Z-score From Wilcoxon rank sum test of Alt vs. Ref number of hard clipped bases">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Combined depth across samples">
##INFO=<ID=ExcessHet,Number=1,Type=Float,Description="Phred-scaled p-value for exact test of excess heterozygosity">
##INFO=<ID=FS,Number=1,Type=Float,Description="Phred-scaled p-value using Fisher's exact test to detect strand bias">
##INFO=<ID=InbreedingCoeff,Number=1,Type=Float,Description="Inbreeding coefficient as estimated from the genotype likelihoods per-sample when compared against the Hardy-Weinberg expectation">
##INFO=<ID=MLEAC,Number=A,Type=Integer,Description="Maximum likelihood expectation (MLE) for the allele counts, for each ALT allele, in the same order as listed">
##INFO=<ID=MLEAF,Number=A,Type=Float,Description="Maximum likelihood expectation (MLE) for the allele frequency, for each ALT allele, in the same order as listed">
##INFO=<ID=MQ,Number=1,Type=Float,Description="RMS mapping quality">
##INFO=<ID=MQRankSum,Number=1,Type=Float,Description="Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities">
##INFO=<ID=QD,Number=1,Type=Float,Description="Variant Confidence/Quality by Depth">
##INFO=<ID=ReadPosRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias">
##INFO=<ID=SOR,Number=1,Type=Float,Description="Symmetric Odds Ratio of 2x2 contingency table to detect strand bias">
##CommandLine.GVCFtyper=<ID=GVCFtyper,Version="gatk.4.0",Date="2023-06-10T14:05:36Z",CommandLine="gatk --java-options "-Xmx4g" GenotypeGVCFs -R ref.fa -V merge.g.vcf.gz -O merge.vcf.gz
##contig=<ID=Chr01,length=112420854,assembly=unknown>
##contig=<ID=Chr02,length=103302290,assembly=unknown>
##contig=<ID=Chr03,length=143109472,assembly=unknown>
##contig=<ID=Chr04,length=128801742,assembly=unknown>
##contig=<ID=Chr05,length=116542366,assembly=unknown>
##contig=<ID=Chr06,length=118975115,assembly=unknown>
##reference=file:///path/to/ref.fa
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT P1 P2 S1
Chr01 6484 . TAAA T 41.66 . AC=1;AF=0.167;AN=6;BaseQRankSum=-1.383;ClippingRankSum=-0.000;DP=31;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.167;MQ=60.00;MQRankSum=-0.000;QD=10.42;ReadPosRankSum=-0.000;SOR=0.105 GT:AD:DP:GQ:PL ./.:0,0:0:.:. 0/0:17,0:17:9:0,9,531 0/0:10,0:10:0:0,0,147 0/1:2,2:4:46:78,0,46
Chr01 15370 . A T 47.64 . AC=1;AF=0.167;AN=6;BaseQRankSum=-1.383;ClippingRankSum=-0.000;DP=25;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.167;MQ=60.00;MQRankSum=-0.000;QD=11.91;ReadPosRankSum=-0.000;SOR=0.693 GT:AD:DP:GQ:PGT:PID:PL ./.:0,0:0:.:.:.:. 0/0:14,0:14:39:.:.:0,39,585 0/1:2,2:4:78:0|1:15370_A_T:78,0,143 0/0:7,0:7:18:.:.:0,18,270
Chr01 15377 . T C 47.60 . AC=1;AF=0.167;AN=6;BaseQRankSum=-1.383;ClippingRankSum=-0.000;DP=25;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.167;MQ=60.00;MQRankSum=-0.000;QD=11.90;ReadPosRankSum=-0.000;SOR=0.693 GT:AD:DP:GQ:PGT:PID:PL ./.:0,0:0:.:.:.:. 0/0:12,0:12:33:.:.:0,33,495 0/1:2,2:4:78:0|1:15370_A_T:78,0,143 0/0:9,0:9:24:.:.:0,24,360
Chr01 28237 . A AT 30.55 . AC=1;AF=0.167;AN=6;BaseQRankSum=-0.000;ClippingRankSum=-0.000;DP=35;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.167;MQ=60.00;MQRankSum=-0.000;QD=10.18;ReadPosRankSum=0.967;SOR=0.223 GT:AD:DP:GQ:PL ./.:0,0:0:.:. 0/0:20,0:20:60:0,60,768 0/0:12,0:12:33:0,33,495 0/1:1,2:3:30:70,0,30
Chr01 28861 . A AT 34.02 . AC=1;AF=0.167;AN=6;BaseQRankSum=-1.408;ClippingRankSum=-0.000;DP=37;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.167;MQ=59.30;MQRankSum=-0.812;QD=3.09;ReadPosRankSum=0.393;SOR=0.883 GT:AD:DP:GQ:PL ./.:0,0:0:.:. 0/0:19,0:19:14:0,14,605 0/0:6,0:6:18:0,18,205 0/1:7,4:12:72:72,0,169
Chr01 30412 . G A 34.35 . AC=1;AF=0.25;AN=4;BaseQRankSum=-0.812;ClippingRankSum=-0.000;DP=22;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.25;MQ=60.00;MQRankSum=-0.000;QD=3.82;ReadPosRankSum=-1.480;SOR=0.223 GT:AD:DP:GQ:PGT:PID:PL ./.:0,0:0:.:.:.:. 0/0:13,0:13:36:.:.:0,36,540 0/1:7,2:9:63:0|1:30403_G_A:63,0,368 ./.:0,0:0:.:.:.:.
Chr01 30863 . C CTAT 111.58 . AC=2;AF=1;AN=2;DP=7;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1;MQ=60.00;QD=27.89;SOR=1.609 GT:AD:DP:GQ:PL ./.:0,0:0:.:. 1/1:0,4:4:12:146,12,0 ./.:0,0:0:.:. ./.:0,0:0:.:.
Chr01 32048 . T C 35.59 . AC=1;AF=0.167;AN=6;BaseQRankSum=-0.619;ClippingRankSum=-0.000;DP=44;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.167;MQ=60.00;MQRankSum=-0.000;QD=4.45;ReadPosRankSum=-1.345;SOR=0.307 GT:AD:DP:GQ:PGT:PID:PL ./.:0,0:0:.:.:.:. 0/0:25,0:25:72:.:.:0,72,1080 0/0:11,0:11:30:.:.:0,30,450 0/1:6,2:8:66:0|1:32048_T_C:66,0,243
Chr01 32072 . T C 35.59 . AC=1;AF=0.167;AN=6;BaseQRankSum=-0.619;ClippingRankSum=-0.000;DP=42;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.167;MQ=60.00;MQRankSum=-0.000;QD=4.45;ReadPosRankSum=0.319;SOR=0.693 GT:AD:DP:GQ:PGT:PID:PL ./.:0,0:0:.:.:.:. 0/0:23,0:23:66:.:.:0,66,990 0/0:11,0:11:33:.:.:0,33,263 0/1:6,2:8:66:0|1:32048_T_C:66,0,243
Chr01 33278 . G GA 37.34 . AC=1;AF=0.167;AN=6;BaseQRankSum=-1.520;ClippingRankSum=-0.000;DP=30;ExcessHet=4.7712;FS=3.136;MLEAC=2;MLEAF=0.333;MQ=60.00;MQRankSum=-0.000;QD=4.67;ReadPosRankSum=-0.272;SOR=1.447 GT:AD:DP:GQ:PL ./.:0,0:0:.:. 0/0:9,0:9:0:0,0,231 0/1:4,4:8:71:71,0,89 0/0:12,0:12:0:0,0,215
Chr01 36552 . C A 36.59 . AC=1;AF=0.167;AN=6;BaseQRankSum=-1.534;ClippingRankSum=-0.000;DP=29;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.167;MQ=60.00;MQRankSum=-0.000;QD=4.57;ReadPosRankSum=0.489;SOR=0.307 GT:AD:DP:GQ:PL ./.:0,0:0:.:. 0/0:9,0:9:27:0,27,346 0/0:12,0:12:36:0,36,455 0/1:5,3:8:67:67,0,173
Chr01 40038 . C CTATA 47.27 . AC=2;AF=0.5;AN=4;DP=11;ExcessHet=0.7918;FS=0.000;MLEAC=2;MLEAF=0.5;MQ=60.00;QD=23.64;SOR=0.693 GT:AD:DP:GQ:PL ./.:0,0:0:.:. 1/1:0,2:2:6:87,6,0 ./.:0,0:0:.:. 0/0:5,0:5:15:0,15,148
Chr01 40740 . A G 32.64 . AC=1;AF=0.167;AN=6;BaseQRankSum=-0.812;ClippingRankSum=-0.000;DP=27;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.167;MQ=60.00;MQRankSum=-0.000;QD=3.63;ReadPosRankSum=-1.221;SOR=0.495 GT:AD:DP:GQ:PGT:PID:PL ./.:0,0:0:.:.:.:. 0/0:12,0:12:30:.:.:0,30,450 0/1:7,2:9:63:0|1:40740_A_G:63,0,409 0/0:6,0:6:18:.:.:0,18,205
Chr01 40753 . T C 32.64 . AC=1;AF=0.167;AN=6;BaseQRankSum=-0.812;ClippingRankSum=-0.000;DP=26;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.167;MQ=60.00;MQRankSum=-0.000;QD=3.63;ReadPosRankSum=-0.000;SOR=0.495 GT:AD:DP:GQ:PGT:PID:PL ./.:0,0:0:.:.:.:. 0/0:11,0:11:33:.:.:0,33,382 0/1:7,2:9:63:0|1:40740_A_G:63,0,409 0/0:6,0:6:18:.:.:0,18,217
Chr01 43274 . A T 34.36 . AC=1;AF=0.25;AN=4;BaseQRankSum=-0.812;ClippingRankSum=-0.000;DP=18;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.25;MQ=60.00;MQRankSum=-0.000;QD=3.82;ReadPosRankSum=-0.431;SOR=0.223 GT:AD:DP:GQ:PGT:PID:PL ./.:0,0:0:.:.:.:. 0/0:9,0:9:27:.:.:0,27,358 0/1:7,2:9:63:0|1:43248_GT_G:63,0,285 ./.:0,0:0:.:.:.:.
"#"开头的都是注释说明行
##FORMAT=... --> FORMAT那一列对应的参数信息及解释，比如GT是字符串，表示基因型
##INFO=... --> INFO那一列对应的参数信息及解释，比如DP是整数，表示样本组合深度
##CommandLine... --> 执行命令行
##contig=... --> 染色体编号，长度信息
##reference=file:///path/to/ref.fa --> 参考基因组
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT P1 P2 S1 S2 --> 标题行
以序列编号开头的即为文件主体内容
每列含义与标题行一致
除了一下注释信息有所区别，其他列的注释信息与GCVF类似
第二列:物理位置（单个位点）
第八列：INFO列，里面具体参数解释在表头的注释行里。这里就没有END=了，所以VCF文件只包含有变异的位点
第十至最后一列：都是样本位点信息列，数据解释同GVCF。这里就凸显了bam文件的标签作用，样本对应名称是bam的标签。在合并GVCF之前要确定样本名称，不同样本名称不能一样，如：P1;P2;S1;S2。且要注意亲本子代对应关系。
关于GT和DP
GT，基因型，包括0/0，1/1，0/1，2/2，0/2等多种形式，分隔符除了'/'，还可以是'|'。0/0的意思是与参考基因组（即REF）一致的基因型，1/1的意思是与变异（即ALT）一致的基因型，0/1则是杂合的基因型，2指的是多向突变的碱基（ALT中存在类似于A,T这种）。"./."为缺失基因型。
DP，样本在这个位置的reads覆盖度。
注意
通过上面对GVCF文件和VCF文件的介绍，大家可能就会明白，为什么在做多样本变异信息文件合并的时候，必须得是用GVCF格式合并之后再转为VCF文件的格式了。
因为VCF包含的位点信息不完整，如果用VCF文件合并，就相当于在原始位点的基础上挑选出各自变异的位点，然后再取不同样本变异位点之间的交集，在这个过程中会丢失大量的计算位点。也不会存在亲本基因型一个为1/1，一个为0/0的情况了。而这个是BSA分析VCF文件过滤的主要步骤之一。
Bcftools同样也可以做变异检测，文件流程差不多，但是不同软件之间的格式会有些许的差别，这里不再详细介绍了，有兴趣的同学可以自行百度。
关于变异检测提速
相信很多同学都遇到过这个问题，非商业版的gatk速度慢的是可以，基因组稍大点就真的不知道要跑到猴年马月了。大家可以采取从bam那一步就按区间（不建议太小，一般按染色体就可以）拆分，然后并行去进行变异检测，这样效率会提高不少。
GVCF和VCF文件压缩和索引建立
软件安装下载
# 软件安装
wget https://sourceforge.net/projects/samtools/files/tabix/tabix-0.2.6.tar.bz2
tar xjvf tabix-0.2.6.tar.bz2
cd tabix-0.2.6/
make
#导入环境变量
echo "PATH=$PATH:~/softwarepath/tabix-0.2.6" >>~/.bashrc
source ~/.bashrc
# 软件使用
## 压缩VCF文件，可用gunzip解压，但压缩方式不一样
bgzip P1.vcf
## vcf文件建立索引
tabix -p vcf P1.vcf.gz
