htseq-count output for DESeq

[frymor]

Hi,

I am having a problem using the output of htseq-count with the DESeq package.

I am using htseq-count on a few files with this command:

for file in `ls -1 *.bam`
do
NEW_FILE=$(echo $file | sed "s/.bam//")
echo $NEW_FILE;
samtools index $file;
samtools view -h -o $NEW_FILE.sam $file;
htseq-count -s no $NEW_FILE.sam ~/genomes/celegans/Caenorhabditis_elegans/UCSC/ce10/Annotation/Genes/genes.gtf  > $NEW_FILE.htseq_count.txt;
echo "done with file $file \n"
done

When htseq-count is running I get the following warning:

/usr/local/lib/python2.6/dist-packages/HTSeq-0.5.3p3-py2.6-linux-x86_64.egg/HTSeq/__init__.py:8: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility
from _HTSeq import *
/usr/local/lib/python2.6/dist-packages/HTSeq-0.5.3p3-py2.6-linux-x86_64.egg/HTSeq/__init__.py:8: RuntimeWarning: numpy.flatiter size changed, may indicate binary incompatibility
from _HTSeq import *

But it runs all the way to the end.

Than i have a list of genes with read counts, BUT the first line has only one item.
This is how the output looks like:

	8060
21ur-1	0
21ur-10	0
21ur-100	0
21ur-1000	0
21ur-10000	0
21ur-10001	0
21ur-10002	0
21ur-10003	0
21ur-10004	0

I just can't figure out what this number is. It doesn't correspond to any parameters I checked (#genes with or w.o. reads, #genes in total, with no alignment etc.).

I can't even figure whether or not it has anything to do with the error massage I am getting.
Each of the Output files gets a number on the first line. They all seems to be random.

The problem is that DESeq accepts a two-columns data frame and when this number in the first line shows up, it return an error massage.

> cds_htseq <- newCountDataSetFromHTSeqCount(countTable, directory="HTSEQ_COUNT/")
Error in `colnames<-`(`*tmp*`, value = 1:27) : 
  attempt to set 'colnames' on an object with less than two dimensions

When I remove the first line in all the files, it works fine.


What is the reason for this output?

Thanks
Assa

[Simon Anders]

You may have a gene with no name in your GTF file, i.e., with empty gene ID.

Your error has little to do with the problem. It is about your numpy installation being incompatible to your HTSeq installation, maybe because one of the twoi is 64bit and the other 32 bit. This is ugly but as numpy is not used in htseq-count, it won't matter.