This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Incompatibility between bbmap and htseq?

I would like to count rnaseq mapped reads on a bacterial genome. Firstly I aligned the reads using bbmap with the command: bbmap.sh ref=P_riograndensis_LN831776.fas in=A1_trim.fastq out=A1_bbmap.sam

However, when I try to count the mapped reads with htseq (htseq-count -i locus_tag -t gene A1_bbmap.sam P_riograndensis_LN831776.gff), the following error occurs: "("Illegal CIGAR string '48='", 'line 5 of file A1_bbmap.sam')"

Here is the line 5 of A1_bbmap.sam. ASFH0:01510:12713 0 LN831776 LN831776 Paenibacillus riograndensis SBR5 genome assembly SBR5(T) ,chromosome : I. 81663 3 48= * 0 0 TGAGAAGCCGGAGCCGCAGCGAAAGCGAGTCTGAATAGGGCGACTGAG ??>;;6;<6<7@@<7<><>=<>=4<=><=<===?7><;;4;;<<;;;; XT:A:R NM:i:0 AM:i:3

I believe something is wrong with bbmap sam output. But I could not realize what is the problem. Maybe the "=" signal is not expected?

bbmap htseq
• 2,929 views • link
updated by GenoMax GenoMax • written by fhsantanna fhsantanna

2 answers

@fhsantanna: You won't have to do your alignments again (since that could take time). Instead you can use reformat.sh from BBMap which can use the "sam=1.3" flag for converting "="/"X" symbols to "M" to reformat your existing BAM file. This was based on some discussion about this with @Brian late last year.
BTW: This "incompatibility" also applies to featureCounts (post #140) which does not support those tags. I am not sure if a new version has been released since.

• 0 views • link
written by GenoMax GenoMax

Hi fhsantanna,

you need to set the parameter sam=1.3 during mapping with bbmap in order to have the sam version which is needed for htseq. In sam v 1.4 you can distinguish in the cigar string between match and missmatch, in version 1.3 you only see "M" as a (miss)match.

Cheers, Michael

• 0 views • link
written by michael.ante michael.ante

Log in to answer this question.

More posts like this