HTSeq not working with Bowtie2 .SAM

[all_your_base]

Hi all,

I am having a weird problem with my Bowtie2 .SAM output for use with HTseq to count reads that correspond to genes in a .gff file.

Usually, I can just feed my Bowtie1 .SAM into HTseq using the following command:

htseq-count -m union -s no -t gene -i ID -o myOutput.sam myInput.sam organism.gff


However, after switching to Bowtie2 and running the same command, I get gigabytes of this:


Warning: Read HWI-ST1234:350WK3ACXX:6:1101:1780:2126/1 claims to have an aligned mate which could not be found. (Is the SAM file properly sorted?)
Warning: Read HWI-ST1234:350WK3ACXX:6:1101:1780:2126/2 claims to have an aligned mate which could not be found. (Is the SAM file properly sorted?)
Warning: Read HWI-ST1234:350WK3ACXX:6:1101:1671:2238/1 claims to have an aligned mate which could not be found. (Is the SAM file properly sorted?)
Warning: Read HWI-ST1234:350WK3ACXX:6:1101:1671:2238/2 claims to have an aligned mate which could not be found. (Is the SAM file properly sorted?)
Warning: Read HWI-ST1234:350WK3ACXX:6:1101:2011:2134/1 claims to have an aligned mate which could not be found. (Is the SAM file properly sorted?)
Warning: Read HWI-ST1234:350WK3ACXX:6:1101:2011:2134/2 claims to have an aligned mate which could not be found. (Is the SAM file properly sorted?)


According to other forums, this usually happens when the SAM isn't sorted by read ID, so that htseq can't find the two halves of a paired-end read. However, I tried sorting my SAM in multiple ways, such as:

sort -k1 myfile.sam > myfile_sorted.sam


I still get the same error! Any help or suggestions are greatly appreciated

[dpryan]

If any of those reads are multimapped, then using the command line sort command will not do what you want. Use samtools sort -n

[all_your_base]

@dpryan,

Thanks for the reply, but can you please explain your answer? How does the samtools sort command differ than unix sort?

Also, since Bowtie2 produces a SAM file by default, to use SAMtools sort, do I have to first convert to BAM, then sort, then convert back to SAM?

Thanks...

[kmcarr]

The problem is the /1 and /2 in your read names. The SAM specification indicates that the names of paired reads be identical. SAM identifies read 1 or read 2 by the FLAG bits. Remove the /1 & /2 from the names in your SAM files and repeat your analysis.

[all_your_base]

@kmcarr

Wonderful, I trimmed the /1 and /2 off my reads and made sure the mates were next to each other after sorting, and HTSEQ runs fine without the previous error messages.

Quick question...
After processing a few thousands reads, HTSEQ reports the following error:

Warning: Malformed SAM line: MRNM != '*' although flag bit &0x0008 set
Warning: Malformed SAM line: RNAME != '*' although flag bit &0x0004 set

This is from raw Bowtie2 output; the only modifications were my /1 and /2 trimming and sorting.

Anyone have an idea where these errors are coming from??

Thanks!

[dpryan]

RNAME and MRNM are the name of the chromosome (or scaffold or whatever) to which the current read and its mate (for MRNM) map. Since the flags indicate that the reads are unmapped, it's just complaining that there's stuff here instead of an *, meaning "Not available". I don't recall ever seeing that with bowtie2, only bwa. You can normally ignore such warnings.

[all_your_base]

@dpryan,

Thanks for all your help. My analysis is working well now