CIGAR error when running htseq-count after BBMap

[sunkid]

I am getting the following types of errors when running htseq-count on a SAM file generated with bbmap.sh

Error occured when reading beginning of SAM/BAM file.
  ("Illegal CIGAR string '66='", 'line 81726 of file bb_mapped.sam')
  [Exception type: ValueError, raised in _HTSeq.pyx:1116]

The above SAM file was generated with a command like

bbmap.sh build=2 in=stdin.fq bamscript=bam.sh maxindel=200000 ambiguous=toss \
 usejni=t outu=bb_unm apped.sam outm=bb_mapped.sam rpkm=rpkm.txt -Xmx28g

The htseq-count command I used was

htseq-count -i gene bb_mapped.sam genome.gff > ! counts.txt

Any pointers greatly appreciated!

[GenoMax]

See this for a solution: https://www.biostars.org/p/182156/ You won't need to re-do the alignment but you would have to reformat your bam files.

[sunkid]

See this for a solution: https://www.biostars.org/p/182156/ You won't need to re-do the alignment but you would have to reformat your bam files.

Perfect! Thank you. For some reason I thought 1.3 was the default and so I had only experimented with the 1.4 flag.

[GenoMax]

I assume reformat.sh command will work (can you confirm, if you try it?). I have been using sam=1.3 to avoid this issue.

[sunkid]

I assume reformat.sh command will work (can you confirm, if you try it?). I have been using sam=1.3 to avoid this issue.

Yes, reformat.sh does work.

[chiayi]

On a related note, I encountered this issue even after reformat.

The sam file was generated by BBMap,

bbmap.sh -Xmx32g in1=r1.BBDuk.fastq.gz in2=r2.BBDuk.fastq.gz maxindel=2000 outm=mapped.sam outu=unmapped.sam ehist=ehist

I used reformat to convert the CIGAR to samtools v1.3,

reformat.sh sam=1.3 in=in.sam out=reformat.sam

However, all the reads were still unassigned when I ran featureCounts,

featureCounts -T 4 -p -C -O -a ref.gtf -R -o output reformat.sam

Here's one example read from the BBMap-generated SAM,

HWI-ST975:91:C49TLACXX:7:1101:2847:2191 1:N:0:CAGATC    83      Chr2 CHROMOSOME dumped from ADB: Jun/20/09 14:54; last updated: 2009-02-02      17851283 44       101M    =       17851177        -207    TGCTGCAAATGCATTTGTTCCTCGTAAACGTCCCAATACGGCTGGTAGAGTTTCAGTGGAACATCCAAACGGTGAACATCCTTCAAGGACATTATTTGTTA   8:DDDECC@@>A=DDDDB=/@=;6BCFFEE?23CA@@@F@JIIEHFE3JJIGHF>HGGGGDCBIG?GEE?GIIIGIEJGGGIJGGJJJHFGGGFFFFFCBB     NM:i:0  AM:i:44
HWI-ST975:91:C49TLACXX:7:1101:2847:2191 1:N:0:CAGATC    163     Chr2 CHROMOSOME dumped from ADB: Jun/20/09 14:54; last updated: 2009-02-02      17851177 44       101M    =       17851283        207     TCGAAGAGTGTGATGTCTTTTGCACGGGAGGGGGTATGGAGTTGGATGTTGAATCCCAAGATAATCATGCTGTTGATGCGTCAGGGATGCAGATTTCTGAT   @;;DD;DDBDBFFEH>FGH>@HHHHHI<AGIIID7;@FHADGA4C=CEEHACFFFFFEECEC@;@@CCCCCCC@C:>>(88><BC<89??>>3>C>ADC:A     NM:i:0  AM:i:44

To makre sure this read was not mapped to the unannotated region, I used to same featureCounts command with a HiSat2-generated BAM and that same read was assigned.

HWI-ST975:91:C49TLACXX:7:1101:2847:2191 83      Chr2    17851283        60      101M    =       17851177        -207    TGCTGCAAATGCATTTGTTCCTCGTAAACGTCCCAATACGGCTGGTAGAGTTTCAGTGGAACATCCAAACGGTGAACATCCTTCAAGGACATTATTTGTTA     8:DDDECC@@>A=DDDDB=/@=;6BCFFEE?23CA@@@F@JIIEHFE3JJIGHF>HGGGGDCBIG?GEE?GIIIGIEJGGGIJGGJJJHFGGGFFFFFCBB     AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:101        YS:i:0  YT:Z:CP NH:i:1
HWI-ST975:91:C49TLACXX:7:1101:2847:2191 163     Chr2    17851177        60      101M    =       17851283        207     TCGAAGAGTGTGATGTCTTTTGCACGGGAGGGGGTATGGAGTTGGATGTTGAATCCCAAGATAATCATGCTGTTGATGCGTCAGGGATGCAGATTTCTGAT     @;;DD;DDBDBFFEH>FGH>@HHHHHI<AGIIID7;@FHADGA4C=CEEHACFFFFFEECEC@;@@CCCCCCC@C:>>(88><BC<89??>>3>C>ADC:A     AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:101        YS:i:0  YT:Z:CP NH:i:1

The CIGARs are identical (101M) in both SAM and BAM, so I wasn't sure what may be the cause of this issue. Has anyone encounted this before? Thank you for any advice and input.

[GenoMax]

@chiyai: What version of BBMap are you using?

[chiayi]

@GenoMax
BBMap version 36.62
featureCounts v1.5.1
Thank you.

[Michael.Ante]

Hi chiayi,

in your BBmap index, the naming of the chromosome is "Chr2 CHROMOSOME dumped from ADB: Jun/20/09 14:54; last updated: 2009-02-02", whilst in the Hisat alignment it's just "Chr2". Either you rename each alignment in the bbamp results, or redo the alignment with a reduced fasta header, or rename the chromosmes in your gtf accordingly.

Cheers,

Michael

[GenoMax]

@Michael: Good catch. Not cleaning up fasta headers comes back and bites each time :-)

@Chiayi: Don't leave any spaces in fasta headers since they always cause problems with one thing or other.

[chiayi]

@Michael: I used sed to clean up the chromosome name in each alignment in the SAM. featureCounts works fine after that. Thank you so much.
@GenoMax: Thanks for the suggestion.
BTW, does reformat.sh only clean fasta header, or also clean sam/bam like my case here?

[GenoMax]

BTW, does reformat.sh only clean fasta header, or also clean sam/bam like my case here?

reformat in your command iteration is not touching the headers it only converts CIGAR strings to SAM 1.3 format (bbmap uses SAM format 1.4 by default). featureCounts and HTSeq-count don't understand 1.4 format (as of last time I was looking at that). BTW: Samtools v.1.3.x is unrelated to issue of SAM format v.1.3.

[chiayi]

reformat in your command iteration is not touching the headers it only converts CIGAR strings to SAM 1.3 format (bbmap uses SAM format 1.4 by default). featureCounts and HTSeq-count don't understand 1.4 format (as of last time I was looking at that). BTW: Samtools v.1.3.x is unrelated to issue of SAM format v.1.3.

I was referring to the trd (trimreaddescription) option in reformat. I thought that would clean the fasta header after the first space, or did I misunderstand?

[GenoMax]

I was referring to the trd (trimreaddescription) option in reformat. I thought that would clean the fasta header after the first space, or did I misunderstand?

Discrepancy in your case was with the genome fasta headers and not read names. You could have done trd on fasta genome file to trim names after the first space and then made the index.

I suspect that option won't work when processing sam/bam files. Even if it did, it is for read names and not genome headers. @Brian will have to confirm.

[Brian Bushnell]

Yep - "trd" can be used on the fasta file prior to mapping, or the parameter "trd" can be added during mapping. However, once the reads are mapped, Reformat will not change the "rname" field of the sam file, just the "qname". That's a good idea, though - I'll change it so that it can do that.