懒人先看：

* 标示的为最常用的参数。

$ java -Xmx20g -jar GenomeAnalysisTK.jar \  #使用GATK主程序
*    -R ref.fasta \  #参考序列
*    -T UnifiedGenotyper \  #使用GATK的该程序
*    -I sample1.bam [-I sample2.bam] [-I ...]\  #输入的bam比对结果
     --dbsnp dbSNP.vcf \  #使用该文件中的variants ID加入到结果文件中
     --genotyping_mode GENOTYPE_GIVEN_ALLELES --allels know.vcf \  #只对已知的variant进行基因分型
*    --genotype_likelihoods_model [SNP/INDEL/BOTH] \  #进行SNP、INDEl或者两者同时的calling
     --min_base_quality_score 17 \  #variant calling时碱基质量的最低要求。但是INDEL calling时，该值无效。而是去除reads的低质量3'端直到Q20。用下面的方法来鉴定
     --min_indel_count_for_genotyping 5 \  #至少有这么多的reads在某一位点表现出了indel
     --min_indel_fraction_per_sample 0.25 \  #至少有这么大比例的reads在某一点表现出了indel
*    --out gatk.cvf \  #输出到文件，否则为标准输出
*    --output_mode [EMIT_VARIANTS_ONLY/EMIT_ALL_CONFIDENT_SITES/EMIT_ALL_SITES] \  #输出模式,默认只输出variant位点
     --sample_ploidy 2 \  #样品的倍型
*    --standard_min_confidence_threshold_for_calling 30.0 \  #设定variant位点的置信阈值，低于该阈值则为low quality
*    --standard_min_confidence_threshold_for_emitting 30.0 \  #在vcf结果文件中，低于该值的位点则不进行报告

1. UnifiedGenotyper简介

UnifiedGenotyper是GATK（Genome Analysis ToolKit)中一个主要工具，用于Variant calling。在GATK网站上这样描述它：A variant caller which unifies the approaches of several disparate callers — Works for single-sample and multi-sample data.

UnifiedGenotyper能对单个或多个sample进行SNP和INDEL calling。使用Beyesian genotype likelihood model来对N个sample进行基因型的判断和等位基因频率的计算。

2. 命令的简单使用和输入要求

一个使用UnifiedGenotyper对多个samples进行SNP calling的例子如下：

$ java -jar GenomeAnalysisTK.jar \
    -R resources/Homo_sapiens_assembly18.fasta \
    -T UnifiedGenotyper \
    -I sample1.bam [-I sample2.bam] \
    --dbsnp dbSNP.vcf \
    -o snps.raw.vcf \
    -stand_call_conf [50.0] \
    -stand_emit_conf 10.0 \
    -dcov [50 for 4x, 200 for >30x WGS or Whole exome] \
    [-L targets.interval_list]

该上述命令中输入的文件有2个：参考序列的fasta文件 和 bam格式的比对结果文件。但是只有这两个文件是不行的，其实额外需要fasta文件的dict文件和fai文件，以及bam文件的bai文件。需要使用picard和samtools软件，命令行如下：

$ java -jar $picardHome/CreateSequenceDictionary.jar R=./resources/Homo_sapiens_assembly18.fasta O=./resources/Homo_sapiens_assembly18.dict
$ samtools faidx ./resources/Homo_sapiens_assembly18.fasta
$ samtools index sample1.bam

3. 附加信息

3.1 Read filters

在运行UnifiedGenotyper的时候，会对reads自动进行过滤，使用的相关命令有：
DuplicateReadFilter
FailsVendorQualityCheckFilter
NotPrimaryAlignmentFilter
MalformedReadFilter
BadMateFilter
MappingQualityUnavailableFilter
UnmappedReadFilter

3.2 并行化运算

该程序能进行多线程运算，使用如下参数即可：
NanoSchedulable (-nct)
CPU数。对每个data的运行能使用的CPU数。即对一个数据执行的一个命令进行计算所能用的CPU数。
TreeReducible (-nt)
线程数。将一个总的data分成多少份，然后分别对每个data使用单独的命令进行运算，最大值为24. 最后合并结果可能消耗额外的时间；消耗的内存也按倍数增加。

3.3 Downsampling 设置

This tool overrides the engine’s default downsampling settings.

Mode: BY_SAMPLE
To coverage: 250

3.4 Window size

This tool uses a sliding window on the reference.

Window start: -200 bp before the locus
Window stop: 200 bp after the locus

4. UnifiedGenotyper的参数

4.1 GATK参数

使用了GATK共有的参数:CommandLineGATK.

4.2 UnifiedGenotyper特有的参数

• --alleles / -alleles RodBinding[VariantContext] default: none

The set of alleles at which to genotype when –genotyping_mode is GENOTYPE_GIVEN_ALLELES. When the UnifiedGenotyper is put into GENOTYPE_GIVEN_ALLELES mode it will genotype the samples using only the alleles provide in this rod binding –alleles binds reference ordered data. This argument supports ROD files of the following types: BCF2, VCF, VCF3

• --annotateNDA / -nda default: false

If provided, we will annotate records with the number of alternate alleles that were discovered (but not necessarily genotyped) at a given site. Depending on the value of the –max_alternate_alleles argument, we may genotype only a fraction of the alleles being sent on for genotyping. Using this argument instructs the genotyper to annotate (in the INFO field) the number of alternate alleles that were originally discovered at the site.

• --annotation / -A List[String] default: none

One or more specific annotations to apply to variant calls. Which annotations to add to the output VCF file. See the VariantAnnotator -list argument to view available annotations.

• --comp / -comp List[RodBinding[VariantContext]] default: none

Comparison VCF file. If a call overlaps with a record from the provided comp track, the INFO field will be annotated as such in the output with the track name (e.g. -comp:FOO will have ‘FOO’ in the INFO field). Records that are filtered in the comp track will be ignored. Note that ‘dbSNP’ has been special-cased (see the –dbsnp argument). –comp binds reference ordered data. This argument supports ROD files of the following types: BCF2, VCF, VCF3

• --computeSLOD / -slod boolean default: false

If provided, we will calculate the SLOD (SB annotation). Note that calculating the SLOD increases the runtime by an appreciable amount.

• -contamination / --contamination_fraction_to_filter double default: 0.05

Fraction of contamination in sequencing data (for all samples) to aggressively remove. If this fraction is greater is than zero, the caller will aggressively attempt to remove contamination through biased down-sampling of reads. Basically, it will ignore the contamination fraction of reads for each alternate allele. So if the pileup contains N total bases, then we will try to remove (N * contamination fraction) bases for each alternate allele.

• -contaminationFile / --contamination_fraction_per_sample_file File

Tab-separated File containing fraction of contamination in sequencing data (per sample) to aggressively remove. Format should be “” (Contamination is double) per line; No header.. This argument specifies a file with two columns “sample” and “contamination” specifying the contamination level for those samples. Samples that do not appear in this file will be processed with CONTAMINATION_FRACTION

• --dbsnp / -D RodBinding[VariantContext] default: none

dbSNP file. rsIDs from this file are used to populate the ID column of the output. Also, the DB INFO flag will be set when appropriate. dbSNP is not used in any way for the calculations themselves. –dbsnp binds reference ordered data. This argument supports ROD files of the following types: BCF2, VCF, VCF3

• --excludeAnnotation / -XA List[String] default: none

One or more specific annotations to exclude. Which annotations to exclude from output in the VCF file. Note that this argument has higher priority than the -A or -G arguments, so annotations will be excluded even if they are explicitly included with the other options.

• --genotype_likelihoods_model / -glm Model default: SNP

Genotype likelihoods calculation model to employ — SNP is the default option, while INDEL is also available for calling indels and BOTH is available for calling both together.
The –genotype_likelihoods_model argument is an enumerated type (Model), which can have one of the following values:
SNP
INDEL
GENERALPLOIDYSNP
GENERALPLOIDYINDEL
BOTH

• --genotyping_mode / -gt_mode GENOTYPING_MODE default: DISCOVERY

Specifies how to determine the alternate alleles to use for genotyping.
The –genotyping_mode argument is an enumerated type (GENOTYPING_MODE), which can have one of the following values:
DISCOVERY
the Unified Genotyper will choose the most likely alternate allele
GENOTYPE_GIVEN_ALLELES
only the alleles passed in from a VCF rod bound to the -alleles argument will be used for genotyping

• --group / -G String[] default: [Standard]

One or more classes/groups of annotations to apply to variant calls. If specified, all available annotations in the group will be applied. See the VariantAnnotator -list argument to view available groups. Keep in mind that RODRequiringAnnotations are not intended to be used as a group, because they require specific ROD inputs.

• --heterozygosity / -hets Double default: 0.0010

Heterozygosity value used to compute prior likelihoods for any locus. The expected heterozygosity value used to compute prior likelihoods for any locus. The default priors are: het = 1e-3, P(hom-ref genotype) = 1 – 3 * het / 2, P(het genotype) = het, P(hom-var genotype) = het / 2

• --indel_heterozygosity / -indelHeterozygosity double default: 1.25E-4

Heterozygosity for indel calling. This argument informs the prior probability of having an indel at a site.

• --indelGapContinuationPenalty / -indelGCP byte default: 10

Indel gap continuation penalty, as Phred-scaled probability. I.e., 30 => 10^-30/10.

• --indelGapOpenPenalty / -indelGOP byte default: 45

Indel gap open penalty, as Phred-scaled probability. I.e., 30 => 10^-30/10.

• --input_prior / -inputPrior List[Double] default: none

Input prior for calls. By default, the prior specified with the argument –heterozygosity/-hets is used for variant discovery at a particular locus, using an infinite sites model, see e.g. Waterson (1975) or Tajima (1996). This model asserts that the probability of having a population of k variant sites in N chromosomes is proportional to theta/k, for 1=1:N There are instances where using this prior might not be desireable, e.g. for population studies where prior might not be appropriate, as for example when the ancestral status of the reference allele is not known. By using this argument, user can manually specify priors to be used for calling as a vector for doubles, with the following restriciotns: a) User must specify 2N values, where N is the number of samples. b) Only diploid calls supported. c) Probability values are specified in double format, in linear space. d) No negative values allowed. e) Values will be added and Pr(AC=0) will be 1-sum, so that they sum up to one. f) If user-defined values add to more than one, an error will be produced. If user wants completely flat priors, then user should specify the same value (=1/(2*N+1)) 2*N times,e.g. -inputPrior 0.33 -inputPrior 0.33 for the single-sample diploid case.

• --max_alternate_alleles / -maxAltAlleles int default: 6

Maximum number of alternate alleles to genotype. If there are more than this number of alternate alleles presented to the genotyper (either through discovery or GENOTYPE_GIVEN ALLELES), then only this many alleles will be used. Note that genotyping sites with many alternate alleles is both CPU and memory intensive and it scales exponentially based on the number of alternate alleles. Unless there is a good reason to change the default value, we highly recommend that you not play around with this parameter. As of GATK 2.2 the genotyper can handle a very large number of events, so the default maximum has been increased to 6.

• --max_deletion_fraction / -deletions Double default: 0.05

Maximum fraction of reads with deletions spanning this locus for it to be callable [to disable, set to < 0 or > 1; default:0.05].

• --min_base_quality_score / -mbq int default: 17

Minimum base quality required to consider a base for calling. The minimum confidence needed in a given base for it to be used in variant calling. Note that the base quality of a base is capped by the mapping quality so that bases on reads with low mapping quality may get filtered out depending on this value. Note too that this argument is ignored in indel calling. In indel calling, low-quality ends of reads are clipped off (with fixed threshold of Q20).

• -minIndelCnt / --min_indel_count_for_genotyping int default: 5

Minimum number of consensus indels required to trigger genotyping run. A candidate indel is genotyped (and potentially called) if there are this number of reads with a consensus indel at a site. Decreasing this value will increase sensitivity but at the cost of larger calling time and a larger number of false positives.

• -minIndelFrac / --min_indel_fraction_per_sample double default: 0.25

Minimum fraction of all reads at a locus that must contain an indel (of any allele) for that sample to contribute to the indel count for alleles. Complementary argument to minIndelCnt. Only samples with at least this fraction of indel-containing reads will contribute to counting and overcoming the threshold minIndelCnt. This parameter ensures that in deep data you don’t end up summing lots of super rare errors up to overcome the 5 read default threshold. Should work equally well for low-coverage and high-coverage samples, as low coverage samples with any indel containing reads should easily over come this threshold.

• --out / -o VariantContextWriter default: stdout

File to which variants should be written. A raw, unfiltered, highly sensitive callset in VCF format.

• --output_mode / -out_mode OUTPUT_MODE default: EMIT_VARIANTS_ONLY

Specifies which type of calls we should output.
The –output_mode argument is an enumerated type (OUTPUT_MODE), which can have one of the following values:
EMIT_VARIANTS_ONLY
produces calls only at variant sites
EMIT_ALL_CONFIDENT_SITES
produces calls at variant sites and confident reference sites
EMIT_ALL_SITES
produces calls at any callable site regardless of confidence; this argument is intended only for point mutations (SNPs) in DISCOVERY mode or generally when running in GENOTYPE_GIVEN_ALLELES mode; it will by no means produce a comprehensive set of indels in DISCOVERY mode

• --pair_hmm_implementation / -pairHMM HMM_IMPLEMENTATION default: ORIGINAL

The PairHMM implementation to use for -glm INDEL genotype likelihood calculations. The PairHMM implementation to use for -glm INDEL genotype likelihood calculations. The various implementations balance a tradeoff of accuracy and runtime.
The –pair_hmm_implementation argument is an enumerated type (HMM_IMPLEMENTATION), which can have one of the following values:
EXACT
ORIGINAL
LOGLESS_CACHING

• --pcr_error_rate / -pcr_error Double default: 1.0E-4

The PCR error rate to be used for computing fragment-based likelihoods. The PCR error rate is independent of the sequencing error rate, which is necessary because we cannot necessarily distinguish between PCR errors vs. sequencing errors. The practical implication for this value is that it effectively acts as a cap on the base qualities.

• --sample_ploidy / -ploidy int default: 2

Plody (number of chromosomes) per sample. For pooled data, set to (Number of samples in each pool * Sample Ploidy)..

• -stand_call_conf / --standard_min_confidence_threshold_for_calling double default: 30.0

The minimum phred-scaled confidence threshold at which variants should be called. The minimum phred-scaled Qscore threshold to separate high confidence from low confidence calls. Only genotypes with confidence >= this threshold are emitted as called sites. A reasonable threshold is 30 for high-pass calling (this is the default).

• -stand_emit_conf / --standard_min_confidence_threshold_for_emitting double default: 30.0

The minimum phred-scaled confidence threshold at which variants should be emitted (and filtered with LowQual if less than the calling threshold). This argument allows you to emit low quality calls as filtered records.