Notice:
If you happen to see a question you know the answer to, please do chime in and help your fellow community members. We encourage our fourm members to be more involved, jump in and help out your fellow researchers with their questions. GATK forum is a community forum and helping each other with using GATK tools and research is the cornerstone of our success as a genomics research community.We appreciate your help!

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.
Attention:
We will be out of the office on November 11th and 13th 2019, due to the U.S. holiday(Veteran's day) and due to a team event(Nov 13th). We will return to monitoring the GATK forum on November 12th and 14th respectively. Thank you for your patience.

ReadbackedPhasing Error

Run ReadBackedPhasing to have following error, anyone has any suggestions?

SamDataSource$SAMReaders - Done initializing. BAM readers: total time 0.01
RMDTrackBuilder - Writing Tribble index to disk for file ...x.vcf.idx

A USER error has occurred

ERROR Message: I/O error loading or writing tribble indexc file for x.vcf

Best Answers

Answers

  • joeljoel USAMember

    you mean x.vcf.idx file? The tool tries to create the file .idx, I check the file, it is empty. Looks like the program can create .idx file but fail to write anything to it,

  • joeljoel USAMember

    Regarding idx file, the line RMDTrackBuilder - Writing Tribble index to disk for file ...x.vcf.idx

    I check the file, it is empty. Is that expected? In which case RMDTrackBuilder fails to write to vcf.idx file?

  • joeljoel USAMember

    Thank you pdexheimer. I will try to skip index creation. I ran on Windows platform, there shouldn't be a permission thing. The file path is not that long, less than 50 characters, and the file name is shorter than 8 chars. I wonder if anyone has experienced this issue before and what the workaround will be.

  • pdexheimerpdexheimer Member ✭✭✭✭

    I ran on Windows platform

    Not to get all "arrogant Mac/Unix user" on you, but that's probably your problem. Windows can/does have some weird little incompatibilities regarding file access in particular, and I'm just about certain that GATK is officially unsupported on Windows

  • joeljoel USAMember

    ok, I can run it on Linux. But for now I bypass the index file, now I got an error message

    Contig A does not have a length field,

    in my cvf file, I have this line

    contig=<ID=A,length=3644>

    So, is the error message about something else?

  • joeljoel USAMember

    It does run on Linux (it creates index file), but got the same error message. Any idea?

  • joeljoel USAMember

    By the way, I tried to search where this error message comes from. It looks like a user exception somewhere in the code that makes the program quit. By searching the source code depot, I couldn't find any place that raises the exception. Search by words "does not have a" comes back nothing. Weird.

  • joeljoel USAMember

    INFO 19:16:17,652 HelpFormatter - --------------------------------------------------------------------------------
    INFO 19:16:17,654 HelpFormatter - The Genome Analysis Toolkit (GATK) v3.5-0-g36282e4, Compiled 2015/11/25 04:03:56
    INFO 19:16:17,654 HelpFormatter - Copyright (c) 2010 The Broad Institute
    INFO 19:16:17,654 HelpFormatter - For support and documentation go to http://www.broadinstitute.org/gatk
    INFO 19:16:17,657 HelpFormatter - Program Args: -T ReadBackedPhasing -I /home/macjs/Downloads/gatk/001.bam -V /home/macjs/Downloads/gatk/001.vcf -R /home/macjs/Downloads/gatk/001_ref.fasta -o /home/macjs/Downloads/gatk/.vcf
    INFO 19:16:17,663 HelpFormatter - Executing as [email protected] on Linux 3.13.0-32-generic amd64; OpenJDK 64-Bit Server VM 1.7.0_75-b13.
    INFO 19:16:17,663 HelpFormatter - Date/Time: 2015/12/14 19:16:17
    INFO 19:16:17,663 HelpFormatter - --------------------------------------------------------------------------------
    INFO 19:16:17,663 HelpFormatter - --------------------------------------------------------------------------------
    INFO 19:16:17,722 GenomeAnalysisEngine - Strictness is SILENT
    INFO 19:16:17,775 GenomeAnalysisEngine - Downsampling Settings: Method: BY_SAMPLE, Target Coverage: 1000
    INFO 19:16:17,782 SAMDataSource$SAMReaders - Initializing SAMRecords in serial
    INFO 19:16:17,803 SAMDataSource$SAMReaders - Done initializing BAM readers: total time 0.02
    INFO 19:16:18,934 GATKRunReport - Uploaded run statistics report to AWS S3

    ERROR ------------------------------------------------------------------------------------------
    ERROR A USER ERROR has occurred (version 3.5-0-g36282e4):
    ERROR
    ERROR This means that one or more arguments or inputs in your command are incorrect.
    ERROR The error message below tells you what is the problem.
    ERROR
    ERROR If the problem is an invalid argument, please check the online documentation guide
    ERROR (or rerun your command with --help) to view allowable command-line arguments for this tool.
    ERROR
    ERROR Visit our website and forum for extensive documentation and answers to
    ERROR commonly asked questions http://www.broadinstitute.org/gatk
    ERROR
    ERROR Please do NOT post this error to the GATK forum unless you have really tried to fix it yourself.
    ERROR
    ERROR MESSAGE: Contig A does not have a length field.
    ERROR ------------------------------------------------------------------------------------------
  • joeljoel USAMember

    is this bam file issue or is it vcf file issue?

  • pdexheimerpdexheimer Member ✭✭✭✭

    I don't know exactly where that error is getting triggered, it's not one that I've seen before. It could be in the indices for the vcf or bam, or even in the indices for the reference. If it were me, I'd probably examine the fai and dict files to make sure they're well-formed (alternatively, recreate them), delete and manually remake the .bai file, and then delete and let GATK auto-generate the .vcf.idx file.

    Another note that's likely not related - your output file doesn't have a base name. It's just the (hidden) file ".vcf"

  • joeljoel USAMember

    Thank for pointing out the output file name, by the way, what is the java version that GATK runs on? 1.8 or 1.7?

  • joeljoel USAMember

    but for creating dict file, we have to use picard, and picard has to run on 1.8, this seems to conflict.

  • SheilaSheila Broad InstituteMember, Broadie admin

    @joel
    Hi,

    Yes, we realize this is frustrating. GATK will soon use Java-1.8. In the meantime, you can easily switch between Java versions by using this command structure:

    /usr/libexec/java_home -v 1.8.0_60 --exec java -jar Picard.jar...

    So, you will first give the location of the Java home, then use -v to specify the Java version to use. --exec tells the terminal to use the previously specified Java version to execute the command.

    There will be documentation on this soon.

    -Sheila

  • joeljoel USAMember

    To dig a bit more deep, the error of Contig does not have a length field comes from htsjdk code line:

    public SAMSequenceRecord getSAMSequenceRecord() {
    final String lengthString = this.getGenericFieldValue("length");
    if (lengthString == null) throw new TribbleException("Contig " + this.getID() + " does not have a length field.");

    which is called by RMDTrackBuilder.java validateVariantAgainstSequenceDictionary(final String name, final String descriptorName, final AbstractFeatureReader reader, final SAMSequenceDictionary dict )

    in vcfContigRecords.add(contig.getSAMSequenceRecord());

    It seems to load contig sequence but fails. I thought the contig sequence is from ref fasta file, anyone has experience on this before?

  • joeljoel USAMember

    I will post it in a new thread about this problem to get more attention from the forum

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @joel Posting the same thing several times will not get you more attention, or not the positive kind anyway. Repeat posts cause clutter and impedes our ability to respond efficiently to questions.

    Have you tried what pdexheimer suggested earlier? It's an excellent suggestion.

    f it were me, I'd probably examine the fai and dict files to make sure they're well-formed (alternatively, recreate them), delete and manually remake the .bai file, and then delete and let GATK auto-generate the .vcf.idx file.

    The most likely explanation is that one of your files is indeed malformed, so digging through the code won't help you half as much as digging through the files themselves.

  • joeljoel USAMember

    The error message is different from the error originated for this thread, that is why I posted in a new thread. Not trying to repeat the same post everywhere. I did fix the output file name, but that is nothing to do with this error.

    Picard runs on a java 1.8 while gatk runs on 1.7. I run everything on 1.7. Picard is used to generate dict file which is a text file, so doesn't really matter where to generate dict file. I did regenerate fai/dict and let gatk auto-generate vcf.idx on Linux OS. Since my files also are generated from software, I don;t think I can figure out what is wrong until the source code that raises the exception can be investigated.

    Anyway, just like to share what i learned and see if anyone on this forum has seen this before.

  • SheilaSheila Broad InstituteMember, Broadie admin

    @joel
    Hi,

    So, your VCF does have contig A with a length. Can you post the BAM header and FASTA dict file?

    Thanks,
    Sheila

  • joeljoel USAMember

    dict file

    @HD VN:1.5 SO:unsorted
    @SQ SN:A LN:3644 M5:81dbc3cf5c46c4dbea5256b3aeee4998 UR:file:/home/macjs/Downloads/gatk/001_ref.fasta
    @SQ SN:B LN:3373 M5:f2257d9c41af8dbc5ad6518f9fdf67a0 UR:file:/home/macjs/Downloads/gatk/001_ref.fasta
    @SQ SN:C LN:3427 M5:af595de3b819d394764ee8d5410ecf99 UR:file:/home/macjs/Downloads/gatk/001_ref.fasta
    @SQ SN:DRB1 LN:18814 M5:88d5b128438dbcb7010d3bf55a8a5795 UR:file:/home/macjs/Downloads/gatk/001_ref.fasta
    @SQ SN:DQB1 LN:7747 M5:b62bf68c2e497dde588ffa05a1c4d59d UR:file:/home/macjs/Downloads/gatk/001_ref.fasta
    @SQ SN:DPB1 LN:11538 M5:76c722c4c3882ba3a12128cf78282acc UR:file:/home/macjs/Downloads/gatk/001_ref.fasta
    @SQ SN:DRB3 LN:16872 M5:27ba6b3bc95fc12ba408f61aa13e062b UR:file:/home/macjs/Downloads/gatk/001_ref.fasta

    and bam header:

    @HD VN:1.4 SO:coordinate
    @RG ID:HTSRead PL:iontorrent SM:TER222
    @SQ SN:A LN:3644
    @SQ SN:B LN:3373
    @SQ SN:C LN:3427
    @SQ SN:DRB1 LN:18814
    @SQ SN:DQB1 LN:7747
    @SQ SN:DPB1 LN:11538
    @SQ SN:DRB3 LN:16872

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    There you go, the headers don't match. You say you regenerated the dictionary file; can you post the command lines you used and the output log? Are you sure the run completed successfully? Consider checking the file creation timestamp.

  • joeljoel USAMember

    java -jar picard.jar CreateSequenceDictionary R=001_ref.fasta O=001_ref.dict

  • joeljoel USAMember

    What it should be? THe dictionary file is generated from fasta file while the bam file header is from bam file with the same reference fasta file. You see the contig name and length are the same, except dict file has more info (M5 must be the checksum and file name itself). Does M5/file name matters?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Ack, sorry -- I misread and was looking at version numbers, which is irrelevant here, because I was thinking about someone else's issue. Never mind my previous comment.

    What does the VCF header look like?

  • joeljoel USAMember

    fileformat=VCFv4.2

    FILTER=<ID=PASS,Description="All filters passed">

    samtoolsVersion=1.2??1.2.1

    samtoolsCommand=samtools mpileup -uf ../result/001references.fasta bams001.bam

    reference=file://../result001_references.fasta

    contig=<ID=A,length=3644>

    contig=<ID=B,length=3373>

    contig=<ID=C,length=3427>

    contig=<ID=DRB1,length=18814>

    contig=<ID=DQB1,length=7747>

    contig=<ID=DPB1,length=11538>

    contig=<ID=DRB3,length=16872>

    ALT=<ID=X,Description="Represents allele(s) other than observed.">

    INFO=<ID=INDEL,Number=0,Type=Flag,Description="Indicates that the variant is an INDEL.">

    INFO=<ID=IDV,Number=1,Type=Integer,Description="Maximum number of reads supporting an indel">

    INFO=<ID=IMF,Number=1,Type=Float,Description="Maximum fraction of reads supporting an indel">

    INFO=<ID=DP,Number=1,Type=Integer,Description="Raw read depth">

    INFO=<ID=VDB,Number=1,Type=Float,Description="Variant Distance Bias for filtering splice-site artefacts in RNA-seq data (bigger is better)",Version="3">

    INFO=<ID=RPB,Number=1,Type=Float,Description="Mann-Whitney U test of Read Position Bias (bigger is better)">

    INFO=<ID=MQB,Number=1,Type=Float,Description="Mann-Whitney U test of Mapping Quality Bias (bigger is better)">

    INFO=<ID=BQB,Number=1,Type=Float,Description="Mann-Whitney U test of Base Quality Bias (bigger is better)">

    INFO=<ID=MQSB,Number=1,Type=Float,Description="Mann-Whitney U test of Mapping Quality vs Strand Bias (bigger is better)">

    INFO=<ID=SGB,Number=1,Type=Float,Description="Segregation based metric.">

    INFO=<ID=MQ0F,Number=1,Type=Float,Description="Fraction of MQ0 reads (smaller is better)">

    FORMAT=<ID=PL,Number=G,Type=Integer,Description="List of Phred-scaled genotype likelihoods">

    FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">

    INFO=<ID=ICB,Number=1,Type=Float,Description="Inbreeding Coefficient Binomial test (bigger is better)">

    INFO=<ID=HOB,Number=1,Type=Float,Description="Bias in the number of HOMs number (smaller is better)">

    INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes for each ALT allele, in the same order as listed">

    INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">

    INFO=<ID=DP4,Number=4,Type=Integer,Description="Number of high-quality ref-forward , ref-reverse, alt-forward and alt-reverse bases">

    INFO=<ID=MQ,Number=1,Type=Integer,Description="Average mapping quality">

    bcftools_callVersion=1.2??1.2.1

    bcftools_callCommand=call -mv

  • joeljoel USAMember

    by the way, is there any instruction how to build and debug GenomeAnalysisTK module? I am new to java development, any help?

    Issue · Github
    by Sheila

    Issue Number
    424
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    vdauwera
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Your file has unexpected characters (like the line samtoolsVersion=1.2??1.2.1 ) which makes me think that it is somehow malformed so the parsing code is choking on it. This could be a side effect of working on Windows, or something else. I'm afraid we can't help much with that kind of problem.

    There are some developerٍ oriented docs in the guide (see DevZone) including compilation instructions. Debugging is not covered.

    Good luck!

  • joeljoel USAMember

    Actually if I removed the line

    contig=<ID=A,length=3644>

    from the vcf file, the program proceeds without the error out. So somehow the program does not like the field of Contig in the header session.

  • joeljoel USAMember

    By the way, could you point me to the document that describe the new HP field for phased info?

Sign In or Register to comment.