call SNP GATK bwa+samtools 参考文档
GATK call snp 流程图
一、原始数据处理
测序得到的RAW DATA，使用NGSQCToolkit进行过滤
$PATH/IlluQC.pl -pe read1.fastq read2.fastq 2 A -l 70 -s 20 -p 2 -z g
- -pe paire-end数据
- 2 = Paired End DNA Library
- A = Automatic detection of FASTQ variant
- -l The cut-off value for percentage of read length that should be of given quality
- -s参数 -cutOffQualScore read质量值cutoff
- -z g=输出结果为压缩文件
二、preGATK处理
GATK call snp 需要排序后的bam文件，bam文件由sam文件转化而来。sam文件bwa、bowtie等软件生成。
1.bwa 生成sam文件
#建立参考基因组索引
$path/bwa index ref.fa
#*find input reads SA coordnate *
bwa aln chlorandrum.fa read1.filtered.fq>read1.fq.sai
bwa aln chlorandrum.fa read2.filtered.fq >read2.fq.sai
#sai 格式 to sam 格式
bwa sampe chlorandrum.fa -r "@RG\tID:<ID>\tLB:<LIBRARY_NAME>\tSM:<SAMPLE_NAME>\tPL:ILLUMINA" read1.fq.sai read2.fq.sai read1.filtered.fq read2.filtered.fq >test_read.sam
如果GATK call SNP 必须用-r 参数，
@RG\tID:<ID>\tLB:<LIBRARY_NAME>\tSM:<SAMPLE_NAME>\tPL:ILLUMINA 如果多样品 call SNP ,<SAMPLE_NAME> 更改成相应样品名称，否则将会被当成一个样品处理
2.对sam文件进行重排序
samtools faidx chlorandrum.fa
java -jar $picard CreateSequenceDictionary R=ref.fa O=chlorandrum.dict
java -jar $picard ReorderSam I=test_reads.sam O=test_reads.reordered.sam R=ref.fa
3.sam文件转化成bam文件
samtools view -bs test_reads.reordered.sam -o test.bam
4.对bam文件进行排序
java -jar $path/picard.jar SortSam INPUT=3.test_reads.reordered.sam2.bam OUTPUT=4.test_reads.reordered.sam2.sorted.bam SORT_ORDER=coordinate
5.Duplicates Marking
测序原理是随机打断，那么理论上出现两条完全相同的read的概率是非常低的，而且建库时PCR扩增存在偏向性，因此标出完全相同的read。
java -jar $picard MarkDuplicates REMOVE_DUPLICATES= false MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000 INPUT=test.bam OUTPUT=test.repeatmark.bam METRICS_FILE=test.bam.metrics
6.bam文件生成索引
samtools index test.repeatmark.bam
三、GATK变异检测
1.生成raw vcf 文件
java -Xmx2G -jar $GATK \
-T HaplotypeCaller \
-R ref.fa \
-I test.repeatmark.bam \ #若多样品，则-I sample1.bam -I sample2.bam
--genotyping_mode DISCOVERY \
-stand_emit_conf 10 \
-stand_call_conf 30 \
-o raw_variants.vcf
2.select SNP
java -Xmx2g -jar $GATK -R $REF -T SelectVariants -o $Slect_SNP --variant $RAW_vcf -selectType SNP 2>select_snp.err
3.select indel
java -Xmx2g -jar $GATK -R $REF -T SelectVariants -o $Slest_INdel --variant $RAW_vcf -selectType INDEL 2>select_indel.err
4.filter SNP
java -Xmx4g -jar $GATK -R $REF -T VariantFiltration --variant $Slect_SNP --clusterSize 4 --clusterWindowSize 10 --maskName aroundIndel --mask $Slest_INdel -maskExtend 3 --filterName "lowMQRankSum" --filterExpression "MQRankSum < -12.5" --filterName "highFS" --filterExpression "FS > 60.0" --filterName "lowReadPosRankSum" --filterExpression "ReadPosRankSum < -8.0" --filterName "lowMQ" --filterExpression "MQ < 40.0" --filterName "lowQD" --filterExpression "QD < 2.0" --out $Filter_SNP --genotypeFilterName "lowDP" --genotypeFilterExpression "DP < 8.0" >filter_snp.err
- --clusterSize --clusterwindowsize 这两个参数一起用可以设定指定长度内的variants数量的最大值，例如-window设成10，-cluster设成3，表示10bp以内的variants的数量不应当超过3个，如果超过了则会用SnpCluster标示出来。一般这种聚集的SNP可以被认为是不可靠的(10bp里面超过3个(Bowen et al., 2011))。
- --mask参数用来筛选与某个指定文件中记录一致的位点，例如筛选repeat区域的位点的时候就可以使用这个参数，这边给的是一个文件，支持的格式很多，所以适用的范围其实很广，这边暂时先简单的介绍一下；--maskName同样是来设置被FILTER掉的位点的标识的。
- --filter --filterName 这两个参数也是一起使用的，可以一个第一个参数设置筛选用的条件表达式，第二个参数设置这个筛选在结果里面标识的名称，两个参数设置的时候必须是成对的，然后可以设置多次
-
--missingValuesInExpressionsShouldEvaluateAsFailing (Boolean with default value false)
最后在讲一下这个参数，当筛选标准比较多的时候，可能有一些位点没有筛选条件当中的一条或几条，例如下面的这个表达式
QUAL < 30.0 || QD < 2.0 || FS > 60.0 || MQ < 40.0 || HaplotypeScore > 13.0
并不一定所有位点都有这些信息，这种情况下GATK运行的时候会报很多WARNING信息，用这个参数可以把这些缺少某些FLAG的位点也给标记成没有通过筛选的。 - FS：使用Fisher’s精确检验来检测strand bias而得到的Fhred格式的p值。该值越小越好。一般进行filter的时候，可以设置 FS < 10～20。
- MQ：RMS Mapping Quality
- QD：Variant Confidence/Quality by Depth
- ReadPosRankSum：Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias
- MQRankSum：Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities
