GATK Best Practices — step1 数据预处
一、数据预处理（Data Pre-processing）介绍
这一步主要是对BWA得到的bam进行一系列处理（包括sort/dup的标记/碱基质量值矫正），最终得到GATK call变异最终输入的bam。
官方给出了三个Pipeline：
1，Prod* germline short variant per-sample calling ：the Broad's Genomic Sequencing Platform平台自测序数据
2，$5 Genome Analysis Pipeline ：全套分析pipeline
3，Generic data pre-processing ：数据前处理
这一步官网也给出了流程的主要步骤（这一步是每个样品单独跑的）：
1，Map to Reference 比对到参考基因组并sort
2，Mark Duplicates 去除数据产生过程中的dup
3，Base (Quality Score) Recalibration 碱基质量重矫正
这次测试用到的是第三个：Generic data pre-processing Github
git clone https://github.com/gatk-workflows/gatk4-data-processing.git
得到：
README.md
LICENSE
processing-for-variant-discovery-gatk4.hg38.wgs.inputs.json Homo hg38参数配置文件
processing-for-variant-discovery-gatk4.b37.wgs.inputs.json human b37参数配置文件
generic.google-papi.options.json
processing-for-variant-discovery-gatk4.wdl WDL文件，包含了全部实际运行的命令
二、WDL中各个task的介绍
学习processing-for-variant-discovery-gatk4.wdl，把这个文件中的命令拆分出来。
主要是看task中的command。里面包含了10个task，一一来看一下功能及命令：
1.GetBwaVersion
功能：Get version of BWA
bwa 2>&1 | \
grep -e '^Version' | \
sed 's/Version: //'
2.SamToFastqAndBwaMem
功能：Read unmapped BAM, convert on-the-fly to FASTQ and stream to BWA MEM for alignment
由于本次测试的输入不会涉及到ubam，这一步省略
3.MergeBamAlignment
功能：Merge original input uBAM file with BWA-aligned BAM file
将uBAM和BWA-aligned合并，由于本次测试的输入不会涉及到ubam，这一步省略
4.SortAndFixTags
功能：Sort BAM file by coordinate order and fix tag values for NM and UQ
gatk \
SortSam \
--INPUT i.bam \
--OUTPUT /dev/stdout \
--SORT_ORDER "coordinate" \
--CREATE_INDEX false \
--CREATE_MD5_FILE false \
| \
gatk \
SetNmMdAndUqTags \
--INPUT /dev/stdin \
--OUTPUT o.bam \
--CREATE_INDEX true \
--CREATE_MD5_FILE true \
--REFERENCE_SEQUENCE ref_fasta
5.MarkDuplicates
功能：Mark duplicate reads to avoid counting non-independent observations
gatk \
MarkDuplicates \
--INPUT i.bam \
--OUTPUT o.bam \
--METRICS_FILE i.duplicate_metrics \
--VALIDATION_STRINGENCY SILENT \
--OPTICAL_DUPLICATE_PIXEL_DISTANCE 2500 \
--ASSUME_SORT_ORDER "queryname" \
--CREATE_MD5_FILE true
6.CreateSequenceGroupingTSV
功能：Generate sets of intervals for scatter-gathering over chromosomes
7.BaseRecalibrator
功能：Generate Base Quality Score Recalibration (BQSR) model
gatk \
BaseRecalibrator \
-R ref_fasta \
-I input.bam \
--use-original-qualities \
-O o.recal_data.csv \
--known-sites dbSNP_vcf \
--known-sites \$(Mills_and_1000G_gold) \
--known-sites $(phase_1000G)
-L Bed
-L参数指定了 bed文件，默认是无。指定这个参数，只会输出指定区域的突变信息
8.GatherBqsrReports
功能：Combine multiple recalibration tables from scattered BaseRecalibrator runs
一个样品分开上机或者一个样品多条lane的情况，本次测试用不到这一步，忽略
9.ApplyBQSR
功能：Apply Base Quality Score Recalibration (BQSR) model
gatk \
ApplyBQSR \
-R ref_fasta \
-I ${input_bam} \
-O ${output_bam_basename}.bam \
-L bed \
-bqsr o.recal_data.csv \
--static-quantized-quals 10 --static-quantized-quals 20 --static-quantized-quals 30 \
--add-output-sam-program-record \
--create-output-bam-md5 \
--use-original-qualities
-L参数指定了 bed文件，默认是无
10.GatherBamFiles
功能：Combine multiple recalibrated BAM files from scattered ApplyRecalibration runs
一个样品分开上机或者一个样品多条lane的情况，一般用不到这一步，忽略
三、执行命令
以上，Generic data pre-processing这一步就学习完了，接下来需要找一些数据进行实际的操作（4/5/7/9步），需要用到：
成对的bam（N/T分别经过BWA比对得到）；
Ref ucsc.hg19.fasta（step0中有下载链接）；
三个数据库:dbsnp \ Mills_and_1000G_gold \ phase_1000G（step0中有下载链接）；
运行脚本(以normal样品为例)：
echo MarkDuplicates is Begin
time gatk MarkDuplicates --INPUT samplename_N.bam --OUTPUT samplename_N.rmdup.bam --METRICS_FILE samplename_N.duplicate_metrics --VALIDATION_STRINGENCY SILENT --OPTICAL_DUPLICATE_PIXEL_DISTANCE 2500 --ASSUME_SORT_ORDER "queryname" --CREATE_MD5_FILE true --CREATE_INDEX true
echo MarkDuplicates is End
echo SortSam is Begin
time gatk SortSam --INPUT samplename_N.rmdup.bam --OUTPUT /dev/stdout --SORT_ORDER "coordinate" --CREATE_INDEX false --CREATE_MD5_FILE false | gatk SetNmMdAndUqTags --INPUT /dev/stdin --OUTPUT samplename_N.sorted.bam --CREATE_INDEX true --CREATE_MD5_FILE true --REFERENCE_SEQUENCE ucsc.hg19.fasta
echo SortSam is End
echo BaseRecalibrator is Begin
time gatk BaseRecalibrator -R ucsc.hg19.fasta -I samplename_N.sorted.bam --use-original-qualities -O samplename_N.recal_data.csv --known-sites dbsnp_138.hg19.vcf --known-sites Mills_and_1000G_gold_standard.indels.hg19.vcf --known-sites 1000G_phase1.indels.hg19.vcf -L Covered.bed
echo BaseRecalibrator is End
echo ApplyBQSR is Begin
time gatk ApplyBQSR -R ucsc.hg19.fasta -I samplename_N.sorted.bam -O samplename_N.recal.bam -bqsr samplename_N.recal_data.csv --static-quantized-quals 10 --static-quantized-quals 20 --static-quantized-quals 30 --add-output-sam-program-record --create-output-bam-md5 --use-original-qualities --create-output-bam-index true -L Covered.bed
echo ApplyBQSR is End
得到：
samplename_N.rmdup.bam
samplename_N.rmdup.bam.md5
samplename_N.duplicate_metrics
samplename_N.sorted.bam
samplename_N.sorted.bam.md5
samplename_N.sorted.bai
samplename_N.recal_data.csv
samplename_N.recal.bam（下面步骤中要用到的）
samplename_N.recal.bam.md5
samplename_N.recal.bai（下面步骤中要用到的）
总共四步，noraml样品每一步的运行时间分别如下：
5m45.381s
6m7.366s
4m53.107s
6m3.656s
对应的tumour运行时间：
7m42.873s
7m54.682s
3m49.228s
8m20.963s
四、注意事项
需要注意的是，这一步中，WDL虽然给出的pipeline中，SORT在前，RMDUP在后，但是，流程中需要将去dup这一步放在sor前面，否则，在运行ApplyBQSR这一步时会报错，显示输入非sort：
WARN ReadUtils - Skipping index file creation for: samplename_N.recal.bam. Index file creation requires reads in coordinate sorted order.
导致最后的samplename_N.recal.bam结果没有index（不会报错，但是会警告）。
对于这一步影响不大，但是Index对于后面步骤call SNP/SNV/INDEL是必须的，如果没有index会报错，提示需要用samtools对bam进行index操作，为了避免增加工作量，需要注意流程中将SORT在RMDUP后面一步。
这一步的学习基本完成。
下一个step是对胚系突变SNP/INDEL（Germline SNPs + Indels）的学习
