Test-drive the GATK tools and Best Practices pipelines on Terra


Check out this blog post to learn how you can get started with GATK and try out the pipelines in preconfigured workspaces (with a user-friendly interface!) without having to install anything.

Optimal batch sizes for GenomeSTRiP

jfarrelljfarrell Member ✭✭

How large a batch of bam files would be recommended for 40x coverage with the latest version of GenomeSTRiP? Presently we would like to run this on about 65 bam files initially. For discovery, is the recommended best practice to run on all of them at once or split them into two smaller groups to keep the coverage below 2000x for a run?

Best Answers

Answers

  • jfarrelljfarrell Member ✭✭

    What is the best set of 1000 genomes SV sites vcf to augment the discovery sites for genotypes?

    I found the Phase 1 (ALL.wgs.BI_genome_strip_hq.20101123.dels.low_coverage.genotypes.vcf.gz). Is there a more recent vcf of Phase 3 SVs on the 1000 Genomes Phase 3 2500 samples available?

  • bhandsakerbhandsaker Member, Broadie, Moderator admin

    The objective is to capture any common deletions that might have been missed by Genome STRiP discovery. So I think Phase 1 is fine, but I would recommend using all of the deletion sites from the integrated call set, which contains calls from 5 or so different methods.

    You can get the site list from the DCC here:

    ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase1/analysis_results/integrated_call_sets/ALL.wgs.integrated_phase1_v3.20101123.snps_indels_sv.sites.vcf.gz

    but you will have to grep out the large deletions (VT=SV in the INFO field).

    Or you can get just the SV sites from our website (should be identical to the records in the above VCF):

    ftp://ftp.broadinstitute.org/pub/svtoolkit/misc/1kg/1kg_phase1_release_v3/svs/ALL.genome.phase1_release_v3.20101123.svs.sites.vcf.gz

  • jfarrelljfarrell Member ✭✭

    Thanks Bob!

    I ran the genotyping step on the 37 bam files with the .ALL.genome.phase1_release_v3.20101123.svs.sites.vcf as the input sites. There were a couple errors reported. As background, GenomeSTRiP genotyping step had successfully completed with the sites previously generated from the GenomeSTRiP Discovery.

    There seems to be 2 VCF SV records triggering the error ( MERGED_DEL_2_82399 and MERGED_DEL_2_24813)

    This first error is logged in SVGenotyper-10.out

    **Program Args**: -T SVGenotyperWalker -R /archive/not-replicated/project/genpro/archive/adsp/ref/Homo_sapiens_assembly19.fasta -O /auto/nfs-archive/ifs/noreplica/project/genpro/archive/adsp/GenomeSTRiP/svtoolkit/adsp/run.Pilot37_kg/P0116.genotypes.vcf -disableGATKTraversal true -md run.Pilot37/metadata -configFile conf/genstrip_parameters.txt -P pairs.uniquifyReadNames:true -runDirectory run.Pilot37_kg -genderMapFile adsp_gender.map -genomeMaskFile /usr3/bustaff/farrell/adsp/GenomeSTRiP/svtoolkit/genomeMaskFile/Homo_sapiens_assembly19.mask.101.fasta -vcf /archive/not-replicated/project/genpro/archive/adsp/GenomeSTRiP/svtoolkit/1000genomes/ALL.genome.phase1_release_v3.20101123.svs.sites.vcf -partitionName P0116 -partition records:11501-11600 -L 1:1-1
    
    
    
        SVGenotyper-10.out:#ERR: SAM record: 482059550  147     16      13944304        46      66S32M104D3M    =       13943981        -462    ACACTTGAACCTGGGGGGCGGGGTTTGCAGTCAGCCAAGATTGTGCCCCTGCCCTCCAGCCTGGGGGACAGAGTAAAACTGCGTCAAAAAAAAAAAAAAAA        ##########################################################################[email protected]=3,,'=5'DIIDDDDDDDDDD???   RG:Z:C27PG.1    NH:i:1  NM:i:104
        SVGenotyper-10.out:#ERR: junction: MERGED_DEL_2_82399_1 401 1-200 16:13944137-13944336 202-401 16:13944440-13944639 CAACTCTACTAAAAATAAAAAAATAAAAATAAAAATAAATTAGCTGGGTGTGATGCCACATGCCTGTAATCCCAGCTACTCAGGAGGCTGAGGCAGGAAAATCACTTGAACTTGGGAGGCGGAGTTTGCAGTCAGCCAAGATTGTGCCACTGCACTCCAGCCTGGGTGACAGAGTAAAACTGCGTCAAAAAAAAAAAAAAAAAAAATGCTTTTTTCAACTGCCTTTCCTATTTAGTCACAGTAGATATAGAGATGGAATCTGAGTTTTCCCAGACAAAAAGTGGCAAAACTGAGGTCATTTCAGGCACTCTACTTGGACATTTTCAAGAACTAGAGAATCCCTGATGTGCTGTAATCAATTCAGAAGGTAGCATTAAGTCCCACAAAAAGTCATCATCAAG
        SVGenotyper-10.out:#ERR: MERGED_DEL_2_82399_1 J CAACTCTACTAAAAATAAAAAAATAAAAATAAAAATAAATTAGCTGGGTGTGATGCCACATGCCTGTAATCCCAGCTACTCAGGAGGCTGAGGCAGGAAAATCACTTGAACTTGGGAGGCGGAGTTTGCAGTCAGCCAAGATTGTGCCACTGCACTCCAGCCTGGGTGACAGAGTAAAACTGCGTCAAAAAAAAAAAAAAAAAAAATGCTTTTTTCAACTGCCTTTCCTATTTAGTCACAGTAGATATAGAGATGGAATCTGAGTTTTCCCAGACAAAAAGTGGCAAAACTGAGGTCATTTCAGGCACTCTACTTGGACATTTTCAAGAACTAGAGAATCCCTGATGTGCTGTAATCAATTCAGAAGGTAGCATTAAGTCCCACAAAAAGTCATCATCAAG
        SVGenotyper-10.out:#ERR: MERGED_DEL_2_82399_1   ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
        SVGenotyper-10.out:#ERR: MERGED_DEL_2_82399_1 L CAACTCTACTAAAAATAAAAAAATAAAAATAAAAATAAATTAGCTGGGTGTGATGCCACATGCCTGTAATCCCAGCTACTCAGGAGGCTGAGGCAGGAAAATCACTTGAACTTGGGAGGCGGAGTTTGCAGTCAGCCAAGATTGTGCCACTGCACTCCAGCCTGGGTGACAGAGTAAAACTGCGTCAAAAAAAAAAAAAAGTTTTAAAAATAAAGTTTATCTTTTGTTCTAAAACAAAGTTTTAAAAAAACTACTACAAATGAAACATAGGACTAGCTGACCACATTTCATTTGCACTTCATCAAAAATGCTTTTTTCAACTGCCTTTCCTATTTAGTCACAGTAGATATAGAGATGGAATCTGAGTTTTCCCAGACAAAAAGTGGCAAAACTGAGGTCAT
        SVGenotyper-10.out:#ERR: MERGED_DEL_2_82399_1 J CAACTCTACTAAAAATAAAAAAATAAAAATAAAAATAAATTAGCTGGGTGTGATGCCACATGCCTGTAATCCCAGCTACTCAGGAGGCTGAGGCAGGAAAATCACTTGAACTTGGGAGGCGGAGTTTGCAGTCAGCCAAGATTGTGCCACTGCACTCCAGCCTGGGTGACAGAGTAAAACTGCGTCAAAAAAAAAAAAAAAAAAAATGCTTTTTTCAACTGCCTTTCCTATTTAGTCACAGTAGATATAGAGATGGAATCTGAGTTTTCCCAGACAAAAAGTGGCAAAACTGAGGTCATTTCAGGCACTCTACTTGGACATTTTCAAGAACTAGAGAATCCCTGATGTGCTGTAATCAATTCAGAAGGTAGCATTAAGTCCCACAAAAAGTCATCATCAAG
        SVGenotyper-10.out:#ERR: MERGED_DEL_2_82399_1                                                                                                                                                                                                            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
        SVGenotyper-10.out:#ERR: MERGED_DEL_2_82399_1 R CACTTGAACTTGGGAGGCGGAGTTTGCAGTCAGCCAAGATTGTGCCACTGCACTCCAGCCTGGGTGACAGAGTAAAACTGCGTCAAAAAAAAAAAAAAGTTTTAAAAATAAAGTTTATCTTTTGTTCTAAAACAAAGTTTTAAAAAAACTACTACAAATGAAACATAGGACTAGCTGACCACATTTCATTTGCACTTCATCAAAAATGCTTTTTTCAACTGCCTTTCCTATTTAGTCACAGTAGATATAGAGATGGAATCTGAGTTTTCCCAGACAAAAAGTGGCAAAACTGAGGTCATTTCAGGCACTCTACTTGGACATTTTCAAGAACTAGAGAATCCCTGATGTGCTGTAATCAATTCAGAAGGTAGCATTAAGTCCCACAAAAAGTCATCATCAAG
        SVGenotyper-10.out:#ERR:            13944337     13944439
        SVGenotyper-10.out:#ERR: AAAAAAAA | GTTTTAAA ... ACTTCATC | AAAAATGC
        SVGenotyper-10.out:#ERR: AAAAAAAA | A                   A | AAAAATGC
        SVGenotyper-10.out:#ERR:      200                           202
        SVGenotyper-10.out:#ERR: junctionEnd: R
        SVGenotyper-10.out:#ERR: refExtentStart: 13944238
        SVGenotyper-10.out:#ERR: refExtentEnd: 13944442
        SVGenotyper-10.out:#ERR: junctionLength: 401
        SVGenotyper-10.out:#ERR: junctionLeftFlankLength: 200
        SVGenotyper-10.out:#ERR: junctionRightFlankLength: 200
        SVGenotyper-10.out:##### ERROR ------------------------------------------------------------------------------------------
        SVGenotyper-10.out:##### ERROR stack trace
        SVGenotyper-10.out:##### ERROR ------------------------------------------------------------------------------------------
        SVGenotyper-10.out:##### ERROR A GATK RUNTIME ERROR has occurred (version 2.5-2-g5671483):
        SVGenotyper-10.out:##### ERROR
        SVGenotyper-10.out:##### ERROR Please check the documentation guide to see if this is a known problem
        SVGenotyper-10.out:##### ERROR If not, please post the error, with stack trace, to the GATK forum
        SVGenotyper-10.out:##### ERROR Visit our website and forum for extensive documentation and answers to
        SVGenotyper-10.out:##### ERROR commonly asked questions http://www.broadinstitute.org/gatk
        SVGenotyper-10.out:##### ERROR
    

    The 2nd error was logged in SVGenotyper-60.out

    Program Args: -T SVGenotyperWalker -R /archive/not-replicated/project/genpro/archive/adsp/ref/Homo_sapiens_assembly19.fasta -O /auto/nfs-archive/ifs/noreplica/project/genpro/archive/adsp/GenomeSTRiP/svtoolkit/adsp/run.Pilot37_kg/P0035.genotypes.vcf -disableGATKTraversal true -md run.Pilot37/metadata -configFile conf/genstrip_parameters.txt -P pairs.uniquifyReadNames:true -runDirectory run.Pilot37_kg -genderMapFile adsp_gender.map -genomeMaskFile /usr3/bustaff/farrell/adsp/GenomeSTRiP/svtoolkit/genomeMaskFile/Homo_sapiens_assembly19.mask.101.fasta -vcf /archive/not-replicated/project/genpro/archive/adsp/GenomeSTRiP/svtoolkit/1000genomes/ALL.genome.phase1_release_v3.20101123.svs.sites.vcf -partitionName P0035 -partition records:3401-3500 -L 1:1-1

    INFO  08:51:23,967 ProgressMeter -        Starting        0.00e+00    5.2 m      512.6 w    100.0%         5.2 m     0.0 s
    Error: Exception processing cnp: 400
    CNP: MERGED_DEL_2_24813 4:59985300-59985588
    INFO  08:51:25,323 GATKRunReport - Uploaded run statistics report to AWS S3
    ##### ERROR ------------------------------------------------------------------------------------------
    ##### ERROR stack trace
    java.lang.ArrayIndexOutOfBoundsException: 400
            at org.broadinstitute.sv.genotyping.SplitReadMapper.getMismatchQualitySum(SplitReadMapper.java:700)
            at org.broadinstitute.sv.genotyping.SplitReadMapper.annotateSplitReadMapped(SplitReadMapper.java:435)
            at org.broadinstitute.sv.genotyping.SplitReadMapper.annotateSplitRead(SplitReadMapper.java:184)
            at org.broadinstitute.sv.genotyping.SplitReadMapper.getSplitReadsInternal(SplitReadMapper.java:152)
            at org.broadinstitute.sv.genotyping.SplitReadMapper.getSplitReads(SplitReadMapper.java:134)
            at org.broadinstitute.sv.genotyping.SplitReadMapper.getSplitReads(SplitReadMapper.java:112)
            at org.broadinstitute.sv.genotyping.GenotypingSplitReadModule.getSplitReads(GenotypingSplitReadModule.java:506)
            at org.broadinstitute.sv.genotyping.GenotypingSplitReadModule.getSplitReadGenotypingReads(GenotypingSplitReadModule.java:371)
            at org.broadinstitute.sv.genotyping.GenotypingSplitReadModule.genotypeSample(GenotypingSplitReadModule.java:119)
            at org.broadinstitute.sv.genotyping.GenotypingSplitReadModule.genotypeCnp(GenotypingSplitReadModule.java:85)
            at org.broadinstitute.sv.genotyping.GenotypingSplitReadModule.genotypeCnp(GenotypingSplitReadModule.java:39)
            at org.broadinstitute.sv.genotyping.GenotypingAlgorithm.genotypeCnpInternal(GenotypingAlgorithm.java:143)
            at org.broadinstitute.sv.genotyping.GenotypingAlgorithm.genotypeCnp(GenotypingAlgorithm.java:108)
            at org.broadinstitute.sv.genotyping.SVGenotyperWalker.processVCFFile(SVGenotyperWalker.java:265)
            at org.broadinstitute.sv.genotyping.SVGenotyperWalker.map(SVGenotyperWalker.java:211)
            at org.broadinstitute.sv.genotyping.SVGenotyperWalker.map(SVGenotyperWalker.java:56)
            at org.broadinstitute.sting.gatk.traversals.TraverseLociNano.traverse(TraverseLociNano.java:107)
            at org.broadinstitute.sting.gatk.traversals.TraverseLociNano.traverse(TraverseLociNano.java:48)
            at org.broadinstitute.sting.gatk.executive.LinearMicroScheduler.execute(LinearMicroScheduler.java:100)
            at org.broadinstitute.sting.gatk.GenomeAnalysisEngine.execute(GenomeAnalysisEngine.java:286)
            at org.broadinstitute.sting.gatk.CommandLineExecutable.execute(CommandLineExecutable.java:113)
            at org.broadinstitute.sv.main.SVCommandLine.execute(SVCommandLine.java:125)
            at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:245)
            at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:152)
            at org.broadinstitute.sv.main.SVCommandLine.main(SVCommandLine.java:79)
            at org.broadinstitute.sv.main.SVGenotyper.main(SVGenotyper.java:21)
    ##### ERROR ------------------------------------------------------------------------------------------
    ##### ERROR A GATK RUNTIME ERROR has occurred (version 2.5-2-g5671483):
    ##### ERROR
    ##### ERROR Please check the documentation guide to see if this is a known problem
    ##### ERROR If not, please post the error, with stack trace, to the GATK forum
    ##### ERROR Visit our website and forum for extensive documentation and answers to
    ##### ERROR commonly asked questions http://www.broadinstitute.org/gatk
    ##### ERROR
    ##### ERROR MESSAGE: 400
    ##### ERROR ------------------------------------------------------------------------------------------
    
  • bhandsakerbhandsaker Member, Broadie, Moderator admin

    Thanks for all of your great bug reports!

    For the first one, it would be nice to have the stack trace so I could see what the actual error is.

    For 1000G phase 1, we didn't use split-read genotyping with the large deletions (because we didn't have good data on the accuracy of the assembled breakpoints). You have a couple of choices:

    a) Turn off split-read genotyping globally: -P genotyping.modules:depth,pairs

    b) Turn off split-read genotyping just for these two sites. This is a little cumbersome. You could remove these two sites from your input file and run them separately (or just ignore them altogether), but then you would have to re-genotype all of the other sites (because the number of sites in the input file has changed, which changes the partitions). Alternatively, you could edit the input file to change the alt alleles to <DEL> (so these sites will not have exact breakpoints) and then Queue should just rerun the failed partitions.

  • jfarrelljfarrell Member ✭✭

    Thanks for the suggested workarounds.I will replace the ALT allele with and give that a try first.

    Below is the stack trace from the first one.

    ##### ERROR ------------------------------------------------------------------------------------------
    ##### ERROR stack trace
    java.lang.RuntimeException: Internal error: junctionOffset = -1
            at org.broadinstitute.sv.genotyping.SplitReadMapper.annotateSplitReadMapped(SplitReadMapper.java:432)
            at org.broadinstitute.sv.genotyping.SplitReadMapper.annotateSplitRead(SplitReadMapper.java:188)
            at org.broadinstitute.sv.genotyping.SplitReadMapper.getSplitReadsInternal(SplitReadMapper.java:152)
            at org.broadinstitute.sv.genotyping.SplitReadMapper.getSplitReads(SplitReadMapper.java:134)
            at org.broadinstitute.sv.genotyping.SplitReadMapper.getSplitReads(SplitReadMapper.java:112)
            at org.broadinstitute.sv.genotyping.GenotypingSplitReadModule.getSplitReads(GenotypingSplitReadModule.java:506)
            at org.broadinstitute.sv.genotyping.GenotypingSplitReadModule.getSplitReadGenotypingReads(GenotypingSplitReadModule.java:371)
            at org.broadinstitute.sv.genotyping.GenotypingSplitReadModule.genotypeSample(GenotypingSplitReadModule.java:119)
            at org.broadinstitute.sv.genotyping.GenotypingSplitReadModule.genotypeCnp(GenotypingSplitReadModule.java:85)
            at org.broadinstitute.sv.genotyping.GenotypingSplitReadModule.genotypeCnp(GenotypingSplitReadModule.java:39)
            at org.broadinstitute.sv.genotyping.GenotypingAlgorithm.genotypeCnpInternal(GenotypingAlgorithm.java:143)
            at org.broadinstitute.sv.genotyping.GenotypingAlgorithm.genotypeCnp(GenotypingAlgorithm.java:108)
            at org.broadinstitute.sv.genotyping.SVGenotyperWalker.processVCFFile(SVGenotyperWalker.java:265)
            at org.broadinstitute.sv.genotyping.SVGenotyperWalker.map(SVGenotyperWalker.java:211)
            at org.broadinstitute.sv.genotyping.SVGenotyperWalker.map(SVGenotyperWalker.java:56)
            at org.broadinstitute.sting.gatk.traversals.TraverseLociNano.traverse(TraverseLociNano.java:107)
            at org.broadinstitute.sting.gatk.traversals.TraverseLociNano.traverse(TraverseLociNano.java:48)
            at org.broadinstitute.sting.gatk.executive.LinearMicroScheduler.execute(LinearMicroScheduler.java:100)
            at org.broadinstitute.sting.gatk.GenomeAnalysisEngine.execute(GenomeAnalysisEngine.java:286)
            at org.broadinstitute.sting.gatk.CommandLineExecutable.execute(CommandLineExecutable.java:113)
            at org.broadinstitute.sv.main.SVCommandLine.execute(SVCommandLine.java:125)
            at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:245)
            at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:152)
            at org.broadinstitute.sv.main.SVCommandLine.main(SVCommandLine.java:79)
            at org.broadinstitute.sv.main.SVGenotyper.main(SVGenotyper.java:21)
    
  • jfarrelljfarrell Member ✭✭

    Bob,

    The replacement of the ALT allele with DEL in the VCF got rid of the error and the genotyping step finished up eventually. However, the genotyping of those 2 deletions each took over 20 hours which about 18 hours longer than any other site. After all that time, no genotypes were generated in the 37 samples for those 2 sites. All samples had LowQual calls with "." for a genotype.

    So the next time, I will go with your other suggestion and just remove them from the input vcf site file.

    Thanks,

    John

  • jfarrelljfarrell Member ✭✭

    For determining the 80% reciprocal overlap of GenomeSTRiP discovery sites with the 1000 genome sites, is there a GATK software tool or other tool available for finding the overlapping sites.

Sign In or Register to comment.