EMIT_ALL_CONFIDENT_SITES option, omitting non-'LowQual' sites

Jakemorris89Jakemorris89 UKMember
edited August 2015 in Ask the GATK team

Hi,

I've run UnifiedGenotyper on a BAM file with both EMIT_ALL_SITES and EMIT_ALL_CONFIDENT_SITES. I've noticed some of the calls that have been omitted with the EMIT_ALL_SITES option seem to be of comparable quality to others that have been emitted. For example, sites like HE670865 368605 are removed as non-confident sites while the site 368591 has not been.

I understand why the positions flagged as "LowQual" are not present when using EMIT_ALL_CONFIDENT_SITES. But I can't see why other positions (such as HE670865 368605 and HE670865 368600) are not being emitted. Of particular importance is the fact that some of these seemingly 'good sites' that have been removed are SNPs that might be of biological importance.

These are the two commands used together with some relevant lines from the resulting vcfs:

java -Xmx4g -jar GenomeAnalysisTK.jar -T UnifiedGenotyper -nt 4 -out_mode EMIT_ALL_SITES -baq CALCULATE_AS_NECESSARY
-hets 0.015 -R genome_ref.fasta -o output_all.vcf -I input_divergentscaffs.realign.bam &>
divergentscaffs_out.GATKlog &

java -Xmx4g -jar GenomeAnalysisTK.jar -T UnifiedGenotyper -nt 4 -out_mode EMIT_ALL_CONFIDENT_SITES -baq CALCULATE_AS_NECESSARY
-hets 0.015 -R genome_ref.fasta -o output_confident.vcf -I input_divergentscaffs.realign.bam &>
divergentscaffs_confident.GATKlog &

These are a sample of lines from the outputted VCF file when using EMIT_ALL_SITES:

HE670865 368589 . A . 46.89 . AN=16;DP=441;MQ=43.26;MQ0=13 GT:DP 0/0:35 0/0:45 0/0:58 0/0:38 0/0:75 0/0:64 0/0:64 0/0:47
HE670865 368590 . A . 50.15 . AN=16;DP=447;MQ=43.41;MQ0=13 GT:DP 0/0:35 0/0:45 0/0:60 0/0:39 0/0:76 0/0:66 0/0:65 0/0:47
HE670865 368591 . A . 39.12 . AN=16;DP=448;MQ=43.32;MQ0=12 GT:DP 0/0:36 0/0:45 0/0:61 0/0:39 0/0:77 0/0:66 0/0:66 0/0:46
HE670865 368592 . C T 5437.28 . AC=16;AF=1.00;AN=16;BaseQRankSum=1.036;DP=445;Dels=0.06;FS=0.000;HaplotypeScore=99.1839;MLEAC=16;MLEAF=1.00;MQ=43.48;MQ0=9;MQRankSum=4.214;QD=12.22;ReadPosRankSum=6.263 GT:AD:DP:GQ:PL 1/1:1,27:29:27:296,27,0 1/1:0,45:45:30:344,30,0 1/1:6,43:52:33:327,33,0 1/1:1,38:39:48:506,48,0 1/1:5,68:76:99:1100,102,0 1/1:4,61:65:68:1068,68,0 1/1:0,65:66:99:1035,99,0 1/1:0,46:46:69:774,69,0
HE670865 368593 . A . 50.76 . AN=16;DP=448;MQ=43.48;MQ0=9 GT:DP 0/0:29 0/0:45 0/0:53 0/0:41 0/0:76 0/0:65 0/0:66 0/0:47
HE670865 368594 . G . 69.75 . AN=16;DP=449;MQ=43.48;MQ0=9 GT:DP 0/0:28 0/0:46 0/0:51 0/0:41 0/0:76 0/0:66 0/0:66 0/0:47
HE670865 368595 . A . 66.75 . AN=16;DP=447;MQ=43.43;MQ0=9 GT:DP 0/0:30 0/0:46 0/0:50 0/0:41 0/0:75 0/0:66 0/0:66 0/0:47
HE670865 368596 . A . 72.75 . AN=16;DP=445;MQ=43.54;MQ0=9 GT:DP 0/0:31 0/0:46 0/0:48 0/0:41 0/0:76 0/0:66 0/0:67 0/0:47
HE670865 368597 . G . 66.75 . AN=16;DP=442;MQ=43.48;MQ0=9 GT:DP 0/0:25 0/0:46 0/0:41 0/0:40 0/0:74 0/0:67 0/0:67 0/0:47
HE670865 368598 . A . 75.75 . AN=16;DP=444;MQ=43.44;MQ0=9 GT:DP 0/0:25 0/0:46 0/0:41 0/0:39 0/0:75 0/0:69 0/0:67 0/0:47
HE670865 368599 . A . 72.75 . AN=16;DP=445;MQ=43.49;MQ0=9 GT:DP 0/0:24 0/0:45 0/0:39 0/0:39 0/0:76 0/0:69 0/0:68 0/0:48
HE670865 368600 . T A 4405.18 . AC=8;AF=0.500;AN=16;BaseQRankSum=2.278;DP=441;Dels=0.08;FS=6.510;HaplotypeScore=100.3222;MLEAC=8;MLEAF=0.500;MQ=43.26;MQ0=9;MQRankSum=0.498;QD=24.75;ReadPosRankSum=1.103 GT:AD:DP:GQ:PL 1/1:3,19:22:48:570,48,0 1/1:2,45:47:99:1386,114,0 1/1:4,33:37:96:1196,96,0 1/1:0,38:38:99:1307,105,0 0/0:76,2:78:99:0,178,2196 0/0:65,0:65:99:0,159,2002 0/0:70,0:70:99:0,168,2042 0/0:49,0:49:99:0,129,1549
HE670865 368601 . T . 63.75 . AN=16;DP=442;MQ=43.37;MQ0=9 GT:DP 0/0:19 0/0:47 0/0:37 0/0:38 0/0:77 0/0:66 0/0:72 0/0:48
HE670865 368602 . A . 66.75 . AN=16;DP=444;MQ=43.44;MQ0=9 GT:DP 0/0:19 0/0:47 0/0:38 0/0:38 0/0:78 0/0:65 0/0:74 0/0:48
HE670865 368603 . C . 66.75 . AN=16;DP=441;MQ=43.50;MQ0=11 GT:DP 0/0:20 0/0:47 0/0:39 0/0:37 0/0:76 0/0:64 0/0:73 0/0:48
HE670865 368604 . A . 63.75 . AN=16;DP=441;MQ=43.24;MQ0=11 GT:DP 0/0:21 0/0:47 0/0:39 0/0:38 0/0:74 0/0:64 0/0:74 0/0:48
HE670865 368605 . G . 60.75 . AN=16;DP=436;MQ=43.17;MQ0=12 GT:DP 0/0:20 0/0:47 0/0:40 0/0:38 0/0:73 0/0:63 0/0:73 0/0:47
HE670865 368606 . A . 60.75 . AN=16;DP=436;MQ=43.05;MQ0=12 GT:DP 0/0:22 0/0:45 0/0:40 0/0:39 0/0:73 0/0:63 0/0:75 0/0:46
HE670865 368607 . T C 3984.32 . AC=15;AF=0.938;AN=16;BaseQRankSum=0.225;DP=433;Dels=0.06;FS=11.987;HaplotypeScore=229.2107;MLEAC=15;MLEAF=0.938;MQ=42.94;MQ0=12;MQRankSum=1.624;QD=9.20;ReadPosRankSum=-0.463 GT:AD:DP:GQ:PL 1/1:7,18:26:12:109,12,0 1/1:2,39:44:57:475,57,0 1/1:6,36:44:51:426,51,0 1/1:1,36:39:39:328,39,0 1/1:6,67:73:44:603,44,0 1/1:0,63:63:99:1016,105,0 1/1:4,69:74:90:817,90,0 0/1:22,23:46:99:232,0,373
HE670865 368608 . A . 24.45 LowQual AN=16;DP=438;MQ=42.86;MQ0=12 GT:DP 0/0:34 0/0:45 0/0:55 0/0:39 0/0:74 0/0:63 0/0:74 0/0:47
HE670865 368609 . A . 29.97 LowQual AN=14;DP=434;MQ=42.67;MQ0=12 GT:DP ./. 0/0:44 0/0:53 0/0:39 0/0:67 0/0:61 0/0:74 0/0:47
HE670865 368610 . A . 20.62 LowQual AN=16;DP=432;MQ=42.47;MQ0=12 GT:DP 0/0:33 0/0:45 0/0:54 0/0:39 0/0:67 0/0:59 0/0:71 0/0:47
HE670865 368611 . C A 434.85 . AC=13;AF=0.813;AN=16;BaseQRankSum=0.922;DP=420;Dels=0.17;FS=0.883;HaplotypeScore=330.1362;MLEAC=14;MLEAF=0.875;MQ=42.34;MQ0=12;MQRankSum=0.871;QD=1.13;ReadPosRankSum=-2.717 GT:AD:DP:GQ:PL 0/0:17,11:28:3:0,3,23 1/1:27,10:38:3:32,3,0 1/1:23,23:46:6:64,6,0 1/1:12,16:29:3:32,3,0 1/1:22,30:53:9:73,9,0 1/1:18,35:53:6:56,6,0 1/1:26,34:60:18:188,18,0 0/1:23,17:40:16:16,0,61
HE670865 368612 . A . 6.62 LowQual AN=16;DP=409;MQ=42.14;MQ0=12 GT:DP 0/0:27 0/0:35 0/0:36 0/0:26 0/0:49 0/0:47 0/0:55 0/0:31
HE670865 368613 . G A 325 . AC=11;AF=0.917;AN=12;BaseQRankSum=3.200;DP=404;Dels=0.46;FS=1.094;HaplotypeScore=347.5139;MLEAC=11;MLEAF=0.917;MQ=41.95;MQ0=12;MQRankSum=1.300;QD=1.02;ReadPosRankSum=1.782 GT:AD:DP:GQ:PL 1/1:4,16:21:3:22,3,0 ./. 1/1:0,25:26:3:23,3,0 1/1:1,15:21:3:26,3,0 0/1:4,28:36:10:62,0,10 1/1:1,35:36:6:52,6,0 1/1:2,35:39:15:156,15,0 ./.
HE670865 368614 . A . 17.30 LowQual AN=12;DP=416;MQ=41.75;MQ0=12 GT:DP ./. 0/0:25 0/0:36 ./. 0/0:47 0/0:41 0/0:50 0/0:25
HE670865 368615 . A . 16.42 LowQual AN=12;DP=411;MQ=41.66;MQ0=11 GT:DP ./. 0/0:25 0/0:35 ./. 0/0:45 0/0:40 0/0:49 0/0:22
HE670865 368616 . A . 17.42 LowQual AN=10;DP=406;MQ=41.71;MQ0=10 GT:DP ./. 0/0:24 ./. ./. 0/0:40 0/0:40 0/0:47 0/0:22
HE670865 368617 . T C 94.41 . AC=4;AF=0.667;AN=6;BaseQRankSum=0.262;DP=385;Dels=0.38;FS=7.489;HaplotypeScore=362.4384;MLEAC=5;MLEAF=0.833;MQ=41.72;MQ0=9;MQRankSum=1.775;QD=0.84;ReadPosRankSum=4.850 GT:AD:DP:GQ:PL ./. 0/0:15,1:20:3:0,3,32 ./. ./. ./. 1/1:5,31:36:9:85,9,0 1/1:14,28:42:3:26,3,0 ./.
HE670865 368618 . C . 17.27 LowQual AN=6;DP=370;MQ=41.65;MQ0=8 GT:DP ./. 0/0:10 ./. ./. ./. 0/0:34 0/0:37 ./.
HE670865 368619 . G A 35.14 . AC=4;AF=0.500;AN=8;BaseQRankSum=-3.320;DP=365;Dels=0.54;FS=0.000;HaplotypeScore=399.7009;MLEAC=4;MLEAF=0.500;MQ=41.71;MQ0=8;MQRankSum=-1.237;QD=0.35;ReadPosRankSum=-4.879 GT:AD:DP:GQ:PL ./. 1/1:0,4:8:3:28,3,0 ./. ./. 0/0:17,3:22:3:0,3,23 0/0:31,3:34:3:0,3,22 1/1:28,8:36:3:25,3,0 ./.
HE670865 368620 . T . 15.82 LowQual AN=2;DP=358;MQ=41.72;MQ0=8 GT:DP ./. 0/0:8 ./. ./. ./. ./. ./. ./.
HE670865 368621 . T . . . DP=358;MQ=41.63;MQ0=8 GT ./. ./. ./. ./. ./. ./. ./. ./.
HE670865 368622 . T . 14.11 LowQual AN=6;DP=358;MQ=41.53;MQ0=8 GT:DP ./. ./. ./. ./. ./. 0/0:11 0/0:19 0/0:17
HE670865 368623 . T . 14.11 LowQual AN=6;DP=365;MQ=41.46;MQ0=7 GT:DP ./. 0/0:24 ./. 0/0:16 ./. 0/0:20 ./. ./.
HE670865 368624 . A . 13.85 LowQual AN=8;DP=366;MQ=41.48;MQ0=7 GT:DP 0/0:22 0/0:26 ./. 0/0:23 ./. 0/0:20 ./. ./.
HE670865 368625 . T . 14.12 LowQual AN=6;DP=353;MQ=41.39;MQ0=7 GT:DP 0/0:20 0/0:26 ./. ./. ./. 0/0:18 ./. ./.
HE670865 368626 . T . . . DP=348;MQ=41.31;MQ0=7 GT ./. ./. ./. ./. ./. ./. ./. ./.
HE670865 368627 . T A 15.18 LowQual AC=2;AF=1.00;AN=2;BaseQRankSum=-1.221;DP=348;Dels=0.55;FS=0.000;HaplotypeScore=403.2395;MLEAC=2;MLEAF=1.00;MQ=41.29;MQ0=7;MQRankSum=-0.971;QD=0.23;ReadPosRankSum=-0.344 GT:AD:DP:GQ:PL ./. ./. ./. ./. 1/1:30,4:34:3:26,3,0 ./. ./. ./.
HE670865 368628 . T . 13.85 LowQual AN=8;DP=357;MQ=41.07;MQ0=8 GT:DP ./. 0/0:22 ./. 0/0:10 0/0:33 0/0:15 ./. ./.
HE670865 368629 . A . 15.14 LowQual AN=10;DP=357;MQ=40.99;MQ0=8 GT:DP 0/0:14 0/0:22 ./. ./. 0/0:34 0/0:15 0/0:21 ./.
HE670865 368630 . T . 18.75 LowQual AN=10;DP=358;MQ=40.89;MQ0=8 GT:DP 0/0:15 0/0:23 ./. 0/0:12 ./. 0/0:16 0/0:21 ./.
HE670865 368631 . C . 18.33 LowQual AN=14;DP=359;MQ=40.70;MQ0=8 GT:DP 0/0:15 0/0:23 0/0:14 ./. 0/0:30 0/0:16 0/0:20 0/0:33
HE670865 368632 . G . 17.87 LowQual AN=16;DP=361;MQ=40.63;MQ0=9 GT:DP 0/0:15 0/0:23 0/0:15 0/0:12 0/0:29 0/0:16 0/0:21 0/0:33
HE670865 368633 . T . 18.90 LowQual AN=14;DP=364;MQ=40.51;MQ0=10 GT:DP 0/0:17 0/0:24 0/0:16 ./. 0/0:32 0/0:17 0/0:22 0/0:34
HE670865 368634 . A . 17 LowQual AN=14;DP=369;MQ=40.67;MQ0=10 GT:DP 0/0:18 0/0:24 0/0:18 ./. 0/0:36 0/0:19 0/0:22 0/0:34
HE670865 368635 . T . 17.22 LowQual AN=14;DP=364;MQ=40.66;MQ0=10 GT:DP 0/0:18 0/0:24 0/0:22 ./. 0/0:37 0/0:20 0/0:21 0/0:32
HE670865 368636 . C . 16.03 LowQual AN=14;DP=362;MQ=40.70;MQ0=10 GT:DP 0/0:16 ./. 0/0:15 0/0:14 0/0:36 0/0:20 0/0:22 0/0:29
HE670865 368637 . A . 0.02 LowQual AN=8;DP=356;MQ=40.91;MQ0=11 GT:DP 0/0:26 ./. 0/0:49 0/0:31 0/0:55 ./. ./. ./.
HE670865 368638 . T G 24.10 LowQual AC=3;AF=0.250;AN=12;BaseQRankSum=-6.075;DP=355;Dels=0.14;FS=8.612;HaplotypeScore=349.9679;MLEAC=3;MLEAF=0.250;MQ=40.80;MQ0=12;MQRankSum=-2.113;QD=0.24;ReadPosRankSum=1.216 GT:AD:DP:GQ:PL ./. 0/0:22,3:25:6:0,6,44 0/0:46,3:49:6:0,6,44 0/0:29,1:30:6:0,6,49 0/0:46,8:54:3:0,3,22 1/1:40,2:42:3:29,3,0 0/1:47,5:52:17:18,0,17 ./.
HE670865 368639 . A . 18.78 LowQual AN=16;DP=352;MQ=40.91;MQ0=13 GT:DP 0/0:25 0/0:26 0/0:48 0/0:30 0/0:52 0/0:42 0/0:50 0/0:26
HE670865 368640 . T . 20.94 LowQual AN=16;DP=351;MQ=40.93;MQ0=13 GT:DP 0/0:25 0/0:26 0/0:48 0/0:30 0/0:51 0/0:42 0/0:50 0/0:26
HE670865 368641 . C . 23.55 LowQual AN=16;DP=347;MQ=40.91;MQ0=13 GT:DP 0/0:24 0/0:25 0/0:49 0/0:26 0/0:45 0/0:39 0/0:47 0/0:18
HE670865 368642 . G T 24 LowQual AC=4;AF=0.250;AN=16;BaseQRankSum=-1.989;DP=349;Dels=0.21;FS=7.313;HaplotypeScore=259.1726;MLEAC=3;MLEAF=0.188;MQ=40.80;MQ0=13;MQRankSum=-2.006;QD=0.14;ReadPosRankSum=1.275 GT:AD:DP:GQ:PL 0/1:21,3:24:11:11,0,88 0/0:23,3:26:18:0,18,152 0/0:47,2:49:30:0,30,252 0/0:25,0:25:30:0,30,256 0/0:38,7:45:18:0,18,164 0/1:37,2:39:27:27,0,177 0/1:43,4:47:2:2,0,144 0/1:13,6:19:13:13,0,79
HE670865 368643 . T . 27.64 LowQual AN=16;DP=345;MQ=40.82;MQ0=13 GT:DP 0/0:24 0/0:26 0/0:49 0/0:26 0/0:44 0/0:38 0/0:46 0/0:20
HE670865 368644 . A . 26.97 LowQual AN=14;DP=343;MQ=40.74;MQ0=13 GT:DP ./. 0/0:22 0/0:47 0/0:26 0/0:36 0/0:37 0/0:42 0/0:17
HE670865 368645 . T . 20.86 LowQual AN=14;DP=340;MQ=40.85;MQ0=12 GT:DP ./. 0/0:21 0/0:47 0/0:26 0/0:36 0/0:36 0/0:41 0/0:17
HE670865 368646 . G C,T 237.20 . AC=5,7;AF=0.417,0.583;AN=12;BaseQRankSum=-1.775;DP=353;Dels=0.25;FS=0.000;HaplotypeScore=196.6140;MLEAC=5,6;MLEAF=0.417,0.500;MQ=40.93;MQ0=12;MQRankSum=-0.594;QD=0.89;ReadPosRankSum=1.392 GT:AD:DP:GQ:PL ./. 2/2:13,14,7:34:3:30,30,30,3,3,0 ./. 1/2:2,13,5:20:15:80,24,15,56,0,53 1/2:10,27,9:46:14:64,20,14,44,0,41 1/2:4,32,2:38:18:61,24,18,37,0,34 2/2:8,30,8:46:9:77,77,77,9,9,0 1/1:18,11,1:30:6:43,6,0,43,6,43
HE670865 368647 . T . 29.30 LowQual AN=16;DP=368;MQ=40.74;MQ0=12 GT:DP 0/0:31 0/0:37 0/0:54 0/0:36 0/0:52 0/0:42 0/0:50 0/0:34
HE670865 368648 . G . 51.60 . AN=16;DP=366;MQ=40.87;MQ0=12 GT:DP 0/0:34 0/0:40 0/0:57 0/0:37 0/0:59 0/0:44 0/0:53 0/0:42

These are the comparable lines from the outputted VCF when using EMIT_ALL_CONFIDENT_SITES:

HE670865 368589 . A . 46.89 . AN=16;DP=441;MQ=43.26;MQ0=13 GT:DP 0/0:35 0/0:45 0/0:58 0/0:38 0/0:75 0/0:64 0/0:64 0/0:47
HE670865 368590 . A . 50.15 . AN=16;DP=447;MQ=43.41;MQ0=13 GT:DP 0/0:35 0/0:45 0/0:60 0/0:39 0/0:76 0/0:66 0/0:65 0/0:47
HE670865 368591 . A . 39.12 . AN=16;DP=448;MQ=43.32;MQ0=12 GT:DP 0/0:36 0/0:45 0/0:61 0/0:39 0/0:77 0/0:66 0/0:66 0/0:46
HE670865 368648 . G . 51.60 . AN=16;DP=366;MQ=40.87;MQ0=12 GT:DP 0/0:34 0/0:40 0/0:57 0/0:37 0/0:59 0/0:44 0/0:53 0/0:42
HE670865 368649 . T . 75.72 . AN=16;DP=368;MQ=40.89;MQ0=12 GT:DP 0/0:34 0/0:41 0/0:57 0/0:38 0/0:59 0/0:44 0/0:53 0/0:42
HE670865 368650 . T . 72.74 . AN=16;DP=369;MQ=40.87;MQ0=12 GT:DP 0/0:34 0/0:41 0/0:57 0/0:38 0/0:59 0/0:44 0/0:53 0/0:43
HE670865 368651 . T . 87.71 . AN=16;DP=369;MQ=40.88;MQ0=12 GT:DP 0/0:33 0/0:42 0/0:58 0/0:37 0/0:59 0/0:44 0/0:53 0/0:43
HE670865 368652 . T . 93.72 . AN=16;DP=373;MQ=40.89;MQ0=12 GT:DP 0/0:32 0/0:43 0/0:59 0/0:37 0/0:59 0/0:45 0/0:54 0/0:43

Any help or insight would be greatly appreciated.

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @Jakemorris89
    Hi,

    I suspect the -nt 4 is causing some differences. Can you try running without -nt?

    -Sheila

  • Hi, thanks.

    I tried running it with EMIT_ALL_CONFIDENT_SITES and without the -nt 4 argument and it still gives me the same output I had before which excludes sites which seem to be 'good':

    HE670865 368589 . A . 46.89 . AN=16;DP=441;MQ=43.26;MQ0=13 GT:DP 0/0:35 0/0:45 0/0:58 0/0:38 0/0:75 0/0:64 0/0:64 0/0:47
    HE670865 368590 . A . 50.15 . AN=16;DP=447;MQ=43.41;MQ0=13 GT:DP 0/0:35 0/0:45 0/0:60 0/0:39 0/0:76 0/0:66 0/0:65 0/0:47
    HE670865 368591 . A . 39.12 . AN=16;DP=448;MQ=43.32;MQ0=12 GT:DP 0/0:36 0/0:45 0/0:61 0/0:39 0/0:77 0/0:66 0/0:66 0/0:46
    HE670865 368648 . G . 51.60 . AN=16;DP=366;MQ=40.87;MQ0=12 GT:DP 0/0:34 0/0:40 0/0:57 0/0:37 0/0:59 0/0:44 0/0:53 0/0:42
    HE670865 368649 . T . 75.72 . AN=16;DP=368;MQ=40.89;MQ0=12 GT:DP 0/0:34 0/0:41 0/0:57 0/0:38 0/0:59 0/0:44 0/0:53 0/0:42
    HE670865 368650 . T . 72.74 . AN=16;DP=369;MQ=40.87;MQ0=12 GT:DP 0/0:34 0/0:41 0/0:57 0/0:38 0/0:59 0/0:44 0/0:53 0/0:43
    HE670865 368651 . T . 87.71 . AN=16;DP=369;MQ=40.88;MQ0=12 GT:DP 0/0:33 0/0:42 0/0:58 0/0:37 0/0:59 0/0:44 0/0:53 0/0:43
    HE670865 368652 . T . 93.72 . AN=16;DP=373;MQ=40.89;MQ0=12 GT:DP 0/0:32 0/0:43 0/0:59 0/0:37 0/0:59 0/0:45 0/0:54 0/0:43

  • tommycarstensentommycarstensen United KingdomMember

    @Jakemorris89 What version are you on? Could you try setting --standard_min_confidence_threshold_for_calling and --standard_min_confidence_threshold_for_emitting to zero and see if the missing variants are brought back to life? Odd...

  • Jakemorris89Jakemorris89 UKMember
    edited August 2015

    Hi @tommycarstensen,

    I was using version 2.7.2. I just tried running it with --standard_min_confidence_threshold_for_calling and --standard_min_confidence_threshold_for_emitting set to 0, but it still doesn't output those sites so I still just get the same output:

    HE670865 368589 . A . 46.89 . AN=16;DP=441;MQ=43.26;MQ0=13 GT:DP 0/0:35 0/0:45 0/0:58 0/0:38 0/0:75 0/0:64 0/0:64 0/0:47
    HE670865 368590 . A . 50.15 . AN=16;DP=447;MQ=43.41;MQ0=13 GT:DP 0/0:35 0/0:45 0/0:60 0/0:39 0/0:76 0/0:66 0/0:65 0/0:47
    HE670865 368591 . A . 39.12 . AN=16;DP=448;MQ=43.32;MQ0=12 GT:DP 0/0:36 0/0:45 0/0:61 0/0:39 0/0:77 0/0:66 0/0:66 0/0:46
    HE670865 368648 . G . 51.60 . AN=16;DP=366;MQ=40.87;MQ0=12 GT:DP 0/0:34 0/0:40 0/0:57 0/0:37 0/0:59 0/0:44 0/0:53 0/0:42
    HE670865 368649 . T . 75.72 . AN=16;DP=368;MQ=40.89;MQ0=12 GT:DP 0/0:34 0/0:41 0/0:57 0/0:38 0/0:59 0/0:44 0/0:53 0/0:42
    HE670865 368650 . T . 72.74 . AN=16;DP=369;MQ=40.87;MQ0=12 GT:DP 0/0:34 0/0:41 0/0:57 0/0:38 0/0:59 0/0:44 0/0:53 0/0:43
    HE670865 368651 . T . 87.71 . AN=16;DP=369;MQ=40.88;MQ0=12 GT:DP 0/0:33 0/0:42 0/0:58 0/0:37 0/0:59 0/0:44 0/0:53 0/0:43
    HE670865 368652 . T . 93.72 . AN=16;DP=373;MQ=40.89;MQ0=12 GT:DP 0/0:32 0/0:43 0/0:59 0/0:37 0/0:59 0/0:45 0/0:54 0/0:43

    rather than the output I had when using EMIT_ALL_SITES minus the LowQual sites, which I agree seems odd.

    I also for good measure ran it with the newest version 3.4-46, again it didn't emit sme of the the sites between HE670865 368589 and HE670865 368652 which were not LowQual, I just had the same output as above.

    I've pasted the full command line options from the top of each VCF below. Perhaps you can see something strange, that I can't.

    For EMIT_ALL_SITES:

    GATKCommandLine=<ID=UnifiedGenotyper,Version=2.7-2-g6bda569,Date="Fri Aug 14 01:12:43 BST 2015",Epoch=1439511163074,CommandLineOptions="analysis_type=UnifiedGenotyper input_file=[all_neruda_scaffmerged_guiana_only.realign.bam] read_buffer_size=null phone_home=AWS gatk_key=null tag=NA read_filter=[] intervals=null excludeIntervals=null interval_set_rule=UNION interval_merging=ALL interval_padding=0 reference_sequence=/data/genome_references/Hmel1-1_primaryScaffolds_mtDNA.fasta nonDeterministicRandomSeed=false disableDithering=false maxRuntime=-1 maxRuntimeUnits=MINUTES downsampling_type=BY_SAMPLE downsample_to_fraction=null downsample_to_coverage=250 baq=CALCULATE_AS_NECESSARY baqGapOpenPenalty=40.0 fix_misencoded_quality_scores=false allow_potentially_misencoded_quality_scores=false useOriginalQualities=false defaultBaseQualities=-1 performanceLog=null BQSR=null quantize_quals=0 disable_indel_quals=false emit_original_quals=false preserve_qscores_less_than=6 globalQScorePrior=-1.0 allow_bqsr_on_reduced_bams_despite_repeated_warnings=false validation_strictness=SILENT remove_program_records=false keep_program_records=false sample_rename_mapping_file=null unsafe=null disable_auto_index_creation_and_locking_when_reading_rods=false num_threads=4 num_cpu_threads_per_data_thread=1 num_io_threads=0 monitorThreadEfficiency=false num_bam_file_handles=null read_group_black_list=null pedigree=[] pedigreeString=[] pedigreeValidationType=STRICT allow_intervals_with_unindexed_bam=false generateShadowBCF=false logging_level=INFO log_to_file=null help=false version=false genotype_likelihoods_model=SNP pcr_error_rate=1.0E-4 computeSLOD=false annotateNDA=false pair_hmm_implementation=LOGLESS_CACHING min_base_quality_score=17 max_deletion_fraction=0.05 allSitePLs=false min_indel_count_for_genotyping=5 min_indel_fraction_per_sample=0.25 indelGapContinuationPenalty=10 indelGapOpenPenalty=45 indelHaplotypeSize=80 indelDebug=false ignoreSNPAlleles=false allReadsSP=false ignoreLaneInfo=false reference_sample_calls=(RodBinding name= source=UNBOUND) reference_sample_name=null sample_ploidy=2 min_quality_score=1 max_quality_score=40 site_quality_prior=20 min_power_threshold_for_calling=0.95 min_reference_depth=100 exclude_filtered_reference_sites=false output_mode=EMIT_ALL_SITES heterozygosity=0.015 indel_heterozygosity=1.25E-4 genotyping_mode=DISCOVERY standard_min_confidence_threshold_for_calling=30.0 standard_min_confidence_threshold_for_emitting=30.0 alleles=(RodBinding name= source=UNBOUND) max_alternate_alleles=6 input_prior=[] contamination_fraction_to_filter=0.0 contamination_fraction_per_sample_file=null p_nonref_model=EXACT_INDEPENDENT exactcallslog=null dbsnp=(RodBinding name= source=UNBOUND) comp=[] out=org.broadinstitute.sting.gatk.io.stubs.VariantContextWriterStub no_cmdline_in_header=org.broadinstitute.sting.gatk.io.stubs.VariantContextWriterStub sites_only=org.broadinstitute.sting.gatk.io.stubs.VariantContextWriterStub bcf=org.broadinstitute.sting.gatk.io.stubs.VariantContextWriterStub onlyEmitSamples=[] debug_file=null metrics_file=null annotation=[] excludeAnnotation=[] filter_reads_with_N_cigar=false filter_mismatching_base_and_quals=false filter_bases_not_stored=false">

    For EMIT_ALL_CONFIDENT_SITES, threholds set to 0:

    GATKCommandLine=<ID=UnifiedGenotyper,Version=2.7-2-g6bda569,Date="Wed Aug 19 14:44:55 BST 2015",Epoch=1439991895568,CommandLineOptions="analysis_type=UnifiedGenotyper input_file=[all_neruda_scaffmerged_guiana_only.realign.bam] read_buffer_size=null phone_home=AWS gatk_key=null tag=NA read_filter=[] intervals=null excludeIntervals=null interval_set_rule=UNION interval_merging=ALL interval_padding=0 reference_sequence=/data/genome_references/Hmel1-1_primaryScaffolds_mtDNA.fasta nonDeterministicRandomSeed=false disableDithering=false maxRuntime=-1 maxRuntimeUnits=MINUTES downsampling_type=BY_SAMPLE downsample_to_fraction=null downsample_to_coverage=250 baq=CALCULATE_AS_NECESSARY baqGapOpenPenalty=40.0 fix_misencoded_quality_scores=false allow_potentially_misencoded_quality_scores=false useOriginalQualities=false defaultBaseQualities=-1 performanceLog=null BQSR=null quantize_quals=0 disable_indel_quals=false emit_original_quals=false preserve_qscores_less_than=6 globalQScorePrior=-1.0 allow_bqsr_on_reduced_bams_despite_repeated_warnings=false validation_strictness=SILENT remove_program_records=false keep_program_records=false sample_rename_mapping_file=null unsafe=null disable_auto_index_creation_and_locking_when_reading_rods=false num_threads=1 num_cpu_threads_per_data_thread=1 num_io_threads=0 monitorThreadEfficiency=false num_bam_file_handles=null read_group_black_list=null pedigree=[] pedigreeString=[] pedigreeValidationType=STRICT allow_intervals_with_unindexed_bam=false generateShadowBCF=false logging_level=INFO log_to_file=null help=false version=false genotype_likelihoods_model=SNP pcr_error_rate=1.0E-4 computeSLOD=false annotateNDA=false pair_hmm_implementation=LOGLESS_CACHING min_base_quality_score=17 max_deletion_fraction=0.05 allSitePLs=false min_indel_count_for_genotyping=5 min_indel_fraction_per_sample=0.25 indelGapContinuationPenalty=10 indelGapOpenPenalty=45 indelHaplotypeSize=80 indelDebug=false ignoreSNPAlleles=false allReadsSP=false ignoreLaneInfo=false reference_sample_calls=(RodBinding name= source=UNBOUND) reference_sample_name=null sample_ploidy=2 min_quality_score=1 max_quality_score=40 site_quality_prior=20 min_power_threshold_for_calling=0.95 min_reference_depth=100 exclude_filtered_reference_sites=false output_mode=EMIT_ALL_CONFIDENT_SITES heterozygosity=0.015 indel_heterozygosity=1.25E-4 genotyping_mode=DISCOVERY standard_min_confidence_threshold_for_calling=0.0 standard_min_confidence_threshold_for_emitting=0.0 alleles=(RodBinding name= source=UNBOUND) max_alternate_alleles=6 input_prior=[] contamination_fraction_to_filter=0.0 contamination_fraction_per_sample_file=null p_nonref_model=EXACT_INDEPENDENT exactcallslog=null dbsnp=(RodBinding name= source=UNBOUND) comp=[] out=org.broadinstitute.sting.gatk.io.stubs.VariantContextWriterStub no_cmdline_in_header=org.broadinstitute.sting.gatk.io.stubs.VariantContextWriterStub sites_only=org.broadinstitute.sting.gatk.io.stubs.VariantContextWriterStub bcf=org.broadinstitute.sting.gatk.io.stubs.VariantContextWriterStub onlyEmitSamples=[] debug_file=null metrics_file=null annotation=[] excludeAnnotation=[] filter_reads_with_N_cigar=false filter_mismatching_base_and_quals=false filter_bases_not_stored=false">

    For EMIT_ALL_CONFIDENT_SITES, new version, default thresholds:

    GATKCommandLine.UnifiedGenotyper=<ID=UnifiedGenotyper,Version=3.4-46-gbc02625,Date="Wed Aug 19 14:58:46 BST 2015",Epoch=1439992726651,CommandLineOptions="analysis_type=UnifiedGenotyper input_file=[all_neruda_scaffmerged_guiana_only.realign.bam] showFullBamList=false read_buffer_size=null phone_home=AWS gatk_key=null tag=NA read_filter=[] disable_read_filter=[] intervals=null excludeIntervals=null interval_set_rule=UNION interval_merging=ALL interval_padding=0 reference_sequence=/data/genome_references/Hmel1-1_primaryScaffolds_mtDNA.fasta nonDeterministicRandomSeed=false disableDithering=false maxRuntime=-1 maxRuntimeUnits=MINUTES downsampling_type=BY_SAMPLE downsample_to_fraction=null downsample_to_coverage=250 baq=CALCULATE_AS_NECESSARY baqGapOpenPenalty=40.0 refactor_NDN_cigar_string=false fix_misencoded_quality_scores=false allow_potentially_misencoded_quality_scores=false useOriginalQualities=false defaultBaseQualities=-1 performanceLog=null BQSR=null quantize_quals=0 disable_indel_quals=false emit_original_quals=false preserve_qscores_less_than=6 globalQScorePrior=-1.0 validation_strictness=SILENT remove_program_records=false keep_program_records=false sample_rename_mapping_file=null unsafe=null disable_auto_index_creation_and_locking_when_reading_rods=false no_cmdline_in_header=false sites_only=false never_trim_vcf_format_field=false bcf=false bam_compression=null simplifyBAM=false disable_bam_indexing=false generate_md5=false num_threads=1 num_cpu_threads_per_data_thread=1 num_io_threads=0 monitorThreadEfficiency=false num_bam_file_handles=null read_group_black_list=null pedigree=[] pedigreeString=[] pedigreeValidationType=STRICT allow_intervals_with_unindexed_bam=false generateShadowBCF=false variant_index_type=DYNAMIC_SEEK variant_index_parameter=-1 logging_level=INFO log_to_file=null help=false version=false genotype_likelihoods_model=SNP pcr_error_rate=1.0E-4 computeSLOD=false pair_hmm_implementation=LOGLESS_CACHING min_base_quality_score=17 max_deletion_fraction=0.05 min_indel_count_for_genotyping=5 min_indel_fraction_per_sample=0.25 indelGapContinuationPenalty=10 indelGapOpenPenalty=45 indelHaplotypeSize=80 indelDebug=false ignoreSNPAlleles=false allReadsSP=false ignoreLaneInfo=false reference_sample_calls=(RodBinding name= source=UNBOUND) reference_sample_name=null min_quality_score=1 max_quality_score=40 site_quality_prior=20 min_power_threshold_for_calling=0.95 annotateNDA=false heterozygosity=0.015 indel_heterozygosity=1.25E-4 standard_min_confidence_threshold_for_calling=30.0 standard_min_confidence_threshold_for_emitting=30.0 max_alternate_alleles=6 input_prior=[] sample_ploidy=2 genotyping_mode=DISCOVERY alleles=(RodBinding name= source=UNBOUND) contamination_fraction_to_filter=0.0 contamination_fraction_per_sample_file=null p_nonref_model=null exactcallslog=null output_mode=EMIT_ALL_CONFIDENT_SITES allSitePLs=false dbsnp=(RodBinding name= source=UNBOUND) comp=[] out=/data/jake_morris/neruda_demeter_analysis/mapping_Neruda_to_melpref/divergentscaff_outG_analysis/all_neruda_scaffmerged_guiana_only_confidentsites_NewV.vcf onlyEmitSamples=[] debug_file=null metrics_file=null annotation=[] excludeAnnotation=[] filter_reads_with_N_cigar=false filter_mismatching_base_and_quals=false filter_bases_not_stored=false">

    Thanks,
    Jake

  • tommycarstensentommycarstensen United KingdomMember

    @Jakemorris89 The current version of GATK is 3.4-0 / 3.4-46. Version 2.7.2 of GATK is like running Windows 3.1. Try with the latest version of GATK. Like Windows 95 the Broad no longer supports version 2.7.2. Best of luck.

  • Hi @tommycarstensen @Sheila,

    I've now rerun with -nt 0 and with --standard_min_confidence_threshold_for_calling and --standard_min_confidence_threshold_for_emitting to zero and with the latest version of GATK (3.4-46). Unfortunately the problem persists just as it did originally. EMIT_ALL_CONFIDENT_SITES still seems to be missing out sites that seem 'good' when they're emitted by EMIT_ALL_SITES.

    Again here are the full command line options from the top of each VCF below:

    Note: these are the runs with default call and emit thresholds, but I tried it with them set to 0 too and same thing, besides I want this to work with thresholds set at default to remove LowQual sites, (otherwise I could just use EMIT_ALL sites) I just don't know why EMIT_ALL_CONFIDENT_SITES is removing what seem like confident/'good' sites too.

    Emit all sites:

    GATKCommandLine.UnifiedGenotyper=<ID=UnifiedGenotyper,Version=3.4-46-gbc02625,Date="Wed Aug 19 18:26:23 BST 2015",Epoch=1440005183113,CommandLineOptions="analysis_type=UnifiedGenotyper input_file=[all_neruda_scaffmerged_guiana_only.realign.bam] showFullBamList=false read_buffer_size=null phone_home=AWS gatk_key=null tag=NA read_filter=[] disable_read_filter=[] intervals=null excludeIntervals=null interval_set_rule=UNION interval_merging=ALL interval_padding=0 reference_sequence=/data/genome_references/Hmel1-1_primaryScaffolds_mtDNA.fasta nonDeterministicRandomSeed=false disableDithering=false maxRuntime=-1 maxRuntimeUnits=MINUTES downsampling_type=BY_SAMPLE downsample_to_fraction=null downsample_to_coverage=250 baq=CALCULATE_AS_NECESSARY baqGapOpenPenalty=40.0 refactor_NDN_cigar_string=false fix_misencoded_quality_scores=false allow_potentially_misencoded_quality_scores=false useOriginalQualities=false defaultBaseQualities=-1 performanceLog=null BQSR=null quantize_quals=0 disable_indel_quals=false emit_original_quals=false preserve_qscores_less_than=6 globalQScorePrior=-1.0 validation_strictness=SILENT remove_program_records=false keep_program_records=false sample_rename_mapping_file=null unsafe=null disable_auto_index_creation_and_locking_when_reading_rods=false no_cmdline_in_header=false sites_only=false never_trim_vcf_format_field=false bcf=false bam_compression=null simplifyBAM=false disable_bam_indexing=false generate_md5=false num_threads=1 num_cpu_threads_per_data_thread=1 num_io_threads=0 monitorThreadEfficiency=false num_bam_file_handles=null read_group_black_list=null pedigree=[] pedigreeString=[] pedigreeValidationType=STRICT allow_intervals_with_unindexed_bam=false generateShadowBCF=false variant_index_type=DYNAMIC_SEEK variant_index_parameter=-1 logging_level=INFO log_to_file=null help=false version=false genotype_likelihoods_model=SNP pcr_error_rate=1.0E-4 computeSLOD=false pair_hmm_implementation=LOGLESS_CACHING min_base_quality_score=17 max_deletion_fraction=0.05 min_indel_count_for_genotyping=5 min_indel_fraction_per_sample=0.25 indelGapContinuationPenalty=10 indelGapOpenPenalty=45 indelHaplotypeSize=80 indelDebug=false ignoreSNPAlleles=false allReadsSP=false ignoreLaneInfo=false reference_sample_calls=(RodBinding name= source=UNBOUND) reference_sample_name=null min_quality_score=1 max_quality_score=40 site_quality_prior=20 min_power_threshold_for_calling=0.95 annotateNDA=false heterozygosity=0.015 indel_heterozygosity=1.25E-4 standard_min_confidence_threshold_for_calling=30.0 standard_min_confidence_threshold_for_emitting=30.0 max_alternate_alleles=6 input_prior=[] sample_ploidy=2 genotyping_mode=DISCOVERY alleles=(RodBinding name= source=UNBOUND) contamination_fraction_to_filter=0.0 contamination_fraction_per_sample_file=null p_nonref_model=null exactcallslog=null output_mode=EMIT_ALL_SITES allSitePLs=false dbsnp=(RodBinding name= source=UNBOUND) comp=[] out=/data/jake_morris/neruda_demeter_analysis/mapping_Neruda_to_melpref/divergentscaff_outG_analysis/all_neruda_scaffmerged_guiana_only_confidentsites_NewV_emitall.vcf onlyEmitSamples=[] debug_file=null metrics_file=null annotation=[] excludeAnnotation=[] filter_reads_with_N_cigar=false filter_mismatching_base_and_quals=false filter_bases_not_stored=false">

    EMIT_ALL_CONFIDENT_SITES:

    GATKCommandLine.UnifiedGenotyper=<ID=UnifiedGenotyper,Version=3.4-46-gbc02625,Date="Wed Aug 19 14:58:46 BST 2015",Epoch=1439992726651,CommandLineOptions="analysis_type=UnifiedGenotyper input_file=[all_neruda_scaffmerged_guiana_only.realign.bam] showFullBamList=false read_buffer_size=null phone_home=AWS gatk_key=null tag=NA read_filter=[] disable_read_filter=[] intervals=null excludeIntervals=null interval_set_rule=UNION interval_merging=ALL interval_padding=0 reference_sequence=/data/genome_references/Hmel1-1_primaryScaffolds_mtDNA.fasta nonDeterministicRandomSeed=false disableDithering=false maxRuntime=-1 maxRuntimeUnits=MINUTES downsampling_type=BY_SAMPLE downsample_to_fraction=null downsample_to_coverage=250 baq=CALCULATE_AS_NECESSARY baqGapOpenPenalty=40.0 refactor_NDN_cigar_string=false fix_misencoded_quality_scores=false allow_potentially_misencoded_quality_scores=false useOriginalQualities=false defaultBaseQualities=-1 performanceLog=null BQSR=null quantize_quals=0 disable_indel_quals=false emit_original_quals=false preserve_qscores_less_than=6 globalQScorePrior=-1.0 validation_strictness=SILENT remove_program_records=false keep_program_records=false sample_rename_mapping_file=null unsafe=null disable_auto_index_creation_and_locking_when_reading_rods=false no_cmdline_in_header=false sites_only=false never_trim_vcf_format_field=false bcf=false bam_compression=null simplifyBAM=false disable_bam_indexing=false generate_md5=false num_threads=1 num_cpu_threads_per_data_thread=1 num_io_threads=0 monitorThreadEfficiency=false num_bam_file_handles=null read_group_black_list=null pedigree=[] pedigreeString=[] pedigreeValidationType=STRICT allow_intervals_with_unindexed_bam=false generateShadowBCF=false variant_index_type=DYNAMIC_SEEK variant_index_parameter=-1 logging_level=INFO log_to_file=null help=false version=false genotype_likelihoods_model=SNP pcr_error_rate=1.0E-4 computeSLOD=false pair_hmm_implementation=LOGLESS_CACHING min_base_quality_score=17 max_deletion_fraction=0.05 min_indel_count_for_genotyping=5 min_indel_fraction_per_sample=0.25 indelGapContinuationPenalty=10 indelGapOpenPenalty=45 indelHaplotypeSize=80 indelDebug=false ignoreSNPAlleles=false allReadsSP=false ignoreLaneInfo=false reference_sample_calls=(RodBinding name= source=UNBOUND) reference_sample_name=null min_quality_score=1 max_quality_score=40 site_quality_prior=20 min_power_threshold_for_calling=0.95 annotateNDA=false heterozygosity=0.015 indel_heterozygosity=1.25E-4 standard_min_confidence_threshold_for_calling=30.0 standard_min_confidence_threshold_for_emitting=30.0 max_alternate_alleles=6 input_prior=[] sample_ploidy=2 genotyping_mode=DISCOVERY alleles=(RodBinding name= source=UNBOUND) contamination_fraction_to_filter=0.0 contamination_fraction_per_sample_file=null p_nonref_model=null exactcallslog=null output_mode=EMIT_ALL_CONFIDENT_SITES allSitePLs=false dbsnp=(RodBinding name= source=UNBOUND) comp=[] out=/data/jake_morris/neruda_demeter_analysis/mapping_Neruda_to_melpref/divergentscaff_outG_analysis/all_neruda_scaffmerged_guiana_only_confidentsites_NewV.vcf onlyEmitSamples=[] debug_file=null metrics_file=null annotation=[] excludeAnnotation=[] filter_reads_with_N_cigar=false filter_mismatching_base_and_quals=false filter_bases_not_stored=false">

    The VCF output from sites HE670865 368589 to HE670865 368652 are still more or less the same as before:

    i.e for EMIT_ALL_SITES:

    HE670865 368589 . A . 46.89 . AN=16;DP=441;MQ=43.26;MQ0=13 GT:DP 0/0:35 0/0:45 0/0:58 0/0:38 0/0:75 0/0:64 0/0:64 0/0:47
    HE670865 368590 . A . 50.15 . AN=16;DP=447;MQ=43.41;MQ0=13 GT:DP 0/0:35 0/0:45 0/0:60 0/0:39 0/0:76 0/0:66 0/0:65 0/0:47
    HE670865 368591 . A . 39.12 . AN=16;DP=448;MQ=43.32;MQ0=12 GT:DP 0/0:36 0/0:45 0/0:61 0/0:39 0/0:77 0/0:66 0/0:66 0/0:46
    HE670865 368592 . C T 5437.28 . AC=16;AF=1.00;AN=16;BaseQRankSum=0.446;DP=445;Dels=0.06;FS=0.000;HaplotypeScore=99.1839;MLEAC=16;MLEAF=1.00;MQ=43.48;MQ0=9;MQRankSum=4.288;QD=12.22;ReadPosRankSum=6.227;SOR=1.105 GT:AD:DP:GQ:PL 1/1:1,27:29:27:296,27,0 1/1:0,45:45:30:344,30,0 1/1:6,43:52:33:327,33,01/1:1,38:39:48:506,48,0 1/1:5,68:76:99:1100,102,0 1/1:4,61:65:68:1068,68,0 1/1:0,65:66:99:1035,99,0 1/1:0,46:46:69:774,69,0
    HE670865 368593 . A . 50.76 . AN=16;DP=448;MQ=43.48;MQ0=9 GT:DP 0/0:29 0/0:45 0/0:53 0/0:41 0/0:76 0/0:65 0/0:66 0/0:47
    HE670865 368594 . G . 69.75 . AN=16;DP=449;MQ=43.48;MQ0=9 GT:DP 0/0:28 0/0:46 0/0:51 0/0:41 0/0:76 0/0:66 0/0:66 0/0:47
    HE670865 368595 . A . 66.75 . AN=16;DP=447;MQ=43.43;MQ0=9 GT:DP 0/0:30 0/0:46 0/0:50 0/0:41 0/0:75 0/0:66 0/0:66 0/0:47
    HE670865 368596 . A . 72.75 . AN=16;DP=445;MQ=43.54;MQ0=9 GT:DP 0/0:31 0/0:46 0/0:48 0/0:41 0/0:76 0/0:66 0/0:67 0/0:47
    HE670865 368597 . G . 66.75 . AN=16;DP=442;MQ=43.48;MQ0=9 GT:DP 0/0:25 0/0:46 0/0:41 0/0:40 0/0:74 0/0:67 0/0:67 0/0:47
    HE670865 368598 . A . 75.75 . AN=16;DP=444;MQ=43.44;MQ0=9 GT:DP 0/0:25 0/0:46 0/0:41 0/0:39 0/0:75 0/0:69 0/0:67 0/0:47
    [jm1293@heliconius divergentscaff_outG_analysis]$ head -n50 grep_HE670865_368589_guiana_only_confidentsites_NewV_emitall
    HE670865 368589 . A . 46.89 . AN=16;DP=441;MQ=43.26;MQ0=13 GT:DP 0/0:35 0/0:45 0/0:58 0/0:38 0/0:75 0/0:64 0/0:64 0/0:47
    HE670865 368590 . A . 50.15 . AN=16;DP=447;MQ=43.41;MQ0=13 GT:DP 0/0:35 0/0:45 0/0:60 0/0:39 0/0:76 0/0:66 0/0:65 0/0:47
    HE670865 368591 . A . 39.12 . AN=16;DP=448;MQ=43.32;MQ0=12 GT:DP 0/0:36 0/0:45 0/0:61 0/0:39 0/0:77 0/0:66 0/0:66 0/0:46
    HE670865 368592 . C T 5437.28 . AC=16;AF=1.00;AN=16;BaseQRankSum=0.446;DP=445;Dels=0.06;FS=0.000;HaplotypeScore=99.1839;MLEAC=16;MLEAF=1.00;MQ=43.48;MQ0=9;MQRankSum=4.288;QD=12.22;ReadPosRankSum=6.227;SOR=1.105 GT:AD:DP:GQ:PL 1/1:1,27:29:27:296,27,0 1/1:0,45:45:30:344,30,0 1/1:6,43:52:33:327,33,01/1:1,38:39:48:506,48,0 1/1:5,68:76:99:1100,102,0 1/1:4,61:65:68:1068,68,0 1/1:0,65:66:99:1035,99,0 1/1:0,46:46:69:774,69,0
    HE670865 368593 . A . 50.76 . AN=16;DP=448;MQ=43.48;MQ0=9 GT:DP 0/0:29 0/0:45 0/0:53 0/0:41 0/0:76 0/0:65 0/0:66 0/0:47
    HE670865 368594 . G . 69.75 . AN=16;DP=449;MQ=43.48;MQ0=9 GT:DP 0/0:28 0/0:46 0/0:51 0/0:41 0/0:76 0/0:66 0/0:66 0/0:47
    HE670865 368595 . A . 66.75 . AN=16;DP=447;MQ=43.43;MQ0=9 GT:DP 0/0:30 0/0:46 0/0:50 0/0:41 0/0:75 0/0:66 0/0:66 0/0:47
    HE670865 368596 . A . 72.75 . AN=16;DP=445;MQ=43.54;MQ0=9 GT:DP 0/0:31 0/0:46 0/0:48 0/0:41 0/0:76 0/0:66 0/0:67 0/0:47
    HE670865 368597 . G . 66.75 . AN=16;DP=442;MQ=43.48;MQ0=9 GT:DP 0/0:25 0/0:46 0/0:41 0/0:40 0/0:74 0/0:67 0/0:67 0/0:47
    HE670865 368598 . A . 75.75 . AN=16;DP=444;MQ=43.44;MQ0=9 GT:DP 0/0:25 0/0:46 0/0:41 0/0:39 0/0:75 0/0:69 0/0:67 0/0:47
    HE670865 368599 . A . 72.75 . AN=16;DP=445;MQ=43.49;MQ0=9 GT:DP 0/0:24 0/0:45 0/0:39 0/0:39 0/0:76 0/0:69 0/0:68 0/0:48
    HE670865 368600 . T A 4405.18 . AC=8;AF=0.500;AN=16;BaseQRankSum=2.667;DP=441;Dels=0.08;FS=6.510;HaplotypeScore=100.3222;MLEAC=8;MLEAF=0.500;MQ=43.26;MQ0=9;MQRankSum=0.535;QD=24.75;ReadPosRankSum=1.174;SOR=1.250 GT:AD:DP:GQ:PL 1/1:3,19:22:48:570,48,0 1/1:2,45:47:99:1386,114,0 1/1:4,33:37:96:1196,96,0 1/1:0,38:38:99:1307,105,0 0/0:76,2:78:99:0,178,2196 0/0:65,0:65:99:0,159,2002 0/0:70,0:70:99:0,168,2042 0/0:49,0:49:99:0,129,1549
    HE670865 368601 . T . 63.75 . AN=16;DP=442;MQ=43.37;MQ0=9 GT:DP 0/0:19 0/0:47 0/0:37 0/0:38 0/0:77 0/0:66 0/0:72 0/0:48
    HE670865 368602 . A . 66.75 . AN=16;DP=444;MQ=43.44;MQ0=9 GT:DP 0/0:19 0/0:47 0/0:38 0/0:38 0/0:78 0/0:65 0/0:74 0/0:48
    HE670865 368603 . C . 66.75 . AN=16;DP=441;MQ=43.50;MQ0=11 GT:DP 0/0:20 0/0:47 0/0:39 0/0:37 0/0:76 0/0:64 0/0:73 0/0:48
    HE670865 368604 . A . 63.75 . AN=16;DP=441;MQ=43.24;MQ0=11 GT:DP 0/0:21 0/0:47 0/0:39 0/0:38 0/0:74 0/0:64 0/0:74 0/0:48
    HE670865 368605 . G . 60.75 . AN=16;DP=436;MQ=43.17;MQ0=12 GT:DP 0/0:20 0/0:47 0/0:40 0/0:38 0/0:73 0/0:63 0/0:73 0/0:47
    HE670865 368606 . A . 60.75 . AN=16;DP=436;MQ=43.05;MQ0=12 GT:DP 0/0:22 0/0:45 0/0:40 0/0:39 0/0:73 0/0:63 0/0:75 0/0:46
    HE670865 368607 . T C 3984.32 . AC=15;AF=0.938;AN=16;BaseQRankSum=0.451;DP=433;Dels=0.06;FS=11.987;HaplotypeScore=229.2107;MLEAC=15;MLEAF=0.938;MQ=42.94;MQ0=12;MQRankSum=1.549;QD=9.20;ReadPosRankSum=-0.389;SOR=1.715 GT:AD:DP:GQ:PL 1/1:7,18:26:12:109,12,0 1/1:2,39:44:57:475,57,0 1/1:6,36:44:51:426,51,01/1:1,36:39:39:328,39,0 1/1:6,67:73:44:603,44,0 1/1:0,63:63:99:1016,105,0 1/1:4,69:74:90:817,90,0 0/1:22,23:46:99:232,0,373
    HE670865 368608 . A . 24.45 LowQual AN=16;DP=438;MQ=42.86;MQ0=12 GT:DP 0/0:34 0/0:45 0/0:55 0/0:39 0/0:74 0/0:63 0/0:74 0/0:47
    HE670865 368609 . A . 30.14 . AN=14;DP=434;MQ=42.67;MQ0=12 GT:DP ./. 0/0:44 0/0:53 0/0:39 0/0:67 0/0:61 0/0:74 0/0:47
    HE670865 368610 . A . 20.62 LowQual AN=16;DP=432;MQ=42.47;MQ0=12 GT:DP 0/0:33 0/0:45 0/0:54 0/0:39 0/0:67 0/0:59 0/0:71 0/0:47
    HE670865 368611 . C A 434.85 . AC=13;AF=0.813;AN=16;BaseQRankSum=0.412;DP=420;Dels=0.17;FS=0.883;HaplotypeScore=330.1362;MLEAC=14;MLEAF=0.875;MQ=42.34;MQ0=12;MQRankSum=0.904;QD=1.13;ReadPosRankSum=-2.718;SOR=0.631 GT:AD:DP:GQ:PL 0/0:17,11:28:3:0,3,23 1/1:27,10:38:3:32,3,0 1/1:23,23:46:6:64,6,0 1/1:12,16:29:3:32,3,0 1/1:22,30:53:9:73,9,0 1/1:18,35:53:6:56,6,0 1/1:26,34:60:18:188,18,0 0/1:23,17:40:16:16,0,61
    HE670865 368612 . A . 6.62 LowQual AN=16;DP=409;MQ=42.14;MQ0=12 GT:DP 0/0:27 0/0:35 0/0:36 0/0:26 0/0:49 0/0:47 0/0:55 0/0:31
    HE670865 368613 . G A 324.61 . AC=11;AF=0.917;AN=12;BaseQRankSum=3.644;DP=404;Dels=0.46;FS=1.094;HaplotypeScore=347.5139;MLEAC=11;MLEAF=0.917;MQ=41.95;MQ0=12;MQRankSum=1.325;QD=1.02;ReadPosRankSum=1.808;SOR=0.476 GT:AD:DP:GQ:PL 1/1:4,16:21:3:22,3,0 ./. 1/1:0,25:26:3:23,3,0 1/1:1,15:21:3:26,3,0 0/1:4,28:36:10:62,0,10 1/1:1,35:36:6:52,6,0 1/1:2,35:39:15:156,15,0 ./.
    HE670865 368614 . A . 17.69 LowQual AN=12;DP=416;MQ=41.75;MQ0=12 GT:DP ./. 0/0:25 0/0:36 ./. 0/0:47 0/0:41 0/0:50 0/0:25
    HE670865 368615 . A . 16.80 LowQual AN=12;DP=411;MQ=41.66;MQ0=11 GT:DP ./. 0/0:25 0/0:35 ./. 0/0:45 0/0:40 0/0:49 0/0:22
    HE670865 368616 . A . 18.07 LowQual AN=10;DP=406;MQ=41.71;MQ0=10 GT:DP ./. 0/0:24 ./. ./. 0/0:40 0/0:40 0/0:47 0/0:22
    HE670865 368617 . T C 92.95 . AC=4;AF=0.667;AN=6;BaseQRankSum=0.785;DP=385;Dels=0.38;FS=7.489;HaplotypeScore=362.4384;MLEAC=5;MLEAF=0.833;MQ=41.72;MQ0=9;MQRankSum=1.791;QD=0.83;ReadPosRankSum=4.858;SOR=0.297 GT:AD:DP:GQ:PL ./. 0/0:15,1:20:3:0,3,32 ./. ./. ./. 1/1:5,31:36:9:85,9,0 1/1:14,28:42:3:26,3,0 ./.
    HE670865 368618 . C . 18.70 LowQual AN=6;DP=370;MQ=41.65;MQ0=8 GT:DP ./. 0/0:10 ./. ./. ./. 0/0:34 0/0:37 ./.
    HE670865 368619 . G A 34.15 . AC=4;AF=0.500;AN=8;BaseQRankSum=-3.350;DP=365;Dels=0.54;FS=0.899;HaplotypeScore=399.7009;MLEAC=4;MLEAF=0.500;MQ=41.71;MQ0=8;MQRankSum=-1.167;QD=0.34;ReadPosRankSum=-4.849;SOR=1.069 GT:AD:DP:GQ:PL ./. 1/1:0,4:8:3:28,3,0 ./. ./. 0/0:17,3:22:3:0,3,23 0/0:31,3:34:3:0,3,22 1/1:28,8:36:3:25,3,0 ./.
    HE670865 368620 . T . 19.41 LowQual AN=2;DP=358;MQ=41.72;MQ0=8 GT:DP ./. 0/0:8 ./. ./. ./. ./. ./. ./.
    HE670865 368621 . T . . . DP=358;MQ=41.63;MQ0=8 GT ./. ./. ./. ./. ./. ./. ./. ./.
    HE670865 368622 . T . 15.52 LowQual AN=6;DP=358;MQ=41.53;MQ0=8 GT:DP ./. ./. ./. ./. ./. 0/0:11 0/0:19 0/0:17
    HE670865 368623 . T . 15.52 LowQual AN=6;DP=365;MQ=41.46;MQ0=7 GT:DP ./. 0/0:24 ./. 0/0:16 ./. 0/0:20 ./. ./.
    HE670865 368624 . A . 14.81 LowQual AN=8;DP=366;MQ=41.48;MQ0=7 GT:DP 0/0:22 0/0:26 ./. 0/0:23 ./. 0/0:20 ./. ./.
    HE670865 368625 . T . 15.53 LowQual AN=6;DP=353;MQ=41.39;MQ0=7 GT:DP 0/0:20 0/0:26 ./. ./. ./. 0/0:18 ./. ./.
    HE670865 368626 . T . . . DP=348;MQ=41.31;MQ0=7 GT ./. ./. ./. ./. ./. ./. ./. ./.
    HE670865 368627 . T A 11.69 LowQual AC=2;AF=1.00;AN=2;BaseQRankSum=-1.284;DP=348;Dels=0.55;FS=0.000;HaplotypeScore=403.2395;MLEAC=2;MLEAF=1.00;MQ=41.29;MQ0=7;MQRankSum=-0.971;QD=0.18;ReadPosRankSum=-0.344;SOR=1.126 GT:AD:DP:GQ:PL ./. ./. ./. ./. 1/1:30,4:34:3:26,3,0 ./. ./. ./.
    HE670865 368628 . T . 14.80 LowQual AN=8;DP=357;MQ=41.07;MQ0=8 GT:DP ./. 0/0:22 ./. 0/0:10 0/0:33 0/0:15 ./. ./.
    HE670865 368629 . A . 15.78 LowQual AN=10;DP=357;MQ=40.99;MQ0=8 GT:DP 0/0:14 0/0:22 ./. ./. 0/0:34 0/0:15 0/0:21 ./.
    HE670865 368630 . T . 19.39 LowQual AN=10;DP=358;MQ=40.89;MQ0=8 GT:DP 0/0:15 0/0:23 ./. 0/0:12 ./. 0/0:16 0/0:21 ./.
    HE670865 368631 . C . 18.51 LowQual AN=14;DP=359;MQ=40.70;MQ0=8 GT:DP 0/0:15 0/0:23 0/0:14 ./. 0/0:30 0/0:16 0/0:20 0/0:33
    HE670865 368632 . G . 17.87 LowQual AN=16;DP=361;MQ=40.63;MQ0=9 GT:DP 0/0:15 0/0:23 0/0:15 0/0:12 0/0:29 0/0:16 0/0:21 0/0:33
    HE670865 368633 . T . 19.08 LowQual AN=14;DP=364;MQ=40.51;MQ0=10 GT:DP 0/0:17 0/0:24 0/0:16 ./. 0/0:32 0/0:17 0/0:22 0/0:34
    HE670865 368634 . A . 17.17 LowQual AN=14;DP=369;MQ=40.67;MQ0=10 GT:DP 0/0:18 0/0:24 0/0:18 ./. 0/0:36 0/0:19 0/0:22 0/0:34
    HE670865 368635 . T . 17.39 LowQual AN=14;DP=364;MQ=40.66;MQ0=10 GT:DP 0/0:18 0/0:24 0/0:22 ./. 0/0:37 0/0:20 0/0:21 0/0:32
    HE670865 368636 . C . 16.20 LowQual AN=14;DP=362;MQ=40.70;MQ0=10 GT:DP 0/0:16 ./. 0/0:15 0/0:14 0/0:36 0/0:20 0/0:22 0/0:29
    HE670865 368637 . A . 0 LowQual AN=8;DP=356;MQ=40.91;MQ0=11 GT:DP 0/0:26 ./. 0/0:49 0/0:31 0/0:55 ./. ./. ./.
    HE670865 368638 . T G 23.71 LowQual AC=3;AF=0.250;AN=12;BaseQRankSum=-6.144;DP=355;Dels=0.14;FS=8.612;HaplotypeScore=349.9679;MLEAC=3;MLEAF=0.250;MQ=40.80;MQ0=12;MQRankSum=-2.191;QD=0.24;ReadPosRankSum=1.265;SOR=1.539 GT:AD:DP:GQ:PL ./. 0/0:22,3:25:6:0,6,44 0/0:46,3:49:6:0,6,44 0/0:29,1:30:6:0,6,49 0/0:46,8:54:3:0,3,22 1/1:40,2:42:3:29,3,0 0/1:47,5:52:17:18,0,17 ./.

    and for EMIT_ALL_CONFIDENT_SITES:

    HE670865 368589 . A . 46.89 . AN=16;DP=441;MQ=43.26;MQ0=13 GT:DP 0/0:35 0/0:45 0/0:58 0/0:38 0/0:75 0/0:64 0/0:64 0/0:47
    HE670865 368590 . A . 50.15 . AN=16;DP=447;MQ=43.41;MQ0=13 GT:DP 0/0:35 0/0:45 0/0:60 0/0:39 0/0:76 0/0:66 0/0:65 0/0:47
    HE670865 368591 . A . 39.12 . AN=16;DP=448;MQ=43.32;MQ0=12 GT:DP 0/0:36 0/0:45 0/0:61 0/0:39 0/0:77 0/0:66 0/0:66 0/0:46
    HE670865 368609 . A . 30.14 . AN=14;DP=434;MQ=42.67;MQ0=12 GT:DP ./. 0/0:44 0/0:53 0/0:39 0/0:67 0/0:61 0/0:74 0/0:47
    HE670865 368648 . G . 51.60 . AN=16;DP=366;MQ=40.87;MQ0=12 GT:DP 0/0:34 0/0:40 0/0:57 0/0:37 0/0:59 0/0:44 0/0:53 0/0:42
    HE670865 368649 . T . 75.72 . AN=16;DP=368;MQ=40.89;MQ0=12 GT:DP 0/0:34 0/0:41 0/0:57 0/0:38 0/0:59 0/0:44 0/0:53 0/0:42
    HE670865 368650 . T . 72.74 . AN=16;DP=369;MQ=40.87;MQ0=12 GT:DP 0/0:34 0/0:41 0/0:57 0/0:38 0/0:59 0/0:44 0/0:53 0/0:43
    HE670865 368651 . T . 87.71 . AN=16;DP=369;MQ=40.88;MQ0=12 GT:DP 0/0:33 0/0:42 0/0:58 0/0:37 0/0:59 0/0:44 0/0:53 0/0:43
    HE670865 368652 . T . 93.72 . AN=16;DP=373;MQ=40.89;MQ0=12 GT:DP 0/0:32 0/0:43 0/0:59 0/0:37 0/0:59 0/0:45 0/0:54 0/0:43
    HE670865 368653 . A G 10069.28 . AC=16;AF=1.00;AN=16;BaseQRankSum=1.361;DP=370;Dels=0.00;FS=0.000;HaplotypeScore=82.1575;MLEAC=16;MLEAF=1.00;MQ=40.92;MQ0=11;MQRankSum=0.132;QD=27.21;ReadPosRankSum=-0.590;SOR=0.652 GT:AD:DP:GQ:PL 1/1:1,31:32:69:797,69,0 1/1:0,42:43:99:1285,102,0 1/1:0,57:57:99:1753,150,0 1/1:0,36:37:84:1051,84,0 1/1:0,58:59:99:1485,126,0 1/1:0,44:44:96:1213,96,0 1/1:0,53:53:99:1479,123,0 1/1:0,43:44:90:1019,90,0

    Again it's still removing non LowQual sites, and I've no idea why.

  • tommycarstensentommycarstensen United KingdomMember
    edited August 2015

    @Jakemorris89 Can you post your command line in addition to what you have already posted? The reason I am asking is because this appears in your log file:

    standard_min_confidence_threshold_for_calling=30.0 standard_min_confidence_threshold_for_emitting=30.0
    

    And you said you had set those values equal to zero. So something doesn't add up. I think we can solve it, if you post the command line you ran.

  • Hi @tommycarstensen ,

    Sorry, I must not have been clear. In the cases above I did have those values set to the default 30, but the output is the same when I've set them to 0 as well. Unfortunately, that doesn't seem to be what makes the difference.

    I tried to explain this above:

    'Note: these are the runs with default call and emit thresholds, but I tried it with them set to 0 too and same thing, besides I want this to work with thresholds set at default to remove LowQual sites, (otherwise I could just use EMIT_ALL sites) I just don't know why EMIT_ALL_CONFIDENT_SITES is removing what seem like confident/'good' sites too.'

    I only didn't paste the command and output with them set to 0, as it was a fairly long post already, and trying what you suggested I can rule those out myself as being the issue. Hope that's clear now, and thanks for your suggestions. Any others would be welcome, but those thresholds and the version don't seem to be the issue.:-)

  • tommycarstensentommycarstensen United KingdomMember

    @Jakemorris89 OK, but can you post your command line? If you indent with four blank spaces, then it will be formatted in a grey box.

  • Jakemorris89Jakemorris89 UKMember
    edited August 2015
    java -Xmx4g -jar /GenomeAnalysisTK-3.4-46/genomeAnalysisTK.jar -T UnifiedGenotyper -out_mode EMIT_ALL_CONFIDENT_SITES --standard_min_confidence_threshold_for_calling 0 --standard_min_confidence_threshold_for_emitting 0 -baq CALCULATE_AS_NECESSARY -hets 0.015 -R genome_ref.fasta -o test_confidentsites_NewV_setthesholds0.vcf -I test.realign.bam &> test_confidentsites_NewV_setthesholds0.GATKlog &
    

    Thanks

  • tommycarstensentommycarstensen United KingdomMember

    @Jakemorris89 I give temporarily up. I will look at it with fresh eyes tomorrow.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator

    @Jakemorris89
    Hi Jake,

    Can you tell me what kind of data you are working with and how you pre-processed it? Did you follow our Best Practices? Can you also post the IGV screenshots of the bam file that correspond to the VCF records you posted above?

    Thanks,
    Sheila

  • Hi @Sheila,

    The data is Illumina Hiseq 125bp paired end, with approx 30x coverage. There are 8 samples mapped to the reference genome in the .bam file.

    The pre-processing was (in detail) as follows:

    1) For each sample trim and remove poor quality ends of reads using cutadapt
    2) Map to reference using STAMPY aligner and --processpart argument.
    3) As this produces many .Sam files, use samtools -view to compress
    4) Samtools sort on each .Bam
    5) Merge --processpart Bams for each sample.
    6) Now sort merged Bams
    7) At this step I used samtools view to make smaller .bams for each sample that contained a single scaffold of interest. This was done for 5 scaffolds.*
    8) The scaffold bams for each sample were then merged using samtools.
    9) I then removed duplicates with MarkDuplicates.jar
    10) Then I merged all the mark duplicates bams
    11) I Indexed this merged bam
    12) Then I realigned around indels

    *I have also run UnifiedGenotyper with just the EMIT_ALL_CONFIDENT_SITES option on a .bam (I no longer have this .bam) where step 7 and 8 were missed out, so that it contained the whole genome. The VCF file output for the region I gave above was the same as that given from the .bam which contained just 5 scaffolds. Note, I did not run this full genome .bam with EMIT_ALL_SITES.

    This 5 scaffold indel realigned .bam I then ran with GATK UnifiedGenotyper using both the EMIT_ALL_SITES and EMIT_ALL_CONFIDENT_SITES option and the latter seems to be omitting sites that are not marked as LowQual and that seem to be 'good quality'. You can see that in the outputs the two options gave (see above posts).

    The IGV image for the region from the VCF file is attached.

    Thanks,
    Jake

    image

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    We're no longer supporting UnifiedGenotyper. GATK Best Practices is now HaplotypeCaller only.

Sign In or Register to comment.