This is a test version of Biostars. For the public version, visit https://www.biostars.org.
GATK based variant calling yields excessive Indels

Dear all

I have a GATK based workflow for my .fq files:

  1. Trimming and adapter clipping by using Trimmomatic
  2. Read alignment against reference with Bowtie2 (all defaults)
  3. Conversion of .bam to .sam and sorting with Samtools
  4. Then I use GATK AddOrReplaceReadGroups:
java -jar picard.jar AddOrReplaceReadGroups \
I=$file \
O=$outfile \
SORT_ORDER=coordinate \
CREATE_INDEX=true VALIDATION_STRINGENCY=LENIENT \
RGID=unique_arbitrary_value \
RGLB= unique_arbitrary_value \
RGPL=ILLUMINA \
RGPU= unique_arbitrary_value \
RGSM= unique_arbitrary_value

Then GATK MarkDuplicates

java -jar $EBROOTPICARD/picard.jar MarkDuplicates \
I=$file \
O=$file.mkd_dup.bam \
M=$file.mkd_dup.txt

Then samtools index $file.mkd_dup.bam Then GATK HaplotypeCaller

gatk --java-options "-Xmx4g" HaplotypeCaller \
   -R REF.fna \
   -I $file.mkd_dup.bam \
   -O $file.g.vcf.gz \
   -ERC GVCF

What I have found was after filtering out loci with too much missing data etc my remaining combined gzipped vcf file had around only 600 SNP markers out of 20,000 remaining – with the rest being INDELS. This seems far too high a percentage of Indels – can anyone suggest improvements to my pipeline?

Incidentally, I merged my vcf.gz files with VCFtools: vcf-merge $gz | bgzip -c > mge.vcf.gz (gave me around 10million identified loci) Because merging with GATK: gatk --java-options "-Xmx4g" CombineGVCFs - created a weird file with ALT alleles called only as “<NON_REF>” (i.e. with no accompanying variant call e.g. “C,<NON_REF>”) and all calls for each sample as “./.”

indels calling gatk variant picard
• 1,823 views • link
updated by Pierre Lindenbaum Pierre Lindenbaum • written by evoclive • evoclive

1 answer

unless I don't understand what you're doing, you shouldn't call your VCF in GVCF mode (-ERC GVCF) if you don't plan to combine many GVCFs file with CombineGVCFs followed by a call to GenotypeGVCFs.

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Hi, well I did initially try to use CombineGVCFs with:

gatk --java-options "-Xmx4g" CombineGVCFs \
   -R REF.fna \
   --variant gatk/AORRG_1.g.vcf.gz \
   --variant gatk/AORRG_10.g.vcf.gz \

but every locus call looks like: CM014959.1 527473 . A <NON_REF> . . . GT:DP:GQ:MIN_DP:PL ./.:0:0:0:0,0,0 ./.:0:0:0:0,0,0 ./.:1:3:1:0,3,11 ./.:1:3:1:0,3,16 ./.:0:0:0:0,0,0 ./.:0:0:0:0,0,0 ./.:0:0:0:0,0,0 ./.:0:0:0:0,0,0 ./.:0:0:0:0,0,0 ./.:0:0:0:0,0,0

Which looks like nonsense to me. Do I need to use GenotypeGVCFs after that? Thus, I used vcftools.

• 0 views • link
written by evoclive • evoclive

Aaah! Very sorry - thanks for your time. Is there a one-document tutorial for this process outlining a full work-flow and what arguments to use under different circumstances (e.g. sequencing platforms)? Or is it just from the webpages?

• 0 views • link
written by evoclive • evoclive

Log in to answer this question.

More posts like this