Asking for Guidance on Variant Calling Pipeline - Biological replicates

Follow

H.T

• November 14, 2023 14:07
• Edited

Hi GATK Team,

 

I am new to variant calling and could use some guidance on a pipeline I've been working on for SNP calling from Hi-C data in our lab.

 

We have five biological replicates, each resulting in a BAM file after alignment. Our objective is to pool information from these five files to enhance SNP calling and avoid losing variants due to low coverage or data sparsity in individual files.

 

After reviewing the GATK documentation, my initial approach involved performing joint variant calling (using `-ERC GVCF`) after processing the BAM files separately. This included steps like marking duplicates, applying BQSR, sorting, indexing, and ensuring consistent SM tags across the five files.

However, when I reached the HaplotypeCaller step, I encountered the error: "A USER ERROR has occurred: Argument emit-ref-confidence has a bad value: Can only be used in single-sample mode currently. Use the --sample-name argument to run on a single sample out of a multi-sample BAM file." I assumed this means that we need more than one sample for joint variant calling, not just multiple files per one sample. Is my understanding correct?

 

In response, I took a step back and merged all the BAM files together. However, I am uncertain whether running HaplotypeCaller in GVCF mode would bring any advantage to the pipeline in this case or if it's unnecessary.

 

For your reference, here are the versions of the tools I've used:

- The Genome Analysis Toolkit (GATK) v4.2.1.0

- HTSJDK Version: 2.24.1

- Picard Version: 2.25.4

 

I have added the final version of my pipeline.

https://crgcnag-my.sharepoint.com/:u:/g/personal/htoukabri_crg_es/EQGRCwe7G_FIrSJ1OXgtLbIBZ_I2bo_Vf4D_8E_u9Rwe4w?e=vUE3cz 

 

I would greatly appreciate your insights, any adjustment suggestions, or clarifications regarding the steps I've followed.

 

Thanks a lot! :)

0

• 

  Gökalp Çelik

  • November 14, 2023 16:56

  Hi H.T

  The message about GVCF mode being incompatible with multiple samples is correct. When HaplotypeCaller complains about multiple samples it means that each bam input has a different sample name or a pooled bam has multiple sample names present in it. If you really wish to distinguish you replicates using different sample names and generate a VCF with these different sample names you can run HaplotypeCaller just without the -ERC GVCF parameter and HaplotypeCaller will still generate a multisample VCF file. 

  You can check the header of your bam file using 

  samtools view -H

  If you observe multiple SM tags in the header then it is the reason for your troubles. 

  If your intention is to have a single sample name but a combined bam file just to increase depth then you need to use

  gatk AddOrReplaceReadGroups 

  tool to fix multiple sample names issue from your bam file. 

  I hope this helps. 

  0

  Comment actions Permalink

Please sign in to leave a comment.