NCBI chr reference name different from dbSNP chr names, is that usable by GATK?

Hello,

I am trying to perform base recalibration on a non-human genome. I've downloaded the dbSNP vcf file from Ensembl (also downloaded another version from NCBI), and for both files the chr naming system is "1, 2, 3...". On the other hand, the NCBI reference file I used named chromosomes differently, chr 1 = gi|966749131|ref|NC_006088.4|, chr2 = gi|966749130|ref|NC_006089.4|...

However, BaseRecalibration seemed to run fine. I am wondering if GATK is capable of internally recognizing the NCBI chromosome names or if their algorithm generated a completely wrong recal table and recalibrated all the scores wrong. I've also attached an output of my AnalyzeCovariates plot in case that helps. Thanks!

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie
    edited February 2016

    That's odd because GATK should notice the mismatching sequence dictionaries and quit with an error. Can you post your command line and the terminal output?

  • angezouangezou Member
    edited February 2016

    Hi @Geraldine_VdAuwera, here is the command line:

    gatk_compute="java -jar /hpf/tools/centos6/gatk/3.5.0/GenomeAnalysisTK.jar"
    $gatk_compute -T BaseRecalibrator -nct 4 -R $location/ref_ensembl/Gallus_gallus.toplevel.fa \
    -I $location/globus/${str}_realigned_all.bam -knownSites $location/rmcontigs.vcf -o $location/${str}_recal.table && \
    $gatk_compute -T BaseRecalibrator -nct 4 -R $location/ref_ensembl/Gallus_gallus.toplevel.fa \
    -I $location/globus/${str}_realigned_all.bam -knownSites $location/rmcontigs.vcf -BQSR $location/${str}_recal.table -o $location/${str}_postrecal.table && \
    $gatk_compute -T AnalyzeCovariates -R $location/ref_ensembl/Gallus_gallus.toplevel.fa \
    -before $location/${str}_recal.table -after $location/${str}_postrecal.table -plots ${str}_recal_plots.pdf

    Portions of the terminal output:

    Picked up _JAVA_OPTIONS: -Xmx40g
    INFO 19:05:45,662 HelpFormatter - --------------------------------------------------------------------------------
    INFO 19:05:45,667 HelpFormatter - The Genome Analysis Toolkit (GATK) v3.5-0-g36282e4, Compiled 2015/11/25 04:03:56
    INFO 19:05:45,668 HelpFormatter - Copyright (c) 2010 The Broad Institute
    INFO 19:05:45,668 HelpFormatter - For support and documentation go to http://www.broadinstitute.org/gatk
    INFO 19:05:45,672 HelpFormatter - Program Args: -T BaseRecalibrator -nct 4 -R /hpf/largeprojects/jparkin/angezou/ref_ensembl/Gallus_gallus.toplevel.fa -I /hpf/largeprojects/jparkin/angezou/151105_realigned_merged.bam -knownSites /hpf/largeprojects/jparkin/angezou/rmcontigs.vcf -o /hpf/largeprojects/jparkin/angezou/151105_recal.table
    INFO 19:05:45,679 HelpFormatter - Executing as angezou@r3b-5 on Linux 2.6.32-358.18.1.el6.x86_64 amd64; OpenJDK 64-Bit Server VM 1.7.0_55-mockbuild_2014_04_16_12_11-b00.
    INFO 19:05:45,679 HelpFormatter - Date/Time: 2016/02/01 19:05:45
    INFO 19:05:45,679 HelpFormatter - --------------------------------------------------------------------------------
    INFO 19:05:45,679 HelpFormatter - --------------------------------------------------------------------------------
    INFO 19:05:46,054 GenomeAnalysisEngine - Strictness is SILENT
    INFO 19:05:46,994 GenomeAnalysisEngine - Downsampling Settings: No downsampling
    INFO 19:05:47,002 SAMDataSource$SAMReaders - Initializing SAMRecords in serial
    INFO 19:05:48,027 SAMDataSource$SAMReaders - Done initializing BAM readers: total time 1.00
    INFO 19:05:48,957 MicroScheduler - Running the GATK in parallel mode with 4 total threads, 4 CPU thread(s) for each of 1 data thread(s), of 40 processors available on this machine
    INFO 19:05:49,284 GenomeAnalysisEngine - Preparing for traversal over 1 BAM files
    INFO 19:05:49,289 GenomeAnalysisEngine - Done preparing for traversal
    INFO 19:05:49,290 ProgressMeter - [INITIALIZATION COMPLETE; STARTING PROCESSING]
    INFO 19:05:49,290 ProgressMeter - | processed | time | per 1M | | total | remaining
    INFO 19:05:49,290 ProgressMeter - Location | reads | elapsed | reads | completed | runtime | runtime
    INFO 19:05:49,312 BaseRecalibrator - The covariates being used here:
    INFO 19:05:49,312 BaseRecalibrator - ReadGroupCovariate
    INFO 19:05:49,313 BaseRecalibrator - QualityScoreCovariate
    INFO 19:05:49,313 BaseRecalibrator - ContextCovariate
    INFO 19:05:49,313 ContextCovariate - Context sizes: base substitution model 2, indel substitution model 3
    INFO 19:05:49,313 BaseRecalibrator - CycleCovariate
    INFO 19:05:49,746 ReadShardBalancer$1 - Loading BAM index data
    INFO 19:05:49,747 ReadShardBalancer$1 - Done loading BAM index data
    INFO 19:06:19,427 ProgressMeter - 1:923466 500007.0 30.0 s 60.0 s 0.1% 9.4 h 9.4 h
    INFO 19:06:49,429 ProgressMeter - 1:1562306 900012.0 60.0 s 66.0 s 0.1% 11.2 h 11.2 h
    INFO 19:07:19,550 ProgressMeter - 1:2193399 1400022.0 90.0 s 64.0 s 0.2% 11.9 h 11.9 h
    .
    .
    .
    .
    .
    INFO 05:06:04,549 ProgressMeter - AADN03012205.1:148 7.51680518E8 10.0 h 47.0 s 100.0% 10.0 h 0.0 s
    INFO 05:06:34,551 ProgressMeter - AADN03012205.1:148 7.51680518E8 10.0 h 47.0 s 100.0% 10.0 h 0.0 s
    INFO 05:07:04,552 ProgressMeter - AADN03012205.1:148 7.51680518E8 10.0 h 47.0 s 100.0% 10.0 h 0.0 s
    INFO 05:07:34,553 ProgressMeter - AADN03012205.1:148 7.51680518E8 10.0 h 48.0 s 100.0% 10.0 h 0.0 s
    INFO 05:07:38,738 BaseRecalibrator - Calculating quantized quality scores...
    INFO 05:07:38,885 BaseRecalibrator - Writing recalibration report...
    INFO 05:07:40,261 BaseRecalibrator - ...done!
    INFO 05:07:40,261 BaseRecalibrator - BaseRecalibrator was able to recalibrate 751496765 reads
    INFO 05:07:40,263 ProgressMeter - done 7.51681338E8 10.0 h 48.0 s 100.0% 10.0 h 0.0 s
    INFO 05:07:40,264 ProgressMeter - Total runtime 36110.97 secs, 601.85 min, 10.03 hours
    INFO 05:07:40,264 MicroScheduler - 182975543 reads were filtered out during the traversal out of approximately 934656881 total reads (19.58%)
    INFO 05:07:40,265 MicroScheduler - -> 0 reads (0.00% of total) failing BadCigarFilter
    INFO 05:07:40,265 MicroScheduler - -> 0 reads (0.00% of total) failing DuplicateReadFilter
    INFO 05:07:40,266 MicroScheduler - -> 0 reads (0.00% of total) failing FailsVendorQualityCheckFilter
    INFO 05:07:40,266 MicroScheduler - -> 0 reads (0.00% of total) failing MalformedReadFilter
    INFO 05:07:40,267 MicroScheduler - -> 0 reads (0.00% of total) failing MappingQualityUnavailableFilter
    INFO 05:07:40,267 MicroScheduler - -> 182192521 reads (19.49% of total) failing MappingQualityZeroFilter
    INFO 05:07:40,267 MicroScheduler - -> 783022 reads (0.08% of total) failing NotPrimaryAlignmentFilter
    INFO 05:07:40,268 MicroScheduler - -> 0 reads (0.00% of total) failing UnmappedReadFilter
    INFO 05:07:42,392 HttpMethodDirector - I/O exception (java.net.ConnectException) caught when processing request: Network is unreachable
    INFO 05:07:42,393 HttpMethodDirector - Retrying request
    INFO 05:07:42,395 HttpMethodDirector - I/O exception (java.net.ConnectException) caught when processing request: Network is unreachable
    INFO 05:07:42,396 HttpMethodDirector - Retrying request
    INFO 05:07:42,398 HttpMethodDirector - I/O exception (java.net.ConnectException) caught when processing request: Network is unreachable
    INFO 05:07:42,398 HttpMethodDirector - Retrying request
    INFO 05:07:42,400 HttpMethodDirector - I/O exception (java.net.ConnectException) caught when processing request: Network is unreachable
    INFO 05:07:42,400 HttpMethodDirector - Retrying request
    INFO 05:07:42,402 HttpMethodDirector - I/O exception (java.net.ConnectException) caught when processing request: Network is unreachable
    INFO 05:07:42,402 HttpMethodDirector - Retrying request
    Picked up _JAVA_OPTIONS: -Xmx40g
    INFO 05:07:45,518 HelpFormatter - --------------------------------------------------------------------------------
    INFO 05:07:45,521 HelpFormatter - The Genome Analysis Toolkit (GATK) v3.5-0-g36282e4, Compiled 2015/11/25 04:03:56
    INFO 05:07:45,521 HelpFormatter - Copyright (c) 2010 The Broad Institute
    INFO 05:07:45,521 HelpFormatter - For support and documentation go to http://www.broadinstitute.org/gatk
    INFO 05:07:45,525 HelpFormatter - Program Args: -T BaseRecalibrator -nct 4 -R /hpf/largeprojects/jparkin/angezou/ref_ensembl/Gallus_gallus.toplevel.fa -I /hpf/largeprojects/jparkin/angezou/151105_realigned_merged.bam -knownSites /hpf/largeprojects/jparkin/angezou/rmcontigs.vcf -BQSR /hpf/largeprojects/jparkin/angezou/151105_recal.table -o /hpf/largeprojects/jparkin/angezou/151105_postrecal.table
    INFO 05:07:45,547 HelpFormatter - Executing as angezou@r3b-5 on Linux 2.6.32-358.18.1.el6.x86_64 amd64; OpenJDK 64-Bit Server VM 1.7.0_55-mockbuild_2014_04_16_12_11-b00.
    INFO 05:07:45,547 HelpFormatter - Date/Time: 2016/02/02 05:07:45
    INFO 05:07:45,547 HelpFormatter - --------------------------------------------------------------------------------
    INFO 05:07:45,548 HelpFormatter - --------------------------------------------------------------------------------
    INFO 05:07:45,685 GenomeAnalysisEngine - Strictness is SILENT
    INFO 05:07:47,040 ContextCovariate - Context sizes: base substitution model 2, indel substitution model 3
    INFO 05:07:47,069 GenomeAnalysisEngine - Downsampling Settings: No downsampling
    INFO 05:07:47,076 SAMDataSource$SAMReaders - Initializing SAMRecords in serial
    INFO 05:07:47,637 SAMDataSource$SAMReaders - Done initializing BAM readers: total time 0.56
    INFO 05:07:48,336 MicroScheduler - Running the GATK in parallel mode with 4 total threads, 4 CPU thread(s) for each of 1 data thread(s), of 40 processors available on this machine
    INFO 05:07:48,495 GenomeAnalysisEngine - Preparing for traversal over 1 BAM files
    INFO 05:07:48,500 GenomeAnalysisEngine - Done preparing for traversal
    INFO 05:07:48,500 ProgressMeter - [INITIALIZATION COMPLETE; STARTING PROCESSING]
    INFO 05:07:48,500 ProgressMeter - | processed | time | per 1M | | total | remaining
    INFO 05:07:48,501 ProgressMeter - Location | reads | elapsed | reads | completed | runtime | runtime
    INFO 05:07:48,526 BaseRecalibrator - The covariates being used here:
    INFO 05:07:48,527 BaseRecalibrator - ReadGroupCovariate
    INFO 05:07:48,527 BaseRecalibrator - QualityScoreCovariate
    INFO 05:07:48,527 BaseRecalibrator - ContextCovariate
    INFO 05:07:48,527 ContextCovariate - Context sizes: base substitution model 2, indel substitution model 3
    INFO 05:07:48,528 BaseRecalibrator - CycleCovariate
    INFO 05:07:48,556 ReadShardBalancer$1 - Loading BAM index data
    INFO 05:07:48,557 ReadShardBalancer$1 - Done loading BAM index data
    INFO 05:08:18,522 ProgressMeter - Starting 0.0 30.0 s 49.6 w 100.0% 30.0 s 0.0 s
    INFO 05:08:48,524 ProgressMeter - 1:189539 0.0 60.0 s 99.2 w 0.0% 92.1 h 92.0 h
    INFO 05:09:48,526 ProgressMeter - 1:546987 200003.0 120.0 s 10.0 m 0.1% 63.8 h 63.8 h
    INFO 05:10:48,527 ProgressMeter - 1:897774 400006.0 3.0 m 7.5 m 0.1% 58.3 h 58.3 h
    .
    .
    .
    .
    .
    INFO 01:53:12,580 ProgressMeter - AADN03011097.1:672 7.51357566E8 44.8 h 3.6 m 99.9% 44.8 h 93.0 s
    INFO 01:54:12,582 ProgressMeter - AADN03012205.1:148 7.51680518E8 44.8 h 3.6 m 100.0% 44.8 h 0.0 s
    INFO 01:55:12,585 ProgressMeter - AADN03012205.1:148 7.51680518E8 44.8 h 3.6 m 100.0% 44.8 h 0.0 s
    INFO 01:56:12,587 ProgressMeter - AADN03012205.1:148 7.51680518E8 44.8 h 3.6 m 100.0% 44.8 h 0.0 s
    INFO 01:56:54,566 BaseRecalibrator - Calculating quantized quality scores...
    INFO 01:56:54,679 BaseRecalibrator - Writing recalibration report...
    INFO 01:56:57,002 BaseRecalibrator - ...done!
    INFO 01:56:57,002 BaseRecalibrator - BaseRecalibrator was able to recalibrate 751496765 reads
    INFO 01:56:57,003 ProgressMeter - done 7.51681338E8 44.8 h 3.6 m 100.0% 44.8 h 0.0 s
    INFO 01:56:57,003 ProgressMeter - Total runtime 161348.50 secs, 2689.14 min, 44.82 hours
    INFO 01:56:57,004 MicroScheduler - 182975543 reads were filtered out during the traversal out of approximately 934656881 total reads (19.58%)
    INFO 01:56:57,004 MicroScheduler - -> 0 reads (0.00% of total) failing BadCigarFilter
    INFO 01:56:57,004 MicroScheduler - -> 0 reads (0.00% of total) failing DuplicateReadFilter
    INFO 01:56:57,004 MicroScheduler - -> 0 reads (0.00% of total) failing FailsVendorQualityCheckFilter
    INFO 01:56:57,004 MicroScheduler - -> 0 reads (0.00% of total) failing MalformedReadFilter
    INFO 01:56:57,005 MicroScheduler - -> 0 reads (0.00% of total) failing MappingQualityUnavailableFilter
    INFO 01:56:57,005 MicroScheduler - -> 182192521 reads (19.49% of total) failing MappingQualityZeroFilter
    INFO 01:56:57,005 MicroScheduler - -> 783022 reads (0.08% of total) failing NotPrimaryAlignmentFilter
    INFO 01:56:57,005 MicroScheduler - -> 0 reads (0.00% of total) failing UnmappedReadFilter
    INFO 01:56:58,424 HttpMethodDirector - I/O exception (java.net.ConnectException) caught when processing request: Network is unreachable
    INFO 01:56:58,424 HttpMethodDirector - Retrying request
    INFO 01:56:58,428 HttpMethodDirector - I/O exception (java.net.ConnectException) caught when processing request: Network is unreachable
    INFO 01:56:58,429 HttpMethodDirector - Retrying request
    INFO 01:56:58,432 HttpMethodDirector - I/O exception (java.net.ConnectException) caught when processing request: Network is unreachable
    INFO 01:56:58,432 HttpMethodDirector - Retrying request
    INFO 01:56:58,434 HttpMethodDirector - I/O exception (java.net.ConnectException) caught when processing request: Network is unreachable
    INFO 01:56:58,434 HttpMethodDirector - Retrying request
    INFO 01:56:58,436 HttpMethodDirector - I/O exception (java.net.ConnectException) caught when processing request: Network is unreachable
    INFO 01:56:58,436 HttpMethodDirector - Retrying request

  • @Geraldine_VdAuwera it's ok, since you say it should have exited with error, then I will just bootstrap it

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Actually my point is that the sequence dictionary checks are very thorough, and I see from your command line that you didn't disable them (some people do). So the fact that it worked suggests to me that the contig naming scheme might be correct. Please don't be offended, but did you check the files' actual sequence dictionaries?

  • Definitely not offended, where would I go to check the files' sequence dictionaries?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    In the vcf and bam it's the @SQ lines in the header, for the reference it's the dict file.

Sign In or Register to comment.