本文中vcf文件的獲取主要是使用gatk的HaplotypeCaller模塊，該模塊的原理大概如下：
- 定義active regions 沿著參考基因組以一定的窗口滑動(dòng)蠢棱，統(tǒng)計(jì)比對(duì)的 mismatches, indels 和 softclips等信息計(jì)算基因組每個(gè)位置的活躍得分又厉，使用平滑算法進(jìn)行校摩，此處相當(dāng)于測(cè)定該區(qū)域熵值几睛。當(dāng)熵值達(dá)到某個(gè)設(shè)定的閾值是即確定該區(qū)域?yàn)閍ctive region败晴，用于后續(xù)組裝浓冒。
- 通過(guò)active regions的組裝確定單倍型 對(duì)于每個(gè)active regions，程序構(gòu)建一個(gè)類(lèi)似 De Bruijn 的圖來(lái)重新組裝active regions并識(shí)別 數(shù)據(jù)中可能存在哪些單倍型尖坤。然后稳懒，程序?qū)⒚總€(gè)單倍型與參考重新對(duì)齊 使用 Smith-Waterman 算法進(jìn)行單倍型分析，以識(shí)別潛在的變異位點(diǎn)慢味。
- 確定每個(gè)read的單倍型的似然值 對(duì)于每個(gè)active regions场梆，該程序使用PairHMM算法對(duì)每個(gè)單倍型的每個(gè)讀數(shù)進(jìn)行成對(duì)比對(duì)墅冷。這就產(chǎn)生了給定讀取數(shù)據(jù)的單倍型可能性矩陣。然后將這些可能性邊緣化辙谜，以獲得給定讀取數(shù)據(jù)的每個(gè)潛在變異位點(diǎn)的等位基因的可能性俺榆。
- 分配樣本基因型 對(duì)于每個(gè)潛在的變異位點(diǎn)，該程序應(yīng)用貝葉斯規(guī)則装哆，使用給定讀取數(shù)據(jù)的等位基因的可能性來(lái)計(jì)算給定該樣本的讀取數(shù)據(jù)的每個(gè)樣本的每個(gè)基因型的可能性罐脊。然后將最有可能的基因型分配給樣本。
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/
#導(dǎo)入環(huán)境變量
echo "PATH=$PATH:~/softwarepath/gatk-4.2.1.0" >>~/.bashrc
source ~/.bashrc
#同樣有懶人安裝方法
conda install -c bioconda gatk4
#測(cè)試是否安裝成功
gatk --help #檢查是否可正常運(yùn)行
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.
使用
因?yàn)檫@邊只是介紹該軟件在BSA分析中如何使用蜕琴，所以介紹相對(duì)單一萍桌，也是為了初學(xué)者不會(huì)被那么多冗余的信息干擾。普通變異檢測(cè)同樣適用凌简，當(dāng)樣本太多的時(shí)候可以試試gatk4的GenomicsDB模塊上炎，代替這里的CombineGVCFs + GenotypeGVCFs策略。
BSA分析需要通過(guò)該軟件獲得最終用于分析的vcf文件雏搂，所以這一步可以說(shuō)是入門(mén)的鑰匙藕施。
首先準(zhǔn)備輸入文件：混池?cái)?shù)據(jù)比對(duì)結(jié)果bam文件和參考基因組，比對(duì)過(guò)程參考前面的推文
這里文件命名分別為：
Ｍ怪！Ｉ咽场！一定要區(qū)分清楚樣本之間的對(duì)應(yīng)關(guān)系芙沥，不能混了
P1.sort.rmdup.bam 親本1混池比對(duì)結(jié)果文件
P2.sort.rmdup.bam 親本2混池比對(duì)結(jié)果文件
S1.sort.rmdup.bam 子代1混池比對(duì)結(jié)果文件（與親本1表型一致）
S2.sort.rmdup.bam 子代2混池比對(duì)結(jié)果文件（與親本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
##參數(shù)解讀
"-Xmx4g" 這個(gè)是指定給予的運(yùn)行內(nèi)存诲祸，可以比4g更大，但是要結(jié)合服務(wù)器本身的資源決定而昨。
-R 參考基因組
-I bam文件
-O 輸出文件
-ERC 輸出GVCF的格式
##其他樣本同樣如此操作救氯，分別得到對(duì)應(yīng)的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
##參數(shù)解讀
-V 不同樣本的gvcf文件
-O 合并后的gvcf文件
# GVCF轉(zhuǎn)VCF
gatk --java-options "-Xmx4g" GenotypeGVCFs \
-R ref.fa \
-V merge.g.vcf.gz \
-O merge.vcf.gz
##參數(shù)解讀
-V GVCF文件
-O VCF文件
結(jié)果文件格式解讀
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
"#"開(kāi)頭的都是注釋說(shuō)明行
##FORMAT=... --> FORMAT那一列對(duì)應(yīng)的參數(shù)信息及解釋?zhuān)热鏕T是字符串，表示基因型
##INFO=... --> INFO那一列對(duì)應(yīng)的參數(shù)信息及解釋?zhuān)热鏒P是整數(shù)歌憨，表示樣本組合深度
##CommandLine... --> 執(zhí)行命令行
##contig=... --> 染色體編號(hào)着憨，長(zhǎng)度信息
##reference=file:///path/to/ref.fa --> 參考基因組
##GVCFBlock0-1=minGQ=0(inclusive),maxGQ=1(exclusive) --> GVCF才有，VCF沒(méi)有
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT P1 --> 標(biāo)題行
以序列編號(hào)開(kāi)頭的即為文件主體內(nèi)容
每列含義與標(biāo)題行一致
第一列：序列編號(hào)
第二列:物理起始位置
第三列：變異位點(diǎn)的rsID號(hào)躺孝，如果沒(méi)有的話(huà)用"."表示享扔。
第四列：參考基因組在該位點(diǎn)的堿基，REF
第五列：樣本在該位點(diǎn)的堿基植袍，ALT
第六列：call出這個(gè)位點(diǎn)的質(zhì)量惧眠。這個(gè)值等于-10log10(p)，p值是call錯(cuò)alt allele錯(cuò)誤的概率于个。也就是QUAL越大出錯(cuò)概率越小氛魁。 第七列：對(duì)變異位點(diǎn)進(jìn)行過(guò)濾，如果通過(guò)則為PASS，如果沒(méi)有進(jìn)行過(guò)濾就是"."秀存。第八列 INFO：這一列是額外信息捶码。可能是平臺(tái)的信息或链，也可以是DP等的信息：
第八列：INFO列惫恼，里面具體參數(shù)解釋在表頭的注釋行里。這里使用第一行數(shù)據(jù)作為例子著重說(shuō)明一下澳盐，"END=5538"表示這一行表示的是從1到5538的變異信息祈纯。所以GVCF是包含所有位點(diǎn)的信息文件，不管位點(diǎn)是否存在變異叼耙。
第九列：FORMAT列腕窥，里面具體參數(shù)解釋在表頭的注釋行里
第十列：樣本位點(diǎn)信息列，數(shù)據(jù)以：為分隔符筛婉，數(shù)據(jù)含義與第九列保持一致簇爆。比如第九列是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:.:.:.:.
"#"開(kāi)頭的都是注釋說(shuō)明行
##FORMAT=... --> FORMAT那一列對(duì)應(yīng)的參數(shù)信息及解釋?zhuān)热鏕T是字符串首尼，表示基因型
##INFO=... --> INFO那一列對(duì)應(yīng)的參數(shù)信息及解釋?zhuān)热鏒P是整數(shù)，表示樣本組合深度
##CommandLine... --> 執(zhí)行命令行
##contig=... --> 染色體編號(hào)言秸，長(zhǎng)度信息
##reference=file:///path/to/ref.fa --> 參考基因組
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT P1 P2 S1 S2 --> 標(biāo)題行
以序列編號(hào)開(kāi)頭的即為文件主體內(nèi)容
每列含義與標(biāo)題行一致
除了一下注釋信息有所區(qū)別软能，其他列的注釋信息與GCVF類(lèi)似
第二列:物理位置（單個(gè)位點(diǎn)）
第八列：INFO列，里面具體參數(shù)解釋在表頭的注釋行里举畸。這里就沒(méi)有END=了查排，所以VCF文件只包含有變異的位點(diǎn)
第十至最后一列：都是樣本位點(diǎn)信息列，數(shù)據(jù)解釋同GVCF抄沮。這里就凸顯了bam文件的標(biāo)簽作用跋核，樣本對(duì)應(yīng)名稱(chēng)是bam的標(biāo)簽。在合并GVCF之前要確定樣本名稱(chēng)叛买，不同樣本名稱(chēng)不能一樣砂代，如：P1;P2;S1;S2。且要注意親本子代對(duì)應(yīng)關(guān)系率挣。
關(guān)于GT和DP
GT刻伊，基因型，包括0/0，1/1捶箱，0/1智什，2/2，0/2等多種形式丁屎，分隔符除了'/'荠锭，還可以是'|'。0/0的意思是與參考基因組（即REF）一致的基因型晨川，1/1的意思是與變異（即ALT）一致的基因型证九，0/1則是雜合的基因型，2指的是多向突變的堿基（ALT中存在類(lèi)似于A(yíng),T這種）础爬。"./."為缺失基因型甫贯。
DP，樣本在這個(gè)位置的reads覆蓋度看蚜。
注意
通過(guò)上面對(duì)GVCF文件和VCF文件的介紹叫搁，大家可能就會(huì)明白，為什么在做多樣本變異信息文件合并的時(shí)候供炎，必須得是用GVCF格式合并之后再轉(zhuǎn)為VCF文件的格式了渴逻。
因?yàn)閂CF包含的位點(diǎn)信息不完整，如果用VCF文件合并音诫，就相當(dāng)于在原始位點(diǎn)的基礎(chǔ)上挑選出各自變異的位點(diǎn)惨奕，然后再取不同樣本變異位點(diǎn)之間的交集，在這個(gè)過(guò)程中會(huì)丟失大量的計(jì)算位點(diǎn)竭钝。也不會(huì)存在親本基因型一個(gè)為1/1梨撞，一個(gè)為0/0的情況了。而這個(gè)是BSA分析VCF文件過(guò)濾的主要步驟之一香罐。
Bcftools同樣也可以做變異檢測(cè)卧波，文件流程差不多，但是不同軟件之間的格式會(huì)有些許的差別庇茫，這里不再詳細(xì)介紹了港粱，有興趣的同學(xué)可以自行百度。
關(guān)于變異檢測(cè)提速
相信很多同學(xué)都遇到過(guò)這個(gè)問(wèn)題旦签，非商業(yè)版的gatk速度慢的是可以查坪，基因組稍大點(diǎn)就真的不知道要跑到猴年馬月了。大家可以采取從bam那一步就按區(qū)間（不建議太小宁炫，一般按染色體就可以）拆分偿曙，然后并行去進(jìn)行變異檢測(cè)，這樣效率會(huì)提高不少淋淀。
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
#導(dǎo)入環(huán)境變量
echo "PATH=$PATH:~/softwarepath/tabix-0.2.6" >>~/.bashrc
source ~/.bashrc
# 軟件使用
## 壓縮VCF文件遥昧，可用gunzip解壓覆醇，但壓縮方式不一樣
bgzip P1.vcf
## vcf文件建立索引
tabix -p vcf P1.vcf.gz
