这步，才正式将比对结果，call出SNP，给出一个非常通用的，大家都熟悉的VCF文件 (variant calling files)

SNP， indel是分开call的，同时会对质量进行过滤。从raw data生成PASS的data

#Generate snps raw vcf file

#Using GATK UnifiedGenotyper to generate snps.raw.vcf

 mkdir $PWDS/variants

 java -Xmx${heap}m  \

    -jar $gatk \

    -R $REF \

    -T UnifiedGenotyper \

    --dbsnp $DBSNP \

    -I $PWDS/${subjectID}.realigned.recal.bam \

    --out $PWDS/variants/${subjectID}.snps.raw.vcf \

    -stand_call_conf 30.0 \

    -stand_emit_conf 10.0 \

    -out_mode EMIT_VARIANTS_ONLY \

    -l INFO \

    -A DepthOfCoverage \

    -A HaplotypeScore \

    -A InbreedingCoeff \

    -glm SNP  \

    -nt 1  \

    -L $ExonFile

 rm -f $PWDS/variants/*.idx

#Generate indels raw vcf file

Using GATK UnifiedGenotyper to generate indels.raw.vcf

 java -Xmx${heap}m  \

    -jar $gatk \

    -R $REF \

    -T UnifiedGenotyper \

    --dbsnp $DBSNP \

    -I $PWDS/${subjectID}.realigned.recal.bam \

    --out $PWDS/variants/${subjectID}.indels.raw.vcf \

    -stand_call_conf 30.0 \

    -stand_emit_conf 10.0 \

    -out_mode EMIT_VARIANTS_ONLY \

    -l INFO \

    -A DepthOfCoverage \

    -A HaplotypeScore \

    -A InbreedingCoeff \

    -glm INDEL  \

    -nt 1  \

    -L $ExonFile

 rm -f $PWDS/variants/*.idx

#Variant filtration generate indels and indels.PASS vcf files

Using GATK VariantFiltration to generate snps.vcf and snps.PASS.vcf files

 java -Xmx${heap}m  \

  -jar $gatk \

  -l INFO -T VariantFiltration \

  -R $REF \

  -o $PWDS/variants/${subjectID}.indels.vcf \

  -V:VCF $PWDS/variants/${subjectID}.indels.raw.vcf \

  --filterExpression "QD < 2.0" \

  --filterName "QDFilter" \

  --filterExpression "ReadPosRankSum < -20.0" \

  --filterName "ReadPosFilter" \

  --filterExpression "FS > 200.0" \

  --filterName "FSFilter" \

  --filterExpression "MQ0 >= 4 && ((MQ0 / (1.0 * DP)) > 0.1)" \

  --filterName "HARD_TO_VALIDATE" \

  --filterExpression "QUAL < 30.0 || DP < 6 || DP > 5000 || HRun > 5" \

  --filterName "QualFilter"

 ##这些参数都是可以改的，比如quanlity默认是大于30以上保留，测序深度是6X，等待，根据自己实际要求修改。

 mline2=`grep -n "#CHROM" $PWDS/variants/${subjectID}.indels.vcf | cut -d':' -f 1`

 head -n $mline2 $PWDS/variants/${subjectID}.indels.vcf > $PWDS/variants/headindels.vcf

 cat $PWDS/variants/${subjectID}.indels.vcf   \

 | grep PASS | cat $PWDS/variants/headindels.vcf - >   \

 $PWDS/variants/${subjectID}.indels.PASS.vcf

 rm -f $PWDS/variants/*.idx

 rm -f $PWDS/variants/headindels.vcf

#Variant filtration generate snps and snps.PASS vcf files

Using GATK VariantFiltration to generate indels.vcf and indels.PASS.vcf

 java -Xmx${heap}m  \

  -jar $gatk \

  -l OFF -T VariantFiltration \

  -R $REF \

  -o $PWDS/variants/${subjectID}.snps.vcf \

  --variant $PWDS/variants/${subjectID}.snps.raw.vcf \

  --mask $PWDS/variants/${subjectID}.indels.PASS.vcf \

  --maskName InDel \

  --clusterSize 3 \

  --clusterWindowSize 10 \

  --filterExpression "QD < 2.0" \

  --filterName "QDFilter" \

  --filterExpression "MQ < 40.0" \

  --filterName "MQFilter" \

  --filterExpression "FS > 60.0" \

  --filterName "FSFilter" \

  --filterExpression "HaplotypeScore > 13.0" \

  --filterName "HaplotypeScoreFilter" \

  --filterExpression "MQRankSum < -12.5" \

  --filterName "MQRankSumFilter" \

  --filterExpression "ReadPosRankSum < -8.0" \

  --filterName "ReadPosRankSumFilter" \

  --filterExpression "QUAL < 30.0 || DP < 6 || DP > 5000 || HRun > 5" \

  --filterName "StandardFilters" \

  --filterExpression "MQ0 >= 4 && ((MQ0 / (1.0 * DP)) > 0.1)" \

  --filterName "HARD_TO_VALIDATE"

 rm -f $PWDS/variants/*.idx

  mline=`grep -n "#CHROM" $PWDS/variants/${subjectID}.snps.vcf | cut -d':' -f 1`

 head -n $mline $PWDS/variants/${subjectID}.snps.vcf > $PWDS/variants/head.vcf

 cat $PWDS/variants/${subjectID}.snps.vcf | grep PASS   \

  | cat $PWDS/variants/head.vcf - >    \

  $PWDS/variants/${subjectID}.snps.PASS.vcf

 rm -f $PWDS/variants/head.vcf

#evaluating snps

Using GATK snps eval to generate snps.PASS.vcf.eval file

 java -Xmx${heap}m  \

   -jar $gatk \

   -T VariantEval \

   --dbsnp $DBSNP \

   -R $REF \

   -eval $PWDS/variants/${subjectID}.snps.PASS.vcf \

   -o $PWDS/variants/${subjectID}.snps.PASS.vcf.eval

 rm -f $PWDS/variants/*.idx