bwa+samtools+picardtools+GATK ca

作者: bio | 来源:发表于2015-09-07 18:07 被阅读16474次

    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

    相关文章

      网友评论

      • ad601e694e74:call Indel之前少了局部重新比对,第一步,通过运行RealignerTargetCreator来确定要进行重新比对的区域。第二步,通过运行IndelRealigner在这些区域内进行重新比对。也没有对bam文件里reads的碱基质量值进行重新校正,
      • 4ed0d00187da:博主您好!我想要请教您一个问题,就是--clusterwindowsize这个参数,为什么10bp内SnpCluster是被认为不可靠的呢,标注的参考文献我没能找到,如果可以的话,能告诉我文章的名字吗?谢谢!

      本文标题:bwa+samtools+picardtools+GATK ca

      本文链接:https://www.haomeiwen.com/subject/upcrcttx.html