Unified Genotyper: --allSitePLs option

prepagamprepagam Member
edited December 2014 in Ask the GATK team

Hi

I've been trying to use the --allSitePLs option with Unified Genotyper with GaTK v3.2-2-gec30cee.

My command is:

java -Xmx${MEM} -jar ${gatk_dir}/GenomeAnalysisTK.jar -R ${ref_dir}/${genome} -T UnifiedGenotyper \
    --min_base_quality_score 30 \
    -I ${in_dir}/"34E_"${REGION}"_ddkbt_RS74408_recalibrated_3.bam" \
    -I ${in_dir}/"34E_"${REGION}"_ddber_RS86405_recalibrated_3.bam" \
    --intervals $sites_dir/$file \
    -o ${out_dir}/"38F_"${REGION}"_UG_allPL.vcf" \
    --output_mode EMIT_ALL_SITES \
   --allSitePLs 

I ran the same command with and without the --allSitePLs option and compared, and the output seems strange.

Specifically, with --allSitePLs - the FILTER column: 5053700 sites = lowqual, 19303 = . and 5059568 had QUAL < 30. By contrast WITHOUT the --allSitePLs 5067411 = lowqual, 5592 = . and 11460 had QUAL < 30. I don't understand why does adding this option changes the QUAL so much when I'm running UG on the exact same data but just requesting that I get a PL for all sites.

Lastly, how is the specific ALT allele selected? Is it random - because they all seem equally unlikely in the example below, where only one allele was seen in both my samples.

## 2 lines from the output
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  ddber_RS86405   ddkbt_RS74408
chr01   1753201 .   C   A   0   LowQual AC=0;AF=0.00;AN=4;DP=43;Dels=0.00;FS=0.000;HaplotypeScore=0.0000;MLEAC=0;MLEAF=0.00;MQ=60.00;MQ0=0 GT:AD:APL:DP:GQ:PL   0/0:17,0:655,51,0,655,51,655,655,51,655,655:17:51:0,51,655  0/0:26,0:1010,78,0,1010,78,1010,1010,78,1010,1010:26:78:0,78,1010
chr01   1753202 .   T   A   0   LowQual AC=0;AF=0.00;AN=4;DP=43;Dels=0.00;FS=0.000;HaplotypeScore=0.0000;MLEAC=0;MLEAF=0.00;MQ=60.00;MQ0=0 GT:AD:APL:DP:GQ:PL   0/0:17,0:630,630,630,630,630,630,48,48,48,0:17:48:0,48,630  0/0:26,0:1027,1027,1027,1027,1027,1027,78,78,78,0:26:78:0,78,1027
Post edited by Geraldine_VdAuwera on

Best Answer

Answers

  • prepagamprepagam Member

    Hi Geraldine
    Thanks - and no worries. I have tried Hap Caller for two recent projects, and it just didn't work out well - when I viewed the calls in IGV it was very odd, whereas UG gave intuitive and expected results. I work on non-model species, and think it is likely related to this, as I know Hap Caller performs well in humans, mice .....

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Ah, I see. It's possible that the results seemed odd when viewed in IGV because HC remaps reads, so sometimes it looks like the calls don't match the data. If you use the -bamout option to view the data as remapped by HC it makes a lot more sense.

Sign In or Register to comment.