This is a test version of Biostars. For the public version, visit https://www.biostars.org.
coverage bed error but produces output file

The below command produces a file (674 MB).

coverageBed \
  -d \
  -sorted \
  -a /home/cmccabe/Desktop/NGS/bed/sorted_unix_5column_xgen_targets.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

Output

chr1    955542    955763    +    AGRN:exon.1    1    0
chr1    955542    955763    +    AGRN:exon.1    2    0
chr1    955542    955763    +    AGRN:exon.1    3    0
chr1    955542    955763    +    AGRN:exon.1    4    1

My question is even-though an output file results I get the below error:

cmccabe@DTV-A5211QLM:~/Desktop/NGS$ coverageBed -d -sorted -a /home/cmccabe/Desktop/NGS/bed/sorted_unix_5column_xgen_targets.bed -b /media/cmccabe/C2F8EFBFF8EFAFB9/pool_I_090215/IonXpress_008_150902.bam > /home/cmccabe/Desktop/NGS/pool_I_090215/IonXpress_008_150902_counts.txt
ERROR: chromomsome sort ordering for file /media/cmccabe/C2F8EFBFF8EFAFB9/pool_I_090215/IonXpress_008_150902.bam is inconsistent with other files. Record was:
chr10   68734   68755   X28LU:04418:11727   6   +

I checked the bam and it is coordinate sorted:

cmccabe@DTV-A5211QLM:~/Desktop/NGS/pool_I_090215$ samtools view -H IonXpress_008_150902_newheader.bam | grep SO
@HD    VN:1.4    GO:none    SO:coordinate

I just want to make sure that I can use the data output of if I am missing something. I can not find that record in the bam file, is that a sequencing error, a bug, or something unique to Ion Torrent bam files. Thank you :).

Sorted bed file

chr1    955542    955763    +    AGRN:exon.1
chr1    957570    957852    +    AGRN:exon.2
chr1    976034    976270    +    AGRN:exon.2;AGRN:exon.3;AGRN:exon.4

I also verified the sort order in the bam

cmccabe@DTV-A5211QLM:~/Desktop/NGS/pool_I_090215$ samtools view -H IonXpress_008_150902_newheader.bam | grep SQ | cut -f 2 | awk '{ sub(/^SN:/, ""); print;}'
chr1
chr2
chr3
chr4
chr5
chr6
chr7
chr8
chr9
chr10
chr11
chr12
chr13
chr14
chr15
chr16
chr17
chr18
chr19
chr20
chr21
chr22
chrX
chrY
chrM
ngs sort bedtools
• 6,112 views • link
updated by Ram Ram • written by bioguy24 bioguy24

Is it a lexical vs numerical sort thing: 1 10 .. 2 20 ..?

• 0 views • link
updated by Ram Ram • written by thackl thackl

You've double-checked the order of your bed file too, especially the chromosome ordering, and the naming convention? The fact that the "mis-ordering" occurs relatively early in chromosome 10 makes me think of the formal possibility that one of your files (maybe the bed file?) is sorted chr1 chr10 chr11 .. chr2 chr20 .. instead of chr1 chr2 .. chr9 chr10..

I know you know awk, but one could awk '{print $1}' aBed.bed | sort -u to quickly check

I'm guessing your output that you showed is not comprehensive, in which case you might check the tail-end of the output file to see if it also ends on chr1, or if it covered the entire genome.

• 0 views • link
updated by Ram Ram • written by Joseph Pearson Joseph Pearson

Looks like I just need to sort -k1,1 -k2,2n

awk '{print $1}' sorted_unix_5column_xgen_targets.bed | sort -u
chr1
chr10
chr11
chr12
chr13
chr14
chr15
chr16
chr17
chr18
chr19
chr2
chr20
chr21
chr22
chr3
chr4
chr5
chr6
chr7
chr8
chr9
chrX
chrY
• 0 views • link
updated by Ram Ram • written by bioguy24 bioguy24

2 answers

According to the bedtools manual -sorted needs to be done with sort -k1,1 -k2,2n, which would do a lexical sort and put chr10,11.. before chr2. You BAM file however, is sorted in a "human" order...

• 0 views • link
written by thackl thackl

-sorted just lets bedtools know that the input files are already sorted by some consistent rule. But you're right, that the bed file (sorted by the recommended -k1,1) would give a lexical sort.

Here's a post describing how to use Unix sort to sort in "human" order: How To Sort Bed Format File

• 0 views • link
updated by Ram Ram • written by Joseph Pearson Joseph Pearson

As the error says it look as though either your BAM or BED file are not sorted correctly. First inspect the sort order of the chromosomes in both the sorted BED and BAM files by running the following commands:

cut -f1 your_file.bed | uniq
samtools view your_file.bam | cut -f3 | uniq

If the order of the chromosomes is different then you will have to provide a genome file to Bedtools with the order your chromosomes are in the BED file.

• 0 views • link
updated by Ram Ram • written by James Ashmore James Ashmore

Thank you everybody :).

• 0 views • link
written by bioguy24 bioguy24

I mis-read the post and my bed was sorted correctly but my bam is in "human" sorting so I need to re-sort the bam in "lexical"? Thank you :).

The output file that results is (with the chromosome ordering):

cmccabe@DTV-A5211QLM:~/Desktop/NGS/bed$ awk '{print $1}' IonXpress_008_150902_counts.txt | sort -u
chr1
chr10
chr11
chr12
chr13
chr14
chr15
chr16
chr17
chr18
chr19
chr2
chr20
chr21
chr22
chr3
chr4
chr5
chr6
chr7
chr8
chr9
• 0 views • link
updated by Ram Ram • written by bioguy24 bioguy24

Yep, a simple samtools sort your.bam your.sorted should suffice.

• 0 views • link
updated by Ram Ram • written by thackl thackl

Beware samtools sort does not sort your alignments by chromosome, it sorts them by position in each chromosome. For instance, if your FASTA index which you aligned your reads to starts with chr12, chr13, chr15 e.t.c then your BAM file with start with chr12, chr13, chr15 e.t.c. Since it is much harder to sort the BAM file then the BED file, I suggest you sort your BED file by the order of the chromosomes in the FASTA index file.

• 0 views • link
updated by Ram Ram • written by James Ashmore James Ashmore

Log in to answer this question.

More posts like this
  • add gene name from bed file to vcf in vcftools
    written by bioguy24 bioguy24
    I am trying to add gene information to a vcf file using vcftools: vcftools annotate \ -a genes.bed.gz \ -c CHROM,FROM,TO,GENE \ -h <(echo '##INFO=<ID=GENE,Number=1,Type=String,Description="Gene …
  • picard error version 1.139 on ubuntu 14.04
    written by bioguy24 bioguy24
    <p>I seem to have trouble with calling picard from the command line but other tools work perfectly.<br /> <br /> path:</p> <pre> /home/cmccabe/Desktop/NGS/bedops-2.14.4/bin:/home/cmccabe/Desktop/NGS/samtools-1.2:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games:/usr/local/games:/home/cmccabe/Desktop/NGS/bedtools2-25.0/bin:/home/cmccabe/Desktop/NGS/bbmap-35.59:/home/cmccabe/Desktop/NGS/FastQC-0.11.4:/home/cmccabe/Desktop/NGS/picard-tools-1.139</pre> <p>The code …
  • 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 …
  • bedtools 2.25 vs bedtools2.20 coverageBed output not consistent
    written by bioguy24 bioguy24
    Hopefully, this edit is more clear in explaining the issue and will help led to a solution. As I am lost as to why there …
  • run picard on all bam files in directory
    written by bioguy24 bioguy24
    I am trying to run a bash loop on all bam files in a directory and having a bit of trouble. The bash so far …
  • samtools-1.2 re-header command on ubuntu 14.04.3
    written by bioguy24 bioguy24
    Why does samtools need to be in `/usr/bin` when it is in a directory (samtools-1.2) downloaded from github, unzipped, and make was run? ``` samtools …