How can I get UnifiedGenotyper to work on triploids

I am using the Gatk v3.1 UnifiedGenotyper for variant calling of some GBS type data in the plant genus Boechera. I have diploids as well as triploids, where I have the ploidy from a prior microsat run. The diploids work out well (I get some 150k variable loci), but I can't seem to get variants for the triploids.

Here's my initial call for the triploids:

Program Args: -T UnifiedGenotyper -R /home/A01768064/assembly/Bstricta_278_v1.fa -I mytripbam.list -o boechera_triploids.vcf -nt 22 -glm SNP -hets 0.001 -mbq 20 -ploidy 3 -pnrm EXACT_GENERAL_PLOIDY -stand_call_conf 50 -maxAltAlleles 2

I was able to call the triploids with the diploid option, but - obviously - this is not what I want.
I additionally ran the following relaxed stringency options.

hets 0.01
mbq 10
hets 0.1 mbq 10

None of them return any variants.

The read group for one example file:

@RG ID:CR1308 LB:CR1308 SM:CR1308 PL:ILLUMINA

(with unique values for each sample; note that in the output vcf, the samples are appearing in the header line)

Additionally, I attached the slurm log file for the last mentioned run.

I am at a complete loss, what am I missing here?

Best Answer

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi Martin,

    Are you able to get variant calls if you get rid of all the non-required parameters except ploidy? Also, have you tried the latest version, 3.5? You should really try doing this with HaplotypeCaller in the new version. UG is deprecated at this point.

  • mschillingmschilling Logan UTMember

    ok, tried the same wth v3.5, with the same results. Will give the HaplotypeCaller a try.

  • mschillingmschilling Logan UTMember

    ok, finally back...

    I tried to run the HaplotypeCaller (v3.5), without success.
    I used the same input data, with the following command

    java -jar -Xmx476g -Djava.io.tmpdir=/pscratch/A01768064/gbs15 /home/A01768064/gatk/GenomeAnalysisTK.jar -T HaplotypeCaller -R /home/A01768064/assembly/Bstricta_278_v1.fa -I mydipbam.list -o boechera_HC_dip_gatk35.vcf -ploidy 2

    I tried a second run with the triploid individuals and -ploidy 3, without success.

    With the UnifiedGenotyper, at least my diploids produced SNVs. Now, none of them seem to work.
    Again, as in the above case, I do get an output vcf with comment lines up to the '#CHROM..." line and that's it. No error or warning message I can find in stdout...

    any help would be greatly appreciated.
    thanks,
    martin

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Can you post some IGV screenshots of the sites called in your triploids as diploids with UG? Do they look reasonable? And can you run HC in debug mode on those same sites? See this article for instructions on how to debug HC calls. See also this article for general recommendations dealing with failure to call expected sites.

  • mschillingmschilling Logan UTMember

    ok, here are two of the lines from the triploids when called as diploids. igv screenshots are attached

    Scaffold116 277 . T A 30.60 LowQual AC=3;AF=0.056;AN=54;DP=1354;Dels=0.00;FS=0.000;HaplotypeScore=0.0459;InbreedingCoeff=-0.1422;MLEAC=3;MLEAF=0.056;MQ=37.00;MQ0=0;QD=0.20 GT:AD:DP:GQ:PL 0/0:40,0:40:9:0,9,120 0/0:24,0:24:3:0,3,40 0/0:24,0:24:12:0,12,160 0/0:83,0:83:15:0,15,200 ./. 0/0:107,0:107:18:0,18,240 0/1:71,1:72:28:28,0,108 ./. 0/0:8,0:8:6:0,6,80 0/0:94,0:94:18:0,18,240 0/0:135,0:135:39:0,39,520 ./. ./. ./. ./. ./. ./.:.:13 ./. 0/0:33,0:33:6:0,6,80 0/1:62,1:63:28:28,0,108 0/0:59,0:59:18:0,18,240 0/0:38,0:38:6:0,6,80 ./. 0/0:67,0:67:9:0,9,120 0/0:15,0:15:3:0,3,40 ./. 0/1:18,1:19:34:34,0,34 ./. ./. ./. 0/0:12,0:12:3:0,3,40 0/0:1,0:1:3:0,3,40 0/0:28,0:28:6:0,6,80 0/0:42,0:42:9:0,9,120 0/0:23,0:23:15:0,15,200 0/0:16,0:16:9:0,9,120 ./. ./. 0/0:127,0:127:27:0,27,360 ./. ./.:.:22 ./. 0/0:16,0:16:3:0,3,40 0/0:45,0:45:9:0,9,120 0/0:18,0:18:6:0,6,80 0/0:2,0:2:3:0,3,40 ./.


    image

    and

    Scaffold116 280 . T A 61.50 . AC=6;AF=0.103;AN=58;DP=1354;Dels=0.00;FS=0.000;HaplotypeScore=0.0459;InbreedingCoeff=0.1138;MLEAC=5;MLEAF=0.086;MQ=37.00;MQ0=0;QD=0.37 GT:AD:DP:GQ:PL 0/0:40,0:40:9:0,9,120 0/0:24,0:24:6:0,6,80 0/0:24,0:24:6:0,6,80 ./. ./. 0/0:107,0:107:12:0,12,160 0/0:72,0:72:9:0,9,120 ./. 0/0:8,0:8:3:0,3,40 0/0:94,0:94:21:0,21,280 0/1:134,1:135:13:13,0,293 0/0:7,0:7:3:0,3,40 0/0:8,0:8:3:0,3,40 0/0:10,0:10:3:0,3,40 ./. ./. 0/0:13,0:13:6:0,6,80 ./. 0/0:33,0:33:6:0,6,80 0/0:63,0:63:6:0,6,80 0/0:59,0:59:6:0,6,80 0/0:38,0:38:3:0,3,40 ./. 0/0:67,0:67:9:0,9,120 1/1:14,1:15:3:40,3,0 ./. 0/0:19,0:19:9:0,9,120 ./. ./. ./. ./. ./. 0/0:28,0:28:9:0,9,120 0/0:42,0:42:9:0,9,120 0/0:23,0:23:6:0,6,80 0/1:15,1:16:34:34,0,34 ./. ./. 0/0:127,0:127:24:0,24,320 ./. 0/0:22,0:22:6:0,6,80 ./. ./. 0/0:45,0:45:15:0,15,200 0/0:18,0:18:3:0,3,40 1/1:1,1:2:3:40,3,0 0/0:15,0:15:3:0,3,40

    image

    I noticed that igv gives base qualities for the two variants as -17 and -4, respectively.
    Negative phred scaled base qualities do not make sense to me, how come I get those?

    thanks,
    m

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @mschilling
    Hi,

    So, you used Illumina to sequence your samples? Illumina should not be giving those negative base quality scores. Can you tell us how exactly you processed your files after you got them from the sequencer? Did you follow the Best Practices?

    -Sheila

  • mschillingmschilling Logan UTMember
    edited December 2015

    Hey Sheila,
    considering that the Best Practices contain variant calling and call-set refinement, I have not been able to finish them (variant calling with the ploidy 3 option does not return any variants, whereas the diploid data work...)

    After initial parsing of barcodes (GBS data) and splitting individuals into separate files, the reads were aligned to the JGI Boechera stricta assembly using bwa samse and then converted to bam using samtools.

    thanks,
    m

    Issue · Github
    by Sheila

    Issue Number
    435
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    vdauwera
  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @mschilling
    Hi m,

    Unfortunately, I'm not sure we can help you out much here. We don't support GBS data, but can you make sure your bam files are valid? You can use Picard's ValidateSamFile.

    As for the VCF not containing any variant calls, it looks like most of the sites don't have enough evidence for a variant call. The few samples that do get called variant show maybe 1 read with a variant allele out of all the reads. That is really not enough evidence to have confidence in the variant call.

    -S

Sign In or Register to comment.