Bug Bulletin: The GenomeLocPArser error in SplitNCigarReads has been fixed; if you encounter it, use the latest nightly build.

Reg. Adding header line to Ensembl vcf file

sswaminathansswaminathan USAPosts: 16Member

Hello

I am working with canine genomes and donwloaded the reference file (ftp://ftp.ensembl.org/pub/release-73/fasta/canis_familiaris/dna/Canis_familiaris.CanFam3.1.73.dna.toplevel.fa.gz) and variation VCF file (ftp://ftp.ensembl.org/pub/release-73/variation/vcf/canis_familiaris/Canis_familiaris.vcf.gz) from Ensembl. I was able to index the reference files and perform the alignment and Mark duplicate steps of the pipeline for the canine genomes.

I extracted SNPs that were marked either deleteions or duplications from the variation VCF file using grep -E "^#|deletion|insertion" Canis_familiaris.vcf > Canis_familiaris.indels.vcf to use in the Indel Realigner steps. The header of the file is as follows:

##fileformat=VCFv4.1
    ##fileDate=20130826
    ##source=ensembl;version=73;url=http://e73.ensembl.org/canis_familiaris
    ##reference=ftp://ftp.ensembl.org/pub/release-73/fasta/canis_familiaris/dna/
    ##INFO=<ID=TSA,Number=0,Type=String,Description="Type of sequence alteration. Child of term sequence_alteration as defined by the sequence ontology project.">
    ##INFO=<ID=E_MO,Number=0,Type=Flag,Description="Multiple_observations.http://www.ensembl.org/info/docs/variation/data_description.html#evidence_status">
    ##INFO=<ID=E_ESP,Number=0,Type=Flag,Description="Exome_Sequencing_Project.http://www.ensembl.org/info/docs/variation/data_description.html#evidence_status">
    ##INFO=<ID=E_1000G,Number=0,Type=Flag,Description="1000Genomes.http://www.ensembl.org/info/docs/variation/data_description.html#evidence_status">
    ##INFO=<ID=E_HM,Number=0,Type=Flag,Description="HapMap.http://www.ensembl.org/info/docs/variation/data_description.html#evidence_status">
    ##INFO=<ID=E_Freq,Number=0,Type=Flag,Description="Frequency.http://www.ensembl.org/info/docs/variation/data_description.html#evidence_status">
    ##INFO=<ID=E_C,Number=0,Type=Flag,Description="Cited.http://www.ensembl.org/info/docs/variation/data_description.html#evidence_status">
    ##INFO=<ID=dbSNP_131,Number=0,Type=Flag,Description="Variants (including SNPs and indels) imported from dbSNP [remapped from build CanFam2.0]">
    1   8252    rs8962840   C   CT  .   .   dbSNP_131;TSA=insertion
    1   9702    rs8457244   CA  C   .   .   dbSNP_131;TSA=deletion
    1   18289   rs8444620   GT  G   .   .   dbSNP_131;TSA=deletion
    1   36381   rs8471229   GT  G   .   .   dbSNP_131;TSA=deletion
    1   52452   rs8469855   AT  A   .   .   dbSNP_131;TSA=deletion
    1   55767   rs8285977   G   GA  .   .   dbSNP_131;TSA=insertion
    1   60065   rs8723538   TC  T   .   .   dbSNP_131;TSA=deletion
    1   62067   rs8395051   TA  T   .   .   dbSNP_131;TSA=deletion

However, when I run the filrst step of Indel Realignment -T RealignerTargetCreator -nt 16 -R canFam3.1_bwa075/Canis_familiaris.CanFam3.1.73.dna.toplevel.fa -I AF23.rg.prmdup.bam -l INFO -S SILENT -o AF23.indel.intervals --known canFam3.1_bwa075/Canis_familiaris.indels.vcf

I get the error message:

##### ERROR A USER ERROR has occurred (version 2.7-2-g6bda569):
    ##### 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: Your input file has a malformed header: We never saw the required CHROM header line (starting with one #) for the input VCF file

I guess it is because the VCF file is missing the standard header (#CHROM POS ID ...). I was wondering if there was any way I could add the header line to the VCF file to be able to use in GATK? Would it also be possible to use the original VCF with all variations in the Indel Realignment steps or would you recommend I subset it for insertions/deletions as I did in this case?

Thank you very much for your help and suggestions!!!

Yours sincerely Shanker Swaminathan

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 6,192Administrator, GATK Developer admin

    You can add in the line manually, but be sure to check there's nothing else wrong with your file. And you can use the original VCF, no need to subset.

    Geraldine Van der Auwera, PhD

  • sswaminathansswaminathan USAPosts: 16Member

    Hi Geraldine Thanks for the reply! Is there an easy way to add the line manually on the command line without having to actually open the file and insert the line? Could you also please explain as what I need to check to make the file is ok? Thanks Shanker

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 6,192Administrator, GATK Developer admin

    Hi Shanker,

    You can do that with some Unix magic, but not through the GATK, so I recommend you look up how to do it somewhere like StackOverflow.com.

    To validate your VCF you can use ValidateVariants.

    Geraldine Van der Auwera, PhD

  • sswaminathansswaminathan USAPosts: 16Member

    Hi Geraldine Thank you for the suggestion! I created a file with the header information and using the command sed '12r header.txt' Canis_familiaris.vcf > Canis_familiaris_header.vcf was able to add the header information. I then attempted to run ValidateVariants using the following command -R canFam3.1_bwa075/Canis_familiaris.CanFam3.1.73.dna.toplevel.fa -T ValidateVariants --variant Canis_familiaris_header.vcf --warnOnErrors but received the error

    ERROR MESSAGE: Input files canFam3.1_bwa075/TEST/Canis_familiaris_header.vcf and reference have incompatible contigs: Relative ordering of overlapping contigs differs, which is unsafe.

    I am not sure why I got the error message since I downloaded the reference and variation files from Ensembl release 73. I then tried to reorder the VCF using a perl program vcfsorter.pl (https://code.google.com/p/vcfsorter/) and reran the above command on the reordered VCF. I received the warnings

    WARN  10:54:12,527 ValidateVariants - ***** the REF allele is incorrect for the record at position 1:8113, fasta says A vs. VCF says T *****
    WARN  10:54:12,528 ValidateVariants - ***** the REF allele is incorrect for the record at position 1:8290, fasta says C vs. VCF says G *****
    WARN  10:54:12,529 ValidateVariants - ***** the REF allele is incorrect for the record at position 1:10374, fasta says T vs. VCF says A *****
    ....
    Found 36208 records with failures.

    I don't understand why there might be many mismatches. I downloaded the reference and variation files from Ensembl release 73. Could you please advise if I am doing anything incorrectly? How would you suggest I proceed?

    Thanks again for the help and suggestions!

    Shanker

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 6,192Administrator, GATK Developer admin

    Hi Shanker,

    It's difficult to comment on this since we don't work with those genomes ourselves. I know that in humans there are several different versions of builds of the reference genome; perhaps there is the same for dogs, and somehow you got mismatching files? I would recommend checking with Ensembl.

    Geraldine Van der Auwera, PhD

  • sswaminathansswaminathan USAPosts: 16Member

    Hi Geraldine

    Thank you for the suggestion. I checked with Ensembl and it appears that the variation file contained variants that should not have been present in the file as they had failed QC. They give me a list of variants to exclude and I removed those variants from the file. On rerunning the ValidateVariants tool, it came back with no failures.

    I then performed the Indel Realignment steps. When I tried to do the first step of Base Recalibration I received the following error:

    ERROR ------------------------------------------------------------------------------------------
    ERROR stack trace

    java.lang.ArrayIndexOutOfBoundsException at java.util.Arrays.copyOfRange(Unknown Source) at org.broadinstitute.sting.gatk.walkers.bqsr.BaseRecalibrator.calculateIsSNP(BaseRecalibrator.java:335) at org.broadinstitute.sting.gatk.walkers.bqsr.BaseRecalibrator.map(BaseRecalibrator.java:253) at org.broadinstitute.sting.gatk.walkers.bqsr.BaseRecalibrator.map(BaseRecalibrator.java:132) at org.broadinstitute.sting.gatk.traversals.TraverseReadsNano$TraverseReadsMap.apply(TraverseReadsNano.java:228) at org.broadinstitute.sting.gatk.traversals.TraverseReadsNano$TraverseReadsMap.apply(TraverseReadsNano.java:216) at org.broadinstitute.sting.utils.nanoScheduler.NanoScheduler.executeSingleThreaded(NanoScheduler.java:274) at org.broadinstitute.sting.utils.nanoScheduler.NanoScheduler.execute(NanoScheduler.java:245) at org.broadinstitute.sting.gatk.traversals.TraverseReadsNano.traverse(TraverseReadsNano.java:102) at org.broadinstitute.sting.gatk.traversals.TraverseReadsNano.traverse(TraverseReadsNano.java:56) at org.broadinstitute.sting.gatk.executive.LinearMicroScheduler.execute(LinearMicroScheduler.java:108) at org.broadinstitute.sting.gatk.GenomeAnalysisEngine.execute(GenomeAnalysisEngine.java:313) at org.broadinstitute.sting.gatk.CommandLineExecutable.execute(CommandLineExecutable.java:113) at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:245) at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:152) at org.broadinstitute.sting.gatk.CommandLineGATK.main(CommandLineGATK.java:91)

    ERROR ------------------------------------------------------------------------------------------
    ERROR A GATK RUNTIME ERROR has occurred (version 2.7-2-g6bda569):
    ERROR
    ERROR This might be a bug. Please check the documentation guide to see if this is a known problem.
    ERROR If not, please post the error message, with stack trace, to the GATK forum.
    ERROR Visit our website and forum for extensive documentation and answers to
    ERROR commonly asked questions http://www.broadinstitute.org/gatk
    ERROR
    ERROR MESSAGE: Code exception (see stack trace for error itself)
    ERROR ------------------------------------------------------------------------------------------

    The command used for the base recalibration is:

    Program Args: -T BaseRecalibrator -I AF23.rg.prmdup.jir.bam -R canFam3.1_bwa075/Canis_familiaris.CanFam3.1.73.dna.toplevel.fa -knownSites canFam3.1_bwa075/Canis_familiaris.CanFam3.1.73.variations.sorted.clean.recode.vcf -o AF23.rg.prmdup.recal_data.grp

    Could you please tell me what could be causing this error and how to correct it?

    Thank you once again from your help!

    Shanker

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 6,192Administrator, GATK Developer admin

    Hi Shanker,

    I'm glad your original problem was resolved. Regarding your new problem, can you please validate your BAM file (with Picard's ValidateSAMFile)? If the file is valid, then we may have a bug; if so please post it as a new discussion thread since it is unrelated to the original problem that is the topic of this thread. Thanks!

    Geraldine Van der Auwera, PhD

  • sswaminathansswaminathan USAPosts: 16Member

    Hi Geraldine

    I ran Picard ValidateSamFile on the file which had been Indel Realigned as follows:

    [Fri Oct 11 07:50:31 MST 2013] net.sf.picard.sam.ValidateSamFile INPUT=AF23.rg.prmdup.jir.bam OUTPUT=AF23.rg.prmdup.jir.bam.validate.summary MODE=SUMMARY REFERENCE_SEQUENCE=canFam3.1_bwa075/Canis_familiaris.CanFam3.1.73.dna.toplevel.fa    MAX_OUTPUT=100 IGNORE_WARNINGS=false VALIDATE_INDEX=true IS_BISULFITE_SEQUENCED=false MAX_OPEN_TEMP_FILES=8000 VERBOSITY=
    INFO QUIET=false VALIDATION_STRINGENCY=STRICT COMPRESSION_LEVEL=5 MAX_RECORDS_IN_RAM=500000 CREATE_INDEX=false CREATE_MD5_FILE=false
    [Fri Oct 11 07:50:31 MST 2013] Executing as sswaminathan@pnap-pe5-s32 on Linux 2.6.32-220.el6.x86_64 amd64; Java HotSpot(TM) 64-Bit Server VM 1.7.0_25-b15; Picard version: 1.99(1563)
    INFO    2013-10-11 07:53:14 SamFileValidator    Validated Read    10,000,000 records.  Elapsed time: 00:02:42s.  Time for last 10,000,000:  162s.  Last read position: 1:115,340,356
    INFO    2013-10-11 07:55:23 SamFileValidator    Validated Read    20,000,000 records.  Elapsed time: 00:04:50s.  Time for last 10,000,000:  128s.  Last read position: 11:44,170,416
    INFO    2013-10-11 07:57:27 SamFileValidator    Validated Read    30,000,000 records.  Elapsed time: 00:06:54s.  Time for last 10,000,000:  124s.  Last read position: 13:5,582,802
    INFO    2013-10-11 07:59:23 SamFileValidator    Validated Read    40,000,000 records.  Elapsed time: 00:08:50s.  Time for last 10,000,000:  115s.  Last read position: 14:48,332,184
    INFO    2013-10-11 08:01:21 SamFileValidator    Validated Read    50,000,000 records.  Elapsed time: 00:10:49s.  Time for last 10,000,000:  118s.  Last read position: 16:37,437,010
    INFO    2013-10-11 08:03:17 SamFileValidator    Validated Read    60,000,000 records.  Elapsed time: 00:12:44s.  Time for last 10,000,000:  115s.  Last read position: 18:28,772,673
    INFO    2013-10-11 08:05:12 SamFileValidator    Validated Read    70,000,000 records.  Elapsed time: 00:14:39s.  Time for last 10,000,000:  114s.  Last read position: 2:30,977,793
    INFO    2013-10-11 08:07:06 SamFileValidator    Validated Read    80,000,000 records.  Elapsed time: 00:16:33s.  Time for last 10,000,000:  114s.  Last read position: 21:19,182,389
    INFO    2013-10-11 08:08:54 SamFileValidator    Validated Read    90,000,000 records.  Elapsed time: 00:18:21s.  Time for last 10,000,000:  107s.  Last read position: 23:14,734,512
    INFO    2013-10-11 08:10:37 SamFileValidator    Validated Read   100,000,000 records.  Elapsed time: 00:20:04s.  Time for last 10,000,000:  103s.  Last read position: 25:32,544,596
    INFO    2013-10-11 08:12:24 SamFileValidator    Validated Read   110,000,000 records.  Elapsed time: 00:21:51s.  Time for last 10,000,000:  107s.  Last read position: 28:25,354,258
    INFO    2013-10-11 08:14:09 SamFileValidator    Validated Read   120,000,000 records.  Elapsed time: 00:23:37s.  Time for last 10,000,000:  105s.  Last read position: 3:50,890,865
    INFO    2013-10-11 08:15:48 SamFileValidator    Validated Read   130,000,000 records.  Elapsed time: 00:25:16s.  Time for last 10,000,000:   98s.  Last read position: 31:33,006,197
    INFO    2013-10-11 08:17:31 SamFileValidator    Validated Read   140,000,000 records.  Elapsed time: 00:26:58s.  Time for last 10,000,000:  102s.  Last read position: 34:30,930,685
    INFO    2013-10-11 08:19:09 SamFileValidator    Validated Read   150,000,000 records.  Elapsed time: 00:28:36s.  Time for last 10,000,000:   98s.  Last read position: 38:9,783,344
    INFO    2013-10-11 08:20:52 SamFileValidator    Validated Read   160,000,000 records.  Elapsed time: 00:30:19s.  Time for last 10,000,000:  102s.  Last read position: 5:13,330,940
    INFO    2013-10-11 08:22:24 SamFileValidator    Validated Read   170,000,000 records.  Elapsed time: 00:31:51s.  Time for last 10,000,000:   91s.  Last read position: 6:61,954,905
    INFO    2013-10-11 08:24:06 SamFileValidator    Validated Read   180,000,000 records.  Elapsed time: 00:33:33s.  Time for last 10,000,000:  102s.  Last read position: 8:15,694,274
    INFO    2013-10-11 08:25:33 SamFileValidator    Validated Read   190,000,000 records.  Elapsed time: 00:35:00s.  Time for last 10,000,000:   87s.  Last read position: MT:6,231
    INFO    2013-10-11 08:27:30 SamFileValidator    Validated Read   200,000,000 records.  Elapsed time: 00:36:57s.  Time for last 10,000,000:  116s.  Last read position: JH373500.1:46,352
    [Fri Oct 11 08:31:19 MST 2013] net.sf.picard.sam.ValidateSamFile done. Elapsed time: 40.79 minutes.
    Runtime.totalMemory()=776273920
    To get help, see http://picard.sourceforge.net/index.shtml#GettingHelp

    Here is the output from the summary file:

    ## HISTOGRAM   java.lang.String
    Error Type  Count
    ERROR:INVALID_CIGAR 497409
    ERROR:MATES_ARE_SAME_END    4316
    ERROR:MISMATCH_FLAG_MATE_NEG_STRAND 6120
    ERROR:MISMATCH_FLAG_MATE_UNMAPPED   78510
    ERROR:MISMATCH_MATE_ALIGNMENT_START 10

    Do you think there is a problem with a BAM file? What would you suggest we do next?

    I can start this post on another thread if you would like, I just wanted to make sure that this may be a different problem.

    Thanks

    Shanker

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 6,192Administrator, GATK Developer admin

    This line

    ERROR:INVALID_CIGAR 497409

    is alarming. It means that many of your reads have a malformed CIGAR string. You need to either fix them or run with -rf BadCigar so that GATK can skip them. If you do that, make sure to check the read filter summary report (in the console output) to see what proportion of data is getting excluded.

    Geraldine Van der Auwera, PhD

  • sswaminathansswaminathan USAPosts: 16Member

    Hi Geraldine

    Could you please tell me what could be causing this problem and how I should go about trying to fix it? Do I have to realign the data again?

    Thanks Shanker

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 6,192Administrator, GATK Developer admin

    This is typically an issue caused by the aligner (BWA etc). It is not necessarily a showstopper, if the proportion of reads affected is not too big.

    Beyond that I can't help you any further since this is not a GATK-specific problem; please get help either from the aligner developers or ask on a more general forum like SeqAnswers.com. Good luck!

    Geraldine Van der Auwera, PhD

  • sswaminathansswaminathan USAPosts: 16Member

    Thanks, Geraldine! I will take a further look and see if we can fix the problem.

Sign In or Register to comment.