The current GATK version is 3.7-0
Examples: Monday, today, last week, Mar 26, 3/26/04

Howdy, Stranger!

It looks like you're new here. If you want to get involved, click one of these buttons!

Did you remember to?


1. Search using the upper-right search box, e.g. using the error message.
2. Try the latest version of tools.
3. Include tool and Java versions.
4. Tell us whether you are following GATK Best Practices.
5. Include relevant details, e.g. platform, DNA- or RNA-Seq, WES (+capture kit) or WGS (PCR-free or PCR+), paired- or single-end, read length, expected average coverage, somatic data, etc.
6. For tool errors, include the error stacktrace as well as the exact command.
7. For format issues, include the result of running ValidateSamFile for BAMs or ValidateVariants for VCFs.
8. For weird results, include an illustrative example, e.g. attach IGV screenshots according to Article#5484.
9. For a seeming variant that is uncalled, include results of following Article#1235.

Did we ask for a bug report?


Then follow instructions in Article#1894.

Formatting tip!


Surround blocks of code, error messages and BAM/VCF snippets--especially content with hashes (#)--with lines with three backticks ( ``` ) each to make a code block.
Powered by Vanilla. Made with Bootstrap.
Picard 2.9.0 is now available. Download and read release notes here.
GATK 3.7 is here! Be sure to read the Version Highlights and optionally the full Release Notes.

GATK IndelRealigner removing too many reads

by0by0 UKMember Posts: 3

Hello,

I'm quite new to SNP calling. I am trying to setup a pipeline which includes GATK IndelRealigner as a final step. My bam file (before realignment) is a little over 1GB. After running the indel realigner however, it's reduced to 18MB! I'm assuming its throwing out way too many reads or something has gone wrong.

I'm calling the indel realigner with the default options as follows:

java -Xmx16g -jar $GATK_DIR/GenomeAnalysisTK.jar \
   -T IndelRealigner \
   -R /path/to/my/ref \
   -I input.bam.intervals \
   -targetIntervals input.bam.intervals \
   -o realn.bam \

I am generating the read groups using AddOrReplaceReadGroups.jar (from picard tools) and interval file using GATK RealignerTargetCreator with default options.

My bam file was generated off the raw reads of experiment SRA181417 fetched from SRA (after cleaning adapters using cutadapt, mapping to reference using bwa-mem, and removing duplicate reads using picard tools)

I have tried this on other reads and do not have the same issue.
Can anyone comment on why indel realigner could be throwing out so many reads.

Thank you

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Administrator, Dev Posts: 11,184 admin

    Can you post the full console log output?

    Geraldine Van der Auwera, PhD

  • by0by0 UKMember Posts: 3

    Ah. I just realised there was an error in the output, with the error message being: We encountered a non-standard non-IUPAC base in the provided reference: '10'
    I would think SRA reads shouldn't contain non-IUPAC bases.
    Any idea?

    Thanks

    INFO  14:41:02,980 HelpFormatter - The Genome Analysis Toolkit (GATK) v3.3-0-g37228af, Compiled 2014/10/24 01:07:22 
    INFO  14:41:02,981 HelpFormatter - Copyright (c) 2010 The Broad Institute 
    INFO  14:41:02,981 HelpFormatter - For support and documentation go to http://www.broadinstitute.org/gatk 
    INFO  14:41:02,985 HelpFormatter - Program Args: -T IndelRealigner -R /reference/ref.fa -I ./aln3_readgroups.bam -targetIntervals ./aln3_readgroups.bam.intervals -o ./aln4_realn.bam 
    INFO  14:41:02,990 HelpFormatter - Executing on Linux 2.6.32-504.1.3.el6.x86_64 amd64; Java HotSpot(TM) 64-Bit Server VM 1.7.0_76-b13. 
    INFO  14:41:02,990 HelpFormatter - Date/Time: 2015/05/19 14:41:02 
    INFO  14:41:02,990 HelpFormatter - -------------------------------------------------------------------------------- 
    INFO  14:41:02,991 HelpFormatter - -------------------------------------------------------------------------------- 
    INFO  14:41:03,540 GenomeAnalysisEngine - Strictness is SILENT 
    INFO  14:41:03,596 GenomeAnalysisEngine - Downsampling Settings: No downsampling 
    INFO  14:41:03,604 SAMDataSource$SAMReaders - Initializing SAMRecords in serial 
    INFO  14:41:03,617 SAMDataSource$SAMReaders - Done initializing BAM readers: total time 0.01 
    INFO  14:41:03,688 GenomeAnalysisEngine - Preparing for traversal over 1 BAM files 
    INFO  14:41:03,693 GenomeAnalysisEngine - Done preparing for traversal 
    INFO  14:41:03,693 ProgressMeter - [INITIALIZATION COMPLETE; STARTING PROCESSING] 
    INFO  14:41:03,693 ProgressMeter -                 | processed |    time |    per 1M |           |   total | remaining 
    INFO  14:41:03,694 ProgressMeter -        Location |     reads | elapsed |     reads | completed | runtime |   runtime 
    INFO  14:41:03,842 ReadShardBalancer$1 - Loading BAM index data 
    INFO  14:41:03,843 ReadShardBalancer$1 - Done loading BAM index data 
    INFO  14:41:18,190 GATKRunReport - Uploaded run statistics report to AWS S3 
    ##### ERROR ------------------------------------------------------------------------------------------
    ##### ERROR A USER ERROR has occurred (version 3.3-0-g37228af): 
    ##### 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: Bad input: We encountered a non-standard non-IUPAC base in the provided reference: '10'
    ##### ERROR ------------------------------------------------------------------------------------------
    [W::bam_hdr_read] EOF marker is absent. The input is probably truncated.
    
    
  • Geraldine_VdAuweraGeraldine_VdAuwera Administrator, Dev Posts: 11,184 admin

    Ah, this is a classic. Your reads are fine, but your reference file is malformed. It was probably stored/edited on a Windows machine at some point, and a Windows-specific character was introduced (the '10' that the error refers to is Windows encoding for new line iirc). On a Unix platform this is then interpreted incorrectly and so the tool fails. Just fix the file encoding or re-download an original copy if the reference and you'll be fine.

    Geraldine Van der Auwera, PhD

  • by0by0 UKMember Posts: 3
    edited May 2015

    Thank you for your reply Geraldine.

    That's odd. My reference / all indexing was done on Linux (Red Hat). My chromosome fasta file comes from this file (at the end)

    Do you see anything wrong with it?

    Here is how I'm indexing.

    # Path to reference fasta (extracted from gff)
    FA_PATH=./reference/ref.fa
    
    # Generate bwa index files for use with samtools
    bwa index -a bwtsw $FA_PATH 
    
    # Generate .fai file for use with gatk 
    samtools faidx $FA_PATH
    
    # Create reference .dict for GATK
    java -jar $PICARD_DIR/CreateSequenceDictionary.jar R=$FA_PATH O=$DICT_PATH.dict GENOME_ASSEMBLY=$GENOME_ASSEMBLY SPECIES=$SPECIES
    
  • Geraldine_VdAuweraGeraldine_VdAuwera Administrator, Dev Posts: 11,184 admin

    Can't look at the file right now (on smartphone) but FYI we've also seen this happen when the sequence is all on one line or when there is an empty line in the sequence. Do a search (search box, top right corner of page) for the error message to see how people have dealt with this in the past.

    Geraldine Van der Auwera, PhD

Sign In or Register to comment.