This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Reordering Bam files based on chromosome number, when the chromosome format is chr1..etc

When trying to call variants I get the error message:

##### ERROR MESSAGE: Lexicographically sorted human genome sequence detected in reads. Please see https://software.broadinstitute.org/gatk/documentation/article?id=1328for more information. Error details: reads contigs = [chr1, chr10, chr11, chr12, chr13, chr14, chr15, chr16, chr17, chr18, chr19, chr2, chr20, chr21, chr22, chr3, chr4, chr5, chr6, chr7, chr8, chr9, chrM, chrX, chrY]

Both my fasta file and my bams have the chromosomes sorted in this order. Therefore I need to re-sort my bams and possibly my fasta to be in numerical order: chr1,chr2 etc.

However, due to the chr in the chromosome name samtools sort does not sort the chromosomes into the desired order. Can anyone help with this sorting?

rna-seq snp
• 4,648 views • link
updated by finswimmer finswimmer • written by jparker4 • jparker4

It would be useful to mention what program you are using to call variants.

• 1 views • link
written by GenoMax GenoMax

Apologies, that was misleading. It is part of a variant calling pipeline but the actual error occurs when trying to run SplitNCigarReads. The full command I am running is posted below the comment underneath this one.

• 1 views • link
written by jparker4 • jparker4

1 answer

Hello,

in the error message there is a link to the documentation about what is wrong and how you can solve this. 

Version-specific alert for GATK 3.5

In version 3.5, we added some beefed-up VCF sequence dictionary validation. Unfortunately, as a side effect of the additional checks, some users have experienced an error that starts with "ERROR MESSAGE: Lexicographically sorted human genome sequence detected in variant." that is due to unintentional activation of a check that is not necessary. This will be fixed in the next release; in the meantime -U ALLOW_SEQ_DICT_INCOMPATIBILITY can be used (with caution) to override the check.

fin swimmer

(Aargh, third edit)

• 1 views • link
written by finswimmer finswimmer

Hi, thanks for the reply. Unfortunately the error is occurring when I run SplitNCigarReads, which requires the -U option ALLOW_N_CIGAR_READS. The tool will not let me pass two -U options in, so I'm a bit stuck.

Below is the error causing command: 

java -Xmx10G -Djava.io.tmpdir=/fastdata/mbp15jdp/tmp -jar ~/Downloads/GenomeAnalysisTK-3.8-0-ge9d806836/GenomeAnalysisTK.jar                    -T SplitNCigarReads                     -R /shared/sudlab1/General/mirror/genomes/hisat/hg19.fa                    -I deduped.dir/GTEX-12WSH-Cerebellum.hisat.bam                     -o split.dir/GTEX-12WSH-Cerebellum.hisat.split.bam                     -rf ReassignOneMappingQuality                     -RMQF 255                     -RMQT 60                     -U ALLOW_N_CIGAR_READS
• 1 views • link
updated by GenoMax GenoMax • written by jparker4 • jparker4

Hello,

ok than we have to resort your contig (that was my first solution before I saw the above information) even if I have no idea why GATK cares about the current order.

1. Sort your reference

$ samtools faidx reference.fa
$ cut -f1 reference.fa.fai|sort -k1V|parallel -k 'samtools faidx reference.fa {}' > ref_sorted.fa

2. Create a dictionary for the sorted ref

$ java -jar picard.jar CreateSequenceDictionary R=ref_sorted.fa O=ref_sorted.dict

3. Reorder bam file

$ java -jar picard.jar ReorderSam I=original.bam O=reordered.bam R=ref_sorted.fa CREATE_INDEX=TRUE

fin swimmer

• 1 views • link
written by finswimmer finswimmer

Log in to answer this question.

More posts like this
  • How can I get this sorted by ascending order? (when the number and character is together)
    written by woojoy14 • woojoy14
    samtools view NA12878_dedup.bam | cut -f 3 | sort -u > chr.list8 > chr1 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr2 …
  • Removing a chromosome from a fasta file
    written by t86dan • t86dan
    Hello, I have a fasta file for human Genome GRCh37 Reference Assembly, but for some reason chromosome 20 is repeated and I would like to …
  • MAFTOOLS HELP - Mutational Signature
    written by David_emir David_emir
    Hi All, I am trying to run mutational signatures on my MAF file, using maftools (https://github.com/PoisonAlien/maftools). while trying to run this i am constantly My …
  • Error with correctGCBias in deepTools
    written by clmin • clmin
    Hello, I'm trying to correct the GC bias with deepTools. I successfully used computeGCBias > computeGCBias -b file.bam --effectiveGenomeSize 2150570000 -g mm9.2bit -l 200 --GCbiasFrequenciesFile …
  • diferent contigs chrM
    written by cristina_sabiers cristina_sabiers
    Hi! I try to call snps to learn how to do it. I got myt bam file, I hadn't made it. I have the snp, …
  • coverageBed with -g option ?
    written by bioguy24 bioguy24
    I am running a coverageBed and getting the below error. ``` cmccabe@DTV-A5211QLM:~/Desktop/NGS$ coverageBed -d -sorted -a /home/cmccabe/Desktop/NGS/bed/bedtools/xgen_targets_sorted.bed -b /home/cmccabe/Desktop/NGS/pool_I_090215/IonXpress_008_150902_newheader.bam > /home/cmccabe/Desktop/NGS/pool_I_090215/IonXpress_008_150902_output.txt ERROR: Database file /home/cmccabe/Desktop/NGS/pool_I_090215/IonXpress_008_150902_newheader.bam …
  • problem with methylKit of the annotation
    written by hua.peng1314 hua.peng1314
    I am dealing with RRBS data follow [the manual][1]. I met some problem when doing annotation.The problem showed below. ``` > read.transcript.features("mm9.refseq.genes.bed",remove.unsual=TRUE,up.flank=1500,down.flank=500,unique.prom=TRUE) > diffAnn <- …
  • CONTRA problem. It does not stop.
    written by mafonso • mafonso
    Hi, I am trying to run an algorithm for CNV detection and was testing CONTRA. The test file works but using my own files the …
  • How To Retrieve Chromosome Lengths Using Bsgenome Package?
    written by James Ashmore • James Ashmore
    <p>Hello, I am following along with the <a href='http://www.ebi.ac.uk/training/online/sites/ebi.ac.uk.training.online/files/user/42/documents/embo_chipseq_course_1101026.pdf'>EBI's online course for ChIP-Seq analysis</a> and have encountered an error when trying to retrieve chromosome lengths:</p> …
  • Error Import File With "Import" Of Rtracklayer
    written by ska007 • ska007
    <p>Hi, This is related to a the question found here (http://biostar.stackexchange.com/questions/3414?sort=newest#sort-top). I am trying to import a gff file using import("file.gff") which contains a set …