GATK concordance

Answered

Follow

Priyadarshini Thirunavukkarasu

• November 11, 2021 11:08
• Edited

Hello

I am running GATK/4.2.2.0-foss-2018b-Java-1.8 in the cluster. After running the given below code, I am not able to find the output file (summary file). This is the link, where the code is given https://gatk.broadinstitute.org/hc/en-us/articles/4405451404699-Concordance#--summary. Please also find the log file below. Is the summary file required as input file to run the below script? Please advice

 

 gatk Concordance \
   -R /scicore/home/cichon/GROUP/memory_optimization/data/reference/gch38.fa \
   -eval /scicore/home/cichon/GROUP/memory_optimization/variants/filtered/sample1_affect.filtered.vcf \
   --truth /scicore/home/cichon/GROUP/memory_optimization/variants/filtered/NA12878.vcf.gz \
   --summary /scicore/home/cichon/GROUP/memory_optimization/variants/filtered/summary.tsv   

11:26:21.545 INFO NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/scicore/soft/apps/GATK/4.2.2.0-foss-2018b-Java-1.8/gatk-package-4.2.2.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
Nov 11, 2021 11:26:21 AM shaded.cloud_nio.com.google.auth.oauth2.ComputeEngineCredentials runningOnComputeEngine
INFO: Failed to detect whether we are running on Google Compute Engine.
11:26:21.681 INFO Concordance - ------------------------------------------------------------
11:26:21.682 INFO Concordance - The Genome Analysis Toolkit (GATK) v4.2.2.0
11:26:21.682 INFO Concordance - For support and documentation go to https://software.broadinstitute.org/gatk/
11:26:21.682 INFO Concordance - Executing as thirun0000@shi85.cluster.bc2.ch on Linux v3.10.0-1062.18.1.el7.x86_64 amd64
11:26:21.682 INFO Concordance - Java runtime: OpenJDK 64-Bit Server VM v1.8.0_212-b03
11:26:21.682 INFO Concordance - Start Date/Time: November 11, 2021 11:26:21 AM CET
11:26:21.682 INFO Concordance - ------------------------------------------------------------
11:26:21.682 INFO Concordance - ------------------------------------------------------------
11:26:21.683 INFO Concordance - HTSJDK Version: 2.24.1
11:26:21.683 INFO Concordance - Picard Version: 2.25.4
11:26:21.683 INFO Concordance - Built for Spark Version: 2.4.5
11:26:21.683 INFO Concordance - HTSJDK Defaults.COMPRESSION_LEVEL : 2
11:26:21.683 INFO Concordance - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
11:26:21.683 INFO Concordance - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
11:26:21.683 INFO Concordance - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
11:26:21.683 INFO Concordance - Deflater: IntelDeflater
11:26:21.684 INFO Concordance - Inflater: IntelInflater
11:26:21.684 INFO Concordance - GCS max retries/reopens: 20
11:26:21.684 INFO Concordance - Requester pays: disabled
11:26:21.684 INFO Concordance - Initializing engine
11:26:22.217 INFO FeatureManager - Using codec VCFCodec to read file file:///scicore/home/cichon/GROUP/memory_optimization/variants/filtered/NA12878.vcf.gz
11:26:22.497 INFO FeatureManager - Using codec VCFCodec to read file file:///scicore/home/cichon/GROUP/memory_optimization/variants/filtered/sample1_affect.filtered.vcf
11:26:22.663 INFO Concordance - Done initializing engine
11:26:22.672 INFO ProgressMeter - Starting traversal
11:26:22.672 INFO ProgressMeter - Current Locus Elapsed Minutes Records Processed Records/Minute
11:26:22.682 INFO Concordance - Shutting down engine
[November 11, 2021 11:26:22 AM CET] org.broadinstitute.hellbender.tools.walkers.validation.Concordance done. Elapsed time: 0.02 minutes.
Runtime.totalMemory()=559939584
java.lang.NullPointerException
at htsjdk.variant.variantcontext.VariantContextComparator.compare(VariantContextComparator.java:87)
at org.broadinstitute.hellbender.engine.AbstractConcordanceWalker$ConcordanceIterator.next(AbstractConcordanceWalker.java:192)
at org.broadinstitute.hellbender.engine.AbstractConcordanceWalker$ConcordanceIterator.next(AbstractConcordanceWalker.java:174)
at java.util.Iterator.forEachRemaining(Iterator.java:116)
at java.util.Spliterators$IteratorSpliterator.forEachRemaining(Spliterators.java:1801)
at java.util.stream.ReferencePipeline$Head.forEach(ReferencePipeline.java:580)
at org.broadinstitute.hellbender.engine.AbstractConcordanceWalker.traverse(AbstractConcordanceWalker.java:132)
at org.broadinstitute.hellbender.engine.GATKTool.doWork(GATKTool.java:1085)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.runTool(CommandLineProgram.java:140)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMainPostParseArgs(CommandLineProgram.java:192)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:211)
at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:160)
at org.broadinstitute.hellbender.Main.mainEntry(Main.java:203)
at org.broadinstitute.hellbender.Main.main(Main.java:289)
Using GATK jar /scicore/soft/apps/GATK/4.2.2.0-foss-2018b-Java-1.8/gatk-package-4.2.2.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /scicore/soft/apps/GATK/4.2.2.0-foss-2018b-Java-1.8/gatk-package-4.2.2.0-local.jar Concordance -R /scicore/home/cichon/GROUP/memory_optimization/data/reference/gch38.fa -eval /scicore/home/cichon/GROUP/memory_optimization/variants/filtered/sample1_affect.filtered.vcf --truth /scicore/home/cichon/GROUP/memory_optimization/variants/filtered/NA12878.vcf.gz --summary /scicore/home/cichon/GROUP/memory_optimization/variants/filtered/summary.tsv

Thanks

 

Priya

1

• 

  Genevieve Brandt (she/her)

  • November 12, 2021 21:24

  Hi Priya,

  There is an error message in your program log starting with "java.lang.NullPointerException". We took a look at the stack trace and found that it is coming from a contig mismatch from your truth and eval VCFs. Specifically, the eval VCF has contigs that the truth VCF doesn't have.

  You can inspect the headers of the VCFs to determine the contig differences and how to proceed. After you take a look, please let me know if you have further questions.

  I also created a GATK ticket so that our developers can have a better error message for this case. Here is the link: https://github.com/broadinstitute/gatk/issues/7562

  Best,

  Genevieve

  0

  Comment actions Permalink
• 

  Priyadarshini Thirunavukkarasu

  • November 14, 2021 20:19
  • Edited

  Hello

  Thank you. I figured that the contigs don't match and for the truth data set to operate over a defined set of regions in the evaluation sample, a interval list has to be specified. But, interval list is in a different format and the script stops running with the below error.

  A USER ERROR has occurred: Badly formed genome unclippedLoc: Contig NC_000001.11 given as location, but this contig isn't present in the Fasta sequence dictionary

  I have also pasted below the interval list which gives the above error

  
  @HD VN:1.5
  @SQ SN:NC_000001.11 LN:248956422
  @SQ SN:NT_187361.1 LN:175055
  @SQ SN:NT_187362.1 LN:32032
  @SQ SN:NT_187363.1 LN:127682
  @SQ SN:NT_187364.1 LN:66860
  @SQ SN:NT_187365.1 LN:40176
  @SQ SN:NT_187366.1 LN:42210
  @SQ SN:NT_187367.1 LN:176043
  @SQ SN:NT_187368.1 LN:40745

  Does anyone know how to format the interval list as shown in the below format:

  
  @HD VN:1.5
  @SQ SN:chr1 LN:248956422
  @SQ SN:chr2 LN:242193529
  @SQ SN:chr3 LN:198295559
  @SQ SN:chr4 LN:190214555
  @SQ SN:chr5 LN:181538259
  @SQ SN:chr6 LN:170805979
  @SQ SN:chr7 LN:159345973
  @SQ SN:chr8 LN:145138636
  @SQ SN:chr9 LN:138394717

  Thanks

  0

  Comment actions Permalink
• 

  Genevieve Brandt (she/her)

  • November 15, 2021 21:49

  Hi Priyadarshini Thirunavukkarasu,

  This is not a formatting issue with your interval list. The error message: "Contig NC_000001.11 given as location, but this contig isn't present in the Fasta sequence dictionary" is indicating that the first contig does not match the fasta file that you provided to the tool. This is another reference mismatch issue.

  Best,

  Genevieve

  0

  Comment actions Permalink
• 

  Priyadarshini Thirunavukkarasu

  • November 16, 2021 16:32

  Hello

  Sorry for the confusion.  I am trying to evaluate the filtered variants in filtered vcf with truth data set. The truth data vcf fies has chromosome number as chr. But, the filtered vcf has accession number in the vcf file to represent the chromosome number. The  "contig NC_000001.11" given in the error message is present in the filtered data set and it is not present in the truth data set (downloaded from GIAB website). Is there way to fix the issue, i.e. replace accession number or chromosome number in the vcf files. I found this issue  by creating intervals from truth data vcf files and filtered vcf data files.

  Thanks

  0

  Comment actions Permalink
• 

  Genevieve Brandt (she/her)

  • November 16, 2021 19:00

  The different contig names indicate that different references were used to map the reads and call variants. You'll need to find a truth set that uses the same reference version as your VCF.

  Please read this article for more information: https://gatk.broadinstitute.org/hc/en-us/articles/360035891131-Errors-about-input-files-having-missing-or-incompatible-contigs

  0

  Comment actions Permalink
• 

  Priyadarshini Thirunavukkarasu

  • November 16, 2021 20:36

  Hello

  I used GRCh38 as reference genome for mapping the exome data and the truth data set was from genome in bottle datasets and was also annotated with GRCh38. So, it is not clear, why there is an error and how to fix the error

  Thanks

   

  0

  Comment actions Permalink
• 

  Genevieve Brandt (she/her)

  • November 16, 2021 20:56

  You can check the contig lines in the header of your VCF: https://gatk.broadinstitute.org/hc/en-us/articles/360035531692-VCF-Variant-Call-Format

  You'll want to make sure that these match if they are truly from the same reference version.

  0

  Comment actions Permalink
• 

  Priyadarshini Thirunavukkarasu

  • November 17, 2021 11:23

  Thank you. So, the header of the filtered vcf data file and truth data vcf file should match?. If not, there are from different reference version. I will check the header of the vcf files. Is there way to visualize the header of vcf file which is gzipped. The truth vcf data set downloaded from the genome in bottle database is gzipped

  Thanks

   

  0

  Comment actions Permalink
• 

  Genevieve Brandt (she/her)

  • November 17, 2021 17:27

  The headers do not need to match exactly, only the contig lines.

  Here's a post on biostars about viewing zipped VCF files: https://www.biostars.org/p/454748/

  0

  Comment actions Permalink
• 

  Priyadarshini Thirunavukkarasu

  • November 17, 2021 21:32

  Thank you. In the truth data set downloaded from database, the contig position is given as chromosome number as Chr 1, Chr 2, Chr 3 and so on. In the filtered data set, the contigs are given as accession number. For example, NC_000001.11 is chromosome 1. The other difference, there are number of contigs in chromosome 1 in the filtered data set but in the truth data set, there is only contig for each chromosome

  Example in the truth data set
  
  @SQ SN:chr1 LN:248956422
  @SQ SN:chr2 LN:242193529
  @SQ SN:chr3 LN:198295559
  @SQ SN:chr4 LN:190214555
  @SQ SN:chr5 LN:181538259
  @SQ SN:chr6 LN:170805979
  @SQ SN:chr7 LN:159345973
  @SQ SN:chr8 LN:145138636
  @SQ SN:chr9 LN:138394717
  
  In filtered data set 
  
  @SQ SN:NC_000001.11 LN:248956422 (Chr1 as given in the above table, length (LN) is same as shown above)
  @SQ SN:NT_187361.1 LN:175055 (additional contig Chr 1, not seen in above table)
  @SQ SN:NT_187362.1 LN:32032 (additional contig Chr 1, not seen in above table)
  @SQ SN:NT_187363.1 LN:127682 (additional contig Chr 1, not seen in above table)
  @SQ SN:NT_187364.1 LN:66860 (additional contig Chr 1, not seen in above table)
  @SQ SN:NT_187365.1 LN:40176 (additional contig Chr 1, not seen in above table)
  @SQ SN:NT_187366.1 LN:42210 (additional contig Chr 1, not seen in above table)
  @SQ SN:NT_187367.1 LN:176043 (additional contig Chr 1, not seen in above table)
  @SQ SN:NT_187368.1 LN:40745 (additional contig Chr 1, not seen in above table)
  @SQ SN:NT_187369.1 LN:41717 (additional contig Chr 1, not seen in above table)
  @SQ SN:NC_000002.12 LN:242193529 (Chromosome 2, please see the table above where the length (LN is same)

  These additional contigs in the filtered data set leads to error "incompatible contigs" or contigs not present in fasta sequence dictionary. I downloaded the reference genome GRCH38, patch 13 from NCBI assembly database to align the sample files. The truth data set was shown to be annotated from genome in bottle database (GRCh38). Not sure, how to evaluate the filtered data sets as there are additional contigs in the filtered data sets compared to the truth data sets. Please advice

  Thanks

  0

  Comment actions Permalink
• 

  Genevieve Brandt (she/her)

  • November 18, 2021 21:24

  We have an article on reference versions here: https://gatk.broadinstitute.org/hc/en-us/articles/360035890951-Human-genome-reference-builds-GRCh38-or-hg38-b37-hg19

  0

  Comment actions Permalink
• 

  nagam surya

  • May 12, 2022 22:13

  I have two vcf files to comapre and understand the concordance between them. I do not have a summary tsv, cause it is of no use to me. So, can't I just compare the two VCF files without the summary table using GATK?

  0

  Comment actions Permalink
• 

  Genevieve Brandt (she/her)

  • July 21, 2022 23:25

  nagam surya did your other post answer this question?

  0

  Comment actions Permalink

Please sign in to leave a comment.