genes counts using htseqcounts

Question

[bolbolman2000]

Hello All,

I am new to these RNA-seq thing, and I really need help.

I have .bam files that I need to get the genes counts from it to use in the edgeR. I don't have information about the reads but iknow that these .bam files is mapping for read of human. The .bam files belong to two group patient and normal.

I tried to use the htseq-count -f bam <.bam file > <.gtf file> > counts.txt

It run then it gives error. How I can know which .gtf file I should use? I know that the reference genome that has been used in mapping is hg19

Also if I have more than one .bam file for patient samples do I suppose to put all of them together in one file, or I can run on every one of them.

Please help

Answer 1

[GouthamAtla]

Ensemble GTF files works well with htseq. By default htseq-count looks for gene_id attribute in gff/gtf file. If you do not have gene_id attribute, you need to specify the attribute name.

Read faq section on htseq-count page. http://www-huber.embl.de/users/anders/HTSeq/doc/count.html

Otherwise post few lines of your gtf/gff file.

[bolbolman2000]

Thanks for you reply,

I did download the gtf file and I run the htseq count a gain using a sam file instead of the .bam after I used the samtools to convert it. but I got the following error

[Malformed SAM line: MRNM == '*' although flag bit &0x0008 cleared", 'line 68 of file sorted_genome_alignments.sam')
  [Exception type: ValueError, raised in _HTSeq.pyx:1323]

Do you have an idea what I can do to fix this?

[GouthamAtla]

HTSeq works pretty well without any complaints. Please make sure you have removed unaligned reads from your sam/bam file. Make sure you set correct option for sort type. The most important step is make sure you are using genome and annotation file from same source and same build.