基因组变异检测(Variance Calling with GATK)

原始来源，仅供备忘：
https://xizhihui.github.io/2018/10/14/%E5%9F%BA%E5%9B%A0%E7%BB%84%E5%8F%98%E5%BC%82%E6%A3%80%E6%B5%8B(Variance-Calling-with-GATK)/

一、基本概念

1.1 名词解释

• 基因组：个体全部DNA序列的无重复集.这里的基因组不仅仅包含了基因在内,由于目前尚有许多DNA序列不编码蛋白,也可能不会转录,反正就是这些序列的功能还没有研究清楚, 这些序列也都包含在基因组这个范畴里面.
• Reads：二代测序中的一个专有名词,表示着测序仪对某个DNA片段的一次测序结果,是该DNA序列的序列组成. 其长度依据测序仪不同而不同.
• 变异：variants, 变异是一个相对的概念,产生于比较之中, 比较是指同耳熟能详的参考基因组相比较. 对于人类基因组的变异来讲,参考基因组是经过“人类基因组”计划测序所得到的最终人类基因组序列.

1.2 基因组变异的种类

个体基因上的变化可以分为两种, 一种是种系上的(germline), 这种变化由遗传所得, 所以还有phasing分析，以区别这个变化来源于哪个亲本。 另一种是体细胞(somatic)上的, 这通常是由癌症所引起的异常变化, 也就是癌变相关的变化。

• 单核苷酸多态性(single nucleotide polymorphism, SNP)
  单核苷酸多态性主要是指在基因组水平上由单个核苷酸的变异所引起的DNA序列多态 性。它是人类可遗传的变异中最常见的一种，占所有已知多态性的80%以上。SNP在人类 基因组中广泛存在，平均每500?1000个碱基对中就有1个，估计其总数可达300万个甚至更多。SNP所表现的多态性只涉及到单个碱基的变异，这种变异可由单个碱基的转换(transition)或颠换(transversion)所引起，也可由碱基的插入或缺失所致。但通常所说的SNP并不包括后两种情况。

• 插入和缺失(insertion and deletion, INDEL)
  插入是指在基因组上的某个位置出现额外的DNA片段, 而缺失则是指某些位置原有的DNA片段发生丢失. 插入和缺失的DNA片段长度一般不是太长, 基本在50bp以内.

• 染色体结构变异(chromosome structure variant)
  当DNA片段缺失的长度大于50bp时,我们则认为发生了染色体结构缺失. 不仅如此, 染色体结构变异还包含其他结构变化, 如染色体的易位(两条染色体)、倒位(自身).有时候, 色体数目变异也被归为是染色体结构变异.

1.3 基因组结构变异的检测策略

• Pair-end Mapping(also named read pair, PEM)
  PEM法可以检测到许多变异, 如deletion, insertion, inversion, intra/inter-chromosome translocation, tandem duplication, interspersed duplications.对应的软件有Breakdancer, variationHunter, Spanner, PEMer.注意, PEM法深受insertSize影响.

• Split read(SR)
  SR法可以检测deletion, insertion, tandem duplication, inversion. SR法需要软件具备soft-clip reads的能力.相关软件有Pindel.

• Read Depth(RD)
  RD检测一些大的deletion或者duplication事件,对于小的变异则无能为力.RD有两种策略,一种是通过检测样本在一个参考基因组上read的深度分布情况来检测CNV,适用于单样本;另一种则是通过和识别出比较两个样本中所存在的丢失和重复倍增区, 以此来获得相对的CNV, 适用于case-control模型的样本.相关的软件有CNVnator, CNV-seq.

• De novo assembly
  从头组装对于变异检测最为有效,特别是long insertion和复杂结构性变异.

1.4 各个检测策略图示

https://s1.ax1x.com/2018/10/14/iU8cPH.png

二、基因组变异检测软件

• bwa: 用于序列比对
• samtools: 处理比对后数据
• picard: 处理比对后数据,功能上与samtools互补
• GATK: variant calling

三、数据预处理：获得可进行分析的bam文件

3.1 比对到参考基因组上,并排序

可用软件：BWA (DNA-seq), Picards MergeBamAlignments(针对unmapped BAM文件), STAR(RNA-seq). 注意需要排序bam文件.

# 1. 下载参考基因组数据 2. bwa建立索引
bwa index ref_genome.fasta # output: ref_genome.amb, .ann, .bwt, .pac, .sa
# 3. bwa比对,使用mem模块
bwa mem -t 4 -R '@RG\tID:foo_lane\tPL:illumina\tLB:library\tSM:sample_name' /path/to/ref_genome.fasta read_1.fq.gz read_2.fq.gz > sample_name.sam
    # -R: 用'\t'分割的Read Group的字符串信息,为read指定分组. 
    # ID是分组ID,设置为测序的laneID;PL是平台名字;SM是样本ID,用于区分样本.LB是测序文库名字.
# 4. 排序
samtools view -S -b sample_name.sam > sample_name.bam
samtools sort -@ 4 -m 4G -O bam -o sample_name.sorted.bam sample_name.bam

3.2 去除重复reads

重复reads是指来自同个DNA模版的reads,其在参考序列上的起始点相同; 对于PE测序, 两个端点的起始点都相同。在测序的过程中，序列是独立的，单个DNA模版的产物是相同，变异检测基于这种DNA模版产物相同而进行；在样品处理、建库测序中产生的错误(比如PCR bias，flowcell lane的邻近cluster污染)将会集中到重复，重复会增大变异检测结果的假阴性和假阳性。
一般来讲，扩增子测序和RNA-seq进行基因表达分析不用去除重复。

# 1.使用samtools去重复
samtools rmdup sample_name.sorted.bam sample_name.sorted.rmdup.bam
# 1.使用picard去重复
java -jar picard.jar MarkDuplicates \ 
 REMOVE_DUPLICATES=true \              # 如果不希望删除重复序列,只是进行标注,本语句不加
 I=sample_name.sorted.bam \
 O=sample_name.sorted.rmdup.bam \
 M=sample_name.markdup_metrics.txt
# 2\. 建立samtools index,方便访问
samtools index sample_name.sorted.rmdup.bam # sample_name.sorted.rmdup.bam.bai

3.3 局部重比对(locate region realignment)

本步骤会去除由多连续碱基引起的INDELs, 对于基于组装的变异检测此步不是必需,因为组装过程中多连续碱基引起的变异检测会被忽略(对于RNA-seq,会进行额外的Split ‘N’ Trim)。本步的目的在于对潜在的序列插入或删除的区域进行校正;因为全局搜索最优匹配算法在存在插入删除的区域附近的比对结果不是很好;再者,不同算法对错配和gap的容忍度不同,同样会导致不同的比对结果。本步执行的效果可见图：https://s1.ax1x.com/2018/07/04/PEHlAx.png。

|

# 1\. 定位需要重比对的目标区域
java -jar /path/to/GenomeAnalysisTK.jar \
 -T RealignerTargetCreator \
 -R /path/to/ref_name.fasta \
 -I sample_name.sorted.rmdup.bam \     # 有多个样本就有多个本语句
 -known /path/to/gatk/bundle/1000G_phase1.indels.b37.vcf \       # 添加已知的候选variant区域;下同; b37的版本对应这ref_name.fasta的版本
 -known /path/to/gatk/bundle/Mills_and_1000G_gold_standard.indels.b37.vcf \
 -o sample_name.IndelRealigner.intervals
# 2\. 进行序列重比对
java -jar /path/to/GenomeAnalysisTK.jar \
 -T IndelRealigner \
 -R /path/to/ref_name.fasta \
 -I sample_name.sorted.markdup.bam \
 -known /path/to/gatk/bundle/1000G_phase1.indels.b37.vcf \
 -known /path/to/gatk/bundle/Mills_and_1000G_gold_standard.indels.b37.vcf \
 -o sample_name.sorted.rmdup.realign.bam \
 --targetIntervals sample_name.IndelRealigner.intervals

|

3.4 碱基质量校正(Base Quality Score Recalibration, BQSR)

系统性误差是变异检测结果不好的主要原因：不同测序技术和测序仪器的系统误差都不一样。再者，后续的分析是基于碱基质量，所以碱基质量分数不能含有PCR bias引起的偏差。为了(尽可能)校正这些错误,BQSR实际上是在执行的时候是按照不同的测序lane或者测序文库来进行，这个时候@RG信息(BWA比对时所设置的)就显得很重要了,算法就是通过@RG中的ID来识别各个独立的测序过程，其过程简要概括为计算lane, original quality score, machine cycle, sequencing context的covariates。

# 1\. build base recalibration model
gatk BaseRecalibrator \
 -R ref.fasta -I sample.bam \
 -knownSites snps.vcf.gz \
 -knownSites indels.vcf.gz \
 -O recal.table
 #-bqsr 1st_recal.table #如果是进行第二次的recallibrator
# 2\. ApplyBQSR
gatk ApplyBQSR \
 -R ref.fasta \
 -I sample.bam \
 -bqsr recal.table \
 -O sample_bqsr.bam
 # parameters below are in GATK3, GATK4 ApplyBQSR omits by default
 # -SQQ 10 -SQQ 20 -SQQ 30 -SQQ 40  # to bin quals
 # --emit-original_quals            # to emit original quals to OQ tag
 # --preserve_qscores_less_than 6   # BQs less than 6 are untouched
# 3\. 可视化校正效果,需要运行1\. build base recalibration model第二次再进行可视化
gatk AnalyzeCovariates -before 1st_recal.table \
 -after 2nd_recal.table -plots plots.pdf

四、Germline variant discovery：SNPs & Indels

4.1 单倍体检测(HaplotypeCaller)

HaplotypeCaller有2个模式可选，一个是默认模式，所有bam文件一起跑；一个GVCF模式，bam文件分开跑，生成中间g.vcf文件，这个模式利于后面添加新的样品联合分析，即所谓的N+1。

# basic mode, no GVCF
gatk HaplotypeCaller \
 -R reference.fasta \
 -I preprocessed_reads.bam -O germline_variants.vcf
# GVCF mode
gatk HaplotypeCaller \
 -R reference.fasta -I preprocessed_reads.bam \
 -O germline_variants.g.vcf -ERC GVCF

4.1.1 HaplottypeCaller的运行逻辑

1. identify activeRegions
   sliding window along the reference => count mismatches, indels and soft-clips => measure of entropy => trim and continue with activeRegions over threshold

2. assemble plusible haplotypes
   local realignment via graph assembly => traverse graph to collect most likely haplotypes => align haplotypes to reference using Smith-Waterman => get likely haplotypes and candidate variant sites

3. score haplotypes using PairHMM
   PairHMM aligns each read to each haplotype => uses base qualities as the estimate of error => get likelihoods of the haplotypes given reads

4. genotype each sample at each potential variant site
   detemine most likely combination of alleles for each site => based on allele likelihoods (from PairHMM) => apply bayes’theorem with ploidy assumption

4.2 合并calling (joint calling, necessary for GVCF mode)

可以进行联合分析，如前所述，有新的样品加入；如果是人的或已经有人分析的同类样品数据，可以加入进来进行联合分析

# GATK3, with CombineGVCFs
gatk CombineGVCFs \
    -R reference.fasta \
    -V sample1.g.vcf \
    -V sample2.g.vcf \
    -O combined.g.vcf

# GATK4, with GenomicsDBImport
gatk GenomicsDBImport \
    -R reference.fasta \
    -V sample1.g.vcf \
    -V sample2.g.vcf \
    -L chr20 \
    --genomicsdb-workspace-path gvcfs_db

# on a single- or multi-sample GVCF, can be the combinned.g.vcf
gatk GenotypeGVCFs \
    -R reference.fasta \
    -V variants.g.vcf \
    -O final_variants.vcf
    
# on a genomicsDB workspace
gatk GenotypeGVCFs \
    -R reference.fasta \
    -V gendb://gvcfs_db \
    -O final_variants.vcf

应该对原始的变异检测集进行过滤,以平衡检测集的敏感性和特异性.由于变异检测的算法通常是十分宽松开放的, 变异检测敏感性升高的同时自然会带来许多假阳性.这时候,我们有两种策略来进行过滤. 其一是使用二分阈值进行硬过滤,十分简单粗暴 但这需要我们有足够的经验去确定这个阈值. 其二则是使用机器学习的方法进行变异校正(variant recalibration，VR),对变异的辨识度较好, 但这个方法需要我们有足够多的已知数据进行训练.采用已知的高置信度的变异位点数据集进行训练, 得到某位点是否为变异位点的概率. VQSLOD就是采用此类方法进行过滤的.不过,可以肯定的是,这两种策略都应该协调好敏感性和特异性的问题,也需要对变异环境进行注释(use variant context annotations)。

VR的处理简要过程为：假定变异位点注释符合高斯聚类，然后基于已知变异位点注释数据集建立高斯混合模型(Gaussian mixure model)，接着根据位点与聚类之间的距离对所有位点进行打分，最后基于打分进行过滤。

当你的样品不是人类或者没有足够的已知数据集时,不用进行本步骤. RNA-seq有独有的过滤操作,也不用进行本步骤,样本数太少也不用进行本步骤。 SNPs和Indels要分开各自进行校正,也就是说要以SNPmode和Indelmode各校正一次.

对于外显子数据,样本量太少的话,数据不足以建立一个健壮的高斯模型,联合分析(jointly)需要至少30个样本.

4.4 表型加强(Genotype refinement)

通过使用额外的数据集改善表型检测和似然情况(genotype calls & likelihoods)。具体过程如下：

# 1.
gatk CalculateGenotypePosteriors \
    -R reference.fasta -V input.vcf \
    -ped family.ped -supporting population.vcf \
    -O output.vcf

# 2. filter low confidence GQs
gatk VariantFiltration \
    -R reference.fasta -V input.vcf \
    --filter-expression 'GQ<20' --filter-name 'lowGQ' \
    -O output.vcf
    
# 3. variantAnnotator
gatk VariantAnnotator \
    -R reference.fasta -V input.vcf \
    -A PossibleDenovo -O output.vcf

4.5 变异集验证(callset evaluation)

callset evaluation是为了验证variant calling结果的好坏，false positive怎么样。这里我们假设基于如果calling set是真实的，那么这个结果是典型且可比较的, 不一致性则表明有错误。(guiding principle: divergence is indicative of error
key assumption: truth set is representative or comparable)。通常用于作为true set的有commonly used truth sets: dbSNP, ExAC and GnomAD, HapMap, OMNI, NIST’s Genomes in a Bottle, Illumina’s Platinum Genomes。一般来讲，TiTv Ratio(Transitions / Transversions)在全基因组测序中的值为2.0-2.1，在全外显子测序中的值为3.0-3.3。

# 1.GATK: get Variant level evaluation
java -jar GenomeAnalysisTK.jar \
 -T VariantEval -R reference.b37.fasta \
 -eval callset.vcf --comp truthset.vcf \
 -o results.eval.grp

# 2.GATK: get Genotype concordance
java -jar GenomeAnalysisTK.jar \
 -T GenotypeConcordance --comp truthset.vcf \
 --eval callset.vcf -o result.grp

# 3.Picard: get variant level evalution
java -jar picard.jar \
 CollectVariantCallingMetrics INPUT=callset.vcf \
 DBSNP=truthset.vcf OUTPUT=results

# 4.Picard: get genotype concordance
java -jar picard.jar \
 GenotypeConcordance CALL_VCF=callset.vcf \
 TRUTH_VCF=truthset.vcf call_SAMPLE=sampleName \
 TRUTH_SAMPLE=sampleName output=results

五、Somatic variant discovery: SNPs & Indels

5.1 用Mutect2检测体细胞小变异,产生bamout文件

|

# 这里在tumor sample和对应的normal sample(a panel of normals, PoN)寻找变异
# 此外,这里还用到了群体种系变异数据集(a population germline variant resource)
gatk --java-options '-Xmx2g' Mutect2 \
 -R hg38/Homo_sapiens_assembly38.fasta \
 -I tumor.bam \
 -I normal.bam \
 -tumor HCC1143_tumor \
 -normal HCC1143_normal \
 -pon resources/chr17_pon.vcf.gz \
 --germline-resource resources/chr17_af-only-gnomad_grch38.vcf.gz \
 --af-of-alleles-not-in-resource 0.0000025 \
 --disable-read-filter MateOnSameConigOrNoMapptedMateReadFilter \
 -L chr17plus.interval_list \
 -O 1_somatic_m2.vcf.gz \
 -bamout 2_tumor_normal_m2.bam

# parameters(下页还有参数)
#-I      提供bam文件,这里分别是肿瘤和正常的bam文件
#-tumor  提供bam文件对应的样本名(SM值),下同
#-normal 同上.可通过samtools view -H tumor.bam | grep '@RG'语句获得
#-pon
#预过滤正常变异位点集(a panel of normals callset)内的变异位点,需要指定PoN VCF文件.PoN不仅包含了常见的种系变异位点，也包含了常见的测序噪音.默认情况下,工具不会对与PoN中匹配的变异位点进行重组装或者判断为变异位点.如果要允许对PoN位点表型分型,使用--genotype-pon-sites设定.

#--germline-resource
#通过指定群体种系资源(population germline resource)来对变异等位基因做注释群体种系资源必需包含等位基因特异频率,如在INFO域的AF注释.本工具使用种群等位基因频率对变异等位基因进行注释.当使用本参数时,需要考虑对--af-of-alleles-not-in-resource参数进行调整到默认的0.001.例如gnomAD resouce af-only-gnomad_grch38.vcf.gz含有200K的外显子和16k的基因组;而-I指定的数据是外显子数据,所以调整值是1/(2*exome)*0.001,即1/(2*200)*0.001=0.0000025.默认的0.001是没有种群资源的人类样品的大概值,这是基于人类平均杂合率计算得到.种群等位基因频率(POP_AF)和af-of-alleles-not-in-resource在计算变异位点时的概率有用到.
#--disable-read-filter MateOnSameConigOrNoMapptedMateReadFilter
#该参数指定保留那些配对read比对到其他不同contig的reads(pair end reads成对的).该参数的效果保留的reads均是比对到alternate contigs和span contigs的reads,对变异检测有帮助作用.
#-L      指定特定的目标基因组区域进行分析.这里指定该参数的目的是为了减少运行时间,加快速度.
#-bamout 该参数指定生成重组装的比对文件.它包含手工的单倍体和normal同Tumor的重组装比对文件.这使得手工检阅获得的变异检测称为可能.本参数可以忽略,但推荐使用,毕竟不会额外花多少时间.

# output => 1_somatic_m2.vcf.gz(raw unfiltered somatic callset), 2_tumor_normal_m2.bam (a reassembled reads bam), 1_somatic_m2.vcf.gz.tbi (indices), 2_tumor_normal_m2.bai (indices)

|

5.2 创建一个只有位点(sites-only)的PoN

Mutect2的tumor-only模式不止可以用于创建PoN,还可以用于样本数很少的变异检测,如线粒体DNA,也可以用于比较2个不同样本之间的变异差异.

|

# 这里面我们使用三个种系样品进行创建, HG00190,NA19771,HG02759
# 1\. 对每个样品运行tumor-only模式的Mutect2, 以HG00190为例
gatk Mutect2 \
 -R path/hg38/Homo_sapiens_assembly38.fasta \
 -I HG00190.bam \
 -tumor HG00190 \
 --disable-read-filter MateOnSameContigOrNoMappedMateReadFilter \
 -L chr17plus.interval_list \
 -O 3_HG00190.vcf.gz
# output => 3_HG00190.vcf.gz, 3_HG00190.bai

# 2\. 使用CreateSomaticPanelOfNormals对所有的normal VCFs进行核对,并合并成一个变异检测集,该方法保留在所有样品vcf中都存在的变异位点
# 注意,这里我们只使用3个样品,正常来讲,应该最少使用40个样本
gatk CreateSomaticPanelOfNormals \
 -vcfs 3_HG00190.vcf.gz \
 -vcfs 4_NA19771.vcf.gz \
 -vcfs 5_HG02759.vcf.gz \
 -O 6_threesamplepon.vcf.gz

#6_threesamplepon_vcf.gz:
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO
#chr6 29942399   .   C   T   .   .   .

5.3 使用GetPileupSummaries和CalculateContamination检测是否存在跨样本污染

# 1\. 对tumor BAM运行GetPileupSummaries得到reads概览,为已知的-系列变异位点做支持
gatk GetPileupSummaries \
 -I tumor.bam \
 -V resources/chr17_small_exac_common_3_grch38.vcf.gz \
 -O 7_tumor_getpileupsummaries.table

# output => 7_tumor_getpileupsummaries.table:
# contig  position    ref_count   alt_count   other_alt_count allele_frequency

#parameters
#--minimum-population-allele-frequency   指定最小的种群等位基因频率,默认为0.01
#--maximum-population-allele-frequency
#默认为0.2.这个工具只考虑等位基因频率在这两个范围内的纯合子可变位点.这样考虑的原因在于,如果纯合子可变位点有罕见的等位基因的话,我们更倾向于将其作为REF.那如果有跨样品污染的话,我们更倾向于将其作为平常的等位基因.
#-L  指定该参数只是为了减少运行时间

# 2\. 使用CalculateContamination检测污染
# 这里有两个模式,简单模式如下. 
# 匹配模式需要对normal sample运行GetPileupSummaries,然后把输出文件在这里指定为-matched参数值
gatk CalculateContamination \
 -I 7_tumor_getpileupsummaries.table \
 -O 8_tumor_calculatecontamination.table

|

5.4 使用FilterMutectCalls过滤出确信的体细胞变异

FilterMutectCalls可确定一个体细胞变异的置信度, 采用变异检测集里面的注释,使用预设的阈值来确定.

# 这里使用全部的变异检测集
# 为了在污染检测的基础上有效进行过滤,可以添加--contamination-table指定污染表(contamination-table)
# contamintation estimate works as a hard cutoff
# 如果过滤所依据的注释缺乏的话,这个过滤步骤会被跳过
gatk FilterMutectCalls \
 -V somatic_m2.vcf.gz \
 --contamination-table tumor_calculatecontamination.table \
 -O 9_somatic_oncefiltered.vcf.gz

5.5 (可选)检测并过滤artifact

FilterByOrientationBias在sequence context artifacts(如OxoG,FFPE)的基础上进行过滤.这步是可选的,如果执行了这一步,还需要额外执行一次FilterMutectCalls.本步所需的metrics可以通过Picard的CollectSequencingArtifactMetrics获取.当然,gatk里面整合了这个方法.

OxoG在自发的种系颠换里面; 在小鼠的生殖细胞的DNA复制时,the oxidized base 8-oxoguanine (8-oxoG)会导致自然且可遗传的G->T颠换(transversion). FFPE是指福尔马林固定的处理方式,它会导致C->T转换(transition).所有在建库时,会可能出现这两种形式的假阳性变异.

# 1\. 使用CollectSequencingArtifactMetrics从sequence context artifacts里面获取metrics
# 这里会将artifacts进行分类,比如分为preadapter(artifacts occur before hybrid selection)和baitbias(artifacts occur during hybrid selection)
gatk CollectSequencingArtifactMetricx \
 -I tumor.bam \
 -O 10_tumor_artifact \
 --FILE_EXTENSION '.txt' \
 -R path/hg38/Homo_sapiens_assembly38.fasta

# 你也可以在picard里面直接使用
java -jar picard.jar \
 CollectSequencingArtifactMetricx \
 ......

# output => five metrics files, pre_adapter_detail_metrics in five is for FilterByOritentationBias

# 2\. 使用FilterByOrientationBias进行定位偏差过滤
# orientation指F1R2, F2R1这类的定位, reads比对到参考序列的定位.
gatk FilterByOrientationBias \
 -A G/T -A C/T \
 -V 9_somatic_oncefiltered.vcf.gz \
 -P tumor_artifact.pre_adapter_detail_metrics.txt \
 -O 11_somatic_twicefiltered.vcf.gz

|

六、Somatic variant discovery: CNV

多数癌症是由CNV引起,如HER2,EGFR,ERBB2等是因为CNV增多引起, APC,BRCA1,BRCA2,PTEN是因为CNV减少引起,这就使得CNV的检测至关重要。

6.1 如何确定CNV的产生CNV的表征方法?

copy number variants会导致覆盖度不平衡(coverage imbalance).但是coverage会因WES中靶标位置和所用kits的不同而不同. 我们可以通过copy number profile来对cnv绝对定量，获得每个基因座的拷贝数(copy number of locus),也可以通过copy ratio profile来对CNV进行相对定量，获得每个基因座的拷贝数在平均倍数性(average ploidy)的占比。虽然copy ratio可以把coverage imbalance的问题解决，但外显子的原始copy ratio数据噪声很大。需要除去噪声。

6.2 获得目标区域的原始coverage counts

在收集coverage counts之前,我们先用一个genomic intervals list确定分析的分辨率(也就是我们的目标区域). 这对wes或者wgs数据分析来说都是必要的.对于外显子数据,我们使用捕获试剂盒进行补充;对于全基因组数据,我们把参考基因组分成等长的intervals或者bins, 也可以不指定这个interval_list.不管怎么样, 都是使用PreprocessIntervals.

6.2.1 根据目标区域interaval对参考序列设定interval

# exome, wes
gatk PreprocessIntervals \
    -R reference.fasta \
    -L intervals.interval_list \
    --bin-length 0 \
    --padding 250 \
    -O preprocessed_intervals.interval_list
# genome, wgs
gatk PreprocessIntervals \
    -R reference.fasta \
    --bin-length 1000 \
    --padding 0 \
    -O preprocessed_intervals.interval_list
# parameters
#-L  可选的参数.指定目标intervals_list
#--bin-length    为不同的数据类型指定适当的值,比如在这里wes是0,wgs是1000
#--interval-merging-rule 设定为OVERLAPPING_ONLY将阻止工具合并相邻的intervals
#--padding   指定intervals的padding(间隔), 默认值250在大部分wes中都效果良好

6.2.2 获得对应interval的coverage counts

Fragment counts in each genomic target/bin allow estimates of segmented copy ratio. CollectFragmentCounts计算的是paired end fragments的覆盖度, 在interval中心重叠的fragment只会被计一次. CollectReadCounts在GATK4.0.3.0中取代了CollectFragmentCounts,而且它还会将在interval重叠的reads计数. Targets/bins + BAM => collectFragmentCounts => integer read counts.

gatk CollectFragmentCounts \
    -I sample.bam \
    -L intervals.interval_list \
    --interval-merging-rule OVERLAPPING_ONLY \
    -O sample.counts.hdf5

#--format    结果格式默认为hdf5.可以指定为TSV,生成tsv格式的结果. hdf5文件可以使用HDFView软件查看.
#filter      这个方法在底层调用了多个过滤方法,如NotDuplicateReadFilter,FirstOfPairReadFilter, etc.

6.3 生成CNV panel of normals(CNV, PoN)

在创建PoN时,CreateReadCountPanelOfNormals使用了Singular Value Decomposition来对counts和interavals进行抽象. PoN用于PCA去噪. PoN里面的正常样品应该与测序方法相对应,这对wes十分重要,因为wes的捕获步骤会引入靶标特异地噪声. Fragment counts + Annotated Intervals => CreateReadCountPanelOfNormals => Panel of Normals

gatk --java-options '-Xmx6500m' CreateReadCountPanelOfNormals \
    -I HG0013.alt_bwamem_GRCh38Dh.20150826.GBR.exome.counts.hdf5 \
    -I HG0013.alt_bwamem_GRCh38Dh.20150826.GBR.exome.counts.hdf5 \
    ... \
    --minimum-interval-median-percentile 5.0 \
    -O cnv.pon.hdf5
# or
gatk CreateReadCountPanelOfNormals \
    -I sample_1.counts.hdf5 \
    -I sample_2.counts.tsv \
    ... \
    --annotated-intervals annotated_intervals.tsv \
    -O cnv.pon.hdf5

# select parameters
#-I  为每个样品提供integer read coverage counts. 格式可以是TSV和HDF5.可以接受单个样品,比如和tumor配对的正常样品(matched-normal)
#--number-of-eigensamples    默认值为20,这个会自动校正.一般来说,在去噪时,这个值越高,敏感性就越低.你可以进行敏感性分析以选取一个合适的.
#--spark-master  如果要使用spark的话,该参数指定.
# filter parameters
#--minimum-interval-median-percentile    默认值为10.0,会将覆盖度的中位数低于该值的targets/bins/intervals给排除.
#--maximum-zeros-in-sample-percentage    默认值为5.0,会将任何具有高于5%的0覆盖率target/bins的样品排除
#--maximum-zeros-in-interval-percentage  默认值为5.0,会将任何具有高于5%的0覆盖率的interval的样品排除
#--extreme-sample-median-percentile      默认值为2.5,校正后的中位数不在2.5-97.5%之间的样品将被排除
#--do-impute-zeros                       默认为true.会把0覆盖率区域的覆盖率值指定为非零值的中位数.
#--extreme-outlier-truncation-percentile 默认值为0.1.把不在0.1-99.9%区域的值分别设为最近的值(<0.1设为0.1)

6.4 Denoise coverage data

gatk --java-options '-Xmx12g' DenoiseReadCounts \
    -I sample.counts.hdf5 \
    --count-panel-of-normals cnv.pon.hdf5 \
    --standardized-copy-ratios sample.standardizedCR.tsv \
    --denoised-copy-ratios sample.denoisedCR.tsv

#--number-of-eigensamples    默认为null,使用PoN中样品的最大可用数.在这里改变该值会改变结果的分辨率.
#--annotated-intervals   可以使用该参数据指定进行GC-bias correlation. 参数的输入文件可以通过AnnotateIntervals产生.
# output => sample.standardizedCR.tsv(标准copy ratio), sample.denoisedCR.tsv(denoised copy ratios)

6.5 可视化去噪的效果

使用PlotDenoisedCopyRatios对得到的标准化copy ratio和去噪的copy ratio作图, 显示去噪的效果。

gatk PlotDenoisedCopyRatios \
    --standardized-copy-ratios sample.standardizedCR.tsv \
    --denoised-copy-ratios sample.denoised.tsv \
    --sequence-dictionary Homo_sapiens_assembly38.dict \
    --minimum-contig-length 46709983 \
    --output plots \
    --output-prefix sample
# parameters
#--sequence-dictionary 指定用于比对的参考序列的字典
#--output-prefix 指定输出文件的前缀
#--output 输出文件路径.会产生6个文件.
#[sample.denoised.png], plots the standardized and denoised read counts. [sample.denoisedLimit4.png], y轴被限定在0-4.
#[sample.standardizedMAD.txt, sample.denoisedMAD.txt], MAD for standardized/denoised copy ratios, 
#[sample.deltaMAD.txt], difference between standardized MAD and denoised MAD.
#[sample.scaledDeltaMAD.txt], fractional difference / standardized MAD.(s-MAD - d-MAD)/s-MAD

6.6 计算常见生殖系变异位点的ref和alt alleles

等位基因特异覆盖率的收集仅仅是原始覆盖率的收集,并没有经过任何统计参考.对于实际样本和匹配的对照比对来说,CollectAllelicCounts都是分别收集allele counts.但是对于匹配的对照进行分析的话,样品和对照一定要匹配,不然后续的ModelSegments会发生错误.

这里没有使用GetPileupSummaries，这是因为他会使用CalculateContamination估测跨样本污染的counts.它通过群体等位基因频率来设定进行counts收集的条件.而CollectAllelicCounts用到的过滤器较少,但二者都会使用MappingQualityReadFilter进行过滤.前者的过滤阈值默认是20,可设置到30.而后者使用的阈值则是50.此外,前者还使用base-quality.

# 为匹配的对照(matched-control)收集位于种系变异位点的counts
gatk --java-options '-Xmx3g' CollectAllelicCounts \
    -L chr17_theta_snps.interval_list -I normal.bam \
    -R path/ref/Homo_sapiens_assembly38.fasta -O normal.allelicCounts.tsv
# 为每个样品在相同的位点收集counts
gatk --java-options '-Xmx3g' CollectAllelicCounts \
    -L chr17_theta_snps.interval_list -I tumor.bam \
    -R path/ref/Homo_sapiens_assembly38.fasta -O tumor.allelicCounts.tsv
# parameters
#-L  必需参数.指定一个或多个genomic intervals,可以是intervals list或VCF格式. 里面包含的sites应该代表了常见的或这样品特异地种系变异位点(而且只是SNP位点).Indel和mixed-variant-type将被忽略.
#--minimum-base-quality  该方法会使用MappingQualityReadFilter进行过滤,默认的过滤值是20, 即MAPQ>=20.

6.7 把contiguousCopyRatios分组成segments
对于allelic copy ratios, 在control in a paired analysis或者 case in case-only analysis, ModelSegments只使用杂合体的位点.这里把allelic copy ratios 定义为alternate-allele fraction(A-AF).

gatk --java-options '-Xmx4g' ModelSegements \
    --denoised-copy-ratios tumor.denoisedCR.tsv \
    --allelic-counts tumor.allelicCounts.tsv \
    --normal-allelic-counts normal.allelicCounts.tsv \
    --ouput-prefix tumor \
    --output output_dir

#output: tumor.modelBegin.seg, tumor.modelFinal.seg, tumor.cr.seg, (data on the segments)
#tumor.modelBegin.af.param, tumor.modelBegin.cr.param, (global parameters for copy ratios,cr, allele fractions,af)
#tumor.modelFinal.af.param, tumor.modelFinal.cr.param,
#tumor.hets.normal.tsv, tumor.hets.tsv (hets:the allelic counts for the control's heterogygous sites, normal: counts for the matched control)

#--denoised-copy-ratios  输入的数据,可以是copy-ratios,也可以是--allelic-counts.下同.
#--allelic-counts        同上.
#--minimum-total-allele-count    默认值为30.对于allelic copy ratios,只考虑read depth coverage超过设定值的位点.
#--genotyping-homozygous-log-ratio-threshold 默认值为-10.0.指定建模时假定的杂合体的位点数目.
#--maximum-number-of-smoothing-iterations    指定建模时的平滑迭代次数,默认为25.

6.7.1 ModelSegments运行的三个阶段

1. 对杂合位点进行分型,依据覆盖深度(depth)和copy-ratio intervals重叠位点进行过滤.对照中杂合的位点写入产生hets.normal.tsv文件,样本中相同的位点写入hets.tsv文件.
2. 执行多维Kernel segmentation.这里使用到allelic counts和denoised copy ratios.
3. 进行Markov-Chain Monte Carlo抽样和平滑分割.初始拟合模型写入modelBegin.seg,及对应的参数设置.最终拟合模型写入modelFinal.seg,及对于的参数设置.

6.8 检测正常拷贝数、拷贝数增多、减少的segments

call copy number events. 此步骤对于画出segmentation结果图并不是必需的. 使用ModelSegments的输出cr.seg作为输入

gatk CallCopyRatioSegments \
    -I tumor.cr.seg -O tumor.called.seg

#--neutral-segment-copy-ratio-lower-bound    默认值为0.9.设置copy-neutral segments的范围.下同.
#--neutral-segment-copy-ratio-upper-bound    默认值为1.1.同上.

6.9 可选.对分割(segmentation result)结果画图

visualizes copy and allelic ratio segmentation results.

gatk PlotModeledSegments \
    --denoised-copy-ratios tumor.denoisedCR.tsv --allelic-counts tumor.hets.tsv \
    --segments tumor.modelFinal.seg --minimum-contig-length 46709983 \
    --sequence-dictionary Homo_sapiens_assembly38.dict \        # also contigs_to_plot.dict
    --output-prefix tumor -O output_dir

#--sequence-dictionary   比对时使用的参考序列的字典
#--minimum-contig-length 对于GRCh38版本的参考序列来说,默认值是1,000,000-46,709,983.设置后,低于设置值的alternate/decoy contigs将会被忽略,不会画到图里面去.
#(-L)    
#你无法通过指定-L来可视化特定染色体或interval区域.你可以选择其他方式：其一, 把sequence-dictionary设置成你感兴趣的contigs, contigs_to_plot.dict。其二, 把感兴趣的CNV数据转化为BED格式,在IGV里面可视化.其三, 把数据导入R,选取感兴趣数据部分进行可视化.