HTSeq_output_Error_occured_when_processing_SAM_input

[Olha]

Dear, SEQanswers community,

I faced issue with creation counting tables for RNA-seq data using HTSeq software.
Before I performed TopHat alignment (accepted_hints.bam) for 28 fastq files (single-end reads) and Cuffmerge assemble (merged.gtf).

This is script, that I used:
htseq-count --format bam --stranded no accepted_hits.bam /merged.gtf > htseq_out.txt

At the end of HTSeq output file I got a message:

33300000 SAM alignment records processed.
Error occured when processing SAM input (record #33301227 in file /media/olha/3CEBA24E475DF793/ALL_olha_2/A15/tophat_out/accepted_hits.bam):
reference_id -1 out of range 0<=tid<84
[Exception type: ValueError, raised in calignmentfile.pyx:642]

Any suggestions will be highly appreciated.

Olha

[dpryan]

A "reference_id -1" should only occur if a read is unmapped. My guess is that you somehow have an unmapped read without the unmapped bit in the flag being set. What's the output of:

samtools view -F 4 /media/olha/3CEBA24E475DF793/ALL_olha_2/A15/tophat_out/accepted_hits.bam | awk '{if($3=="*") print $0}'

That should produce nothing, though I suspect it'll print the problematic line in the BAM file.

[Olha]

Dear, Devon,

Thank you for reply.
I run your script and get this output:

HWI-ST538:357:D2BKUACXX:1:1105:13318:13823	16	*	1204092	50	7M600N25M	*	0	0	CGCACGGACGCCCCCAAAACGCATATGACTCG	HE=AJJJJJJJJIGIJIGHHHHHHFFFFF@CB	AS:i:0	XM:i:2	XO:i:0	XG:i:0	MD:Z:17G13A0	NM:i:2	XS:A:+	NH:i:1
HWI-ST538:357:D2BKUACXX:1:1107:19710:10717	0	*	3860981	50	2M516N22M	*	0	0	CCGCCACTCCACGATGATGGGGTT	<?@DDFFDFFAFHHII?FHGGHIG	AS:i:0	XM:i:2	XO:i:0	XG:i:0	MD:Z:13G8C1	NM:i:2	XS:A:-	NH:i:1
HWI-ST538:357:D2BKUACXX:1:2314:13745:61117	0	*	4385042	50	12M1572N11M	*	0	0	TGGGGTGGGGTCCTGGCTTTGTG	CCCFFBDFHHHHHJJJJJJJJHI	AS:i:0	XM:i:2	XO:i:0	XG:i:0	MD:Z:18C2G1	NM:i:2	XS:A:-	NH:i:1

How can I extract these reads from my accepted_hints.bam file?

Thank You!

[Olha]

Also I tried this script for removing aforemention reads. I put them all into txt file (read_ids_to_remove.txt).

samtools view -h accepted_hits.bam [grep -vf read_ids_to_remove.txt] samtools view -bS -o accepted_hits_filter.bam

But I got error message:

view: invalid option -- 'v'
[samopen] no @SQ lines in the header.
[main_samview] random alignment retrieval only works for indexed BAM files.

Does it means, that I need to sorted and indexing accepted_hits.bam file before removing unmapped reads?