bwa+samtools+picardtools+GATK ca

2015-09-07 本文已影响16474人 bio

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