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.

questions about downsampling

xiang_zuoxiang_zuo shanghaiMember

Hi, I am using HaplotypeCaller to call variants on targer region sequencing data, the mean coverage is about 600, I want to downsample my data to test the minumum mean coverage required for variant calling.So I use the parameter -dcov, the script is pasted below:
-T HaplotypeCaller -A StrandAlleleCountsBySample -L target_region.bed -R hg19.fasta -I sample.bam --genotyping_mode DISCOVERY -stand_emit_conf 10 -stand_call_conf 30 -dcov 100 -o sample.vcf

However, error message appears:
ERROR MESSAGE: Invalid command line: Cannot use -dcov or --downsample_to_coverage for ActiveRegionWalkers, use another downsampling argument

How can I solve the problem? Thanks a lot!

Best,

Answers

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    Hi @xiang_zuo,

    This thread says that we don't recommend downsampling with HaplotypeCaller. In addition to the solution in the thread, you can use PrintReads to downsample first then run HaplotypeCaller on the downsampled BAM.

  • xiang_zuoxiang_zuo shanghaiMember
    edited July 2017

    @shlee said:
    Hi @xiang_zuo,

    This thread says that we don't recommend downsampling with HaplotypeCaller. In addition to the solution in the thread, you can use PrintReads to downsample first then run HaplotypeCaller on the downsampled BAM.

    Hi,shlee.Thanks a lot!I used PrintReads to downsample first. The script is :smile:

    GenomeAnalysisTK.jar -T PrintReads \
        -R hg19.fasta \
        -I input.bam \
        -o output.bam \
        -dcov 100
    

    However, the size of the 2 bam files is the same, that is to say, no reads are substracted during the process. I am confused.

    Post edited by shlee on
  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    @xiang_zuo,

    That's odd given you have mean coverage ~600. Are you using the latest stable release (v3.7)?

  • xiang_zuoxiang_zuo shanghaiMember

    @shlee said:
    @xiang_zuo,

    That's odd given you have mean coverage ~600. Are you using the latest stable release (v3.7)?

    No, I am using v3.5. HC downsample to a mean coverae of 500 defaultly. If I want to use all the data instead of substracting it to 500X, how can I specify the parameters? Thanks!

  • shleeshlee CambridgeMember, Broadie ✭✭✭✭✭

    @xiang_zuo, can you make sure, via IGV visualization, that a region of ~600 coverage still remains ~600 in coverage after the PrintReads command?

    To disable the downsampling, I believe you include -dt NONE to the HaplotypeCaller command. The -dt stands for --downsampling_type. You can read more about this argument in the engine parameters.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @xiang_zuo
    Hi,

    You can also try setting --maxReadsInRegionPerSample to a lower number. The downsampling in GATK3 is confusing and we don't recommend messing with it. It should be easier to change the settings in GATK4.

    -Sheila

Sign In or Register to comment.