This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Getting different results from bedtools and samtools

Hi, I am trying to extract expression counts from some bam files. When I used bedtools to do so, I got all 0 for every specified region (in bed file). But when I used samtools, the counts were there.

For bedtools, the code is following:

bedtools multicov -bams in.bam -bed annotation.bed | head -3

0   1   17578815    17609595    LAP3    0
1   1   201980827   202041410   CFLAR   0
2   1   45429998    45590913    LARS2   0

For samtools, I used

samtools view -c in.bam chr1:17578815-17609595

and the result I got is 24, and for the other two regions, I got 11 and 5147 respectively.

So I am wondering that what I was doing wrong? Or which result should I trust?

bam rna-seq bedtools samtools
• 2,214 views • link
updated by Ram Ram • written by wstla27 • wstla27

Check that your chromosome names in the bed file have the same name like in your bam file, e.g. chr1 is different to 1.

• 0 views • link
written by finswimmer finswimmer

I tried adding chr to the chromosome number in the bed file, but I think bedtools can not recognize such form. It is giving me the error: Unexpected file format. Please use tab-delimited BED, GFF, or VCF. Perhaps you have non-integer starts or ends at line 1?

• 0 views • link
written by wstla27 • wstla27

Please show us some lines of your bed file. Also make sure that the columns are tab delimited.

• 0 views • link
written by finswimmer finswimmer

1 answer

because multicov has some default filters: https://github.com/arq5x/bedtools2/blob/master/src/multiBamCov/multiBamCovMain.cpp#L41

bool keepDuplicates    = false;
bool keepFailedQC = false;

see also option -F and -f of samtools view

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Log in to answer this question.

More posts like this