# Running htseq-count over BAM files

I am trying to derive an expression matrix from BAM files using htseq-count. These are bulk RNASeq BAM's by the way.

I have read the htseq-count documentation as well as samtools and figured that the following command should work:

samtools sort -on file.bam | samtools view -h - | htseq-count --type=exon --idattr=gene_id --additional-attr=gene_name - Homo_sapiens.GRCh38.77.gtf > counts.txt


The .gtf file here is from ensembl.

I have tried many variations of the above command but have failed to generate a count matrix. The error message states:

"0 SAM alignments processed"

Here are the first few lines of my .gtf file:

#!genome-build GRCh37.p13
#!genome-version GRCh37
#!genome-date 2009-02
#!genome-build-accession NCBI:GCA_000001405.14
#!genebuild-last-updated 2013-09
1       pseudogene      gene    11869   14412   .       +       .       gene_id "ENSG00000223972"; gene_name "DDX11L1";
gene_source "ensembl_havana"; gene_biotype "pseudogene";
1       processed_transcript    transcript      11869   14409   .       +       .       gene_id "ENSG00000223972"; trans
cript_id "ENST00000456328"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; transcript_nam
e "DDX11L1-002"; transcript_source "havana";
1       processed_transcript    exon    11869   12227   .       +       .       gene_id "ENSG00000223972"; transcript_id
"ENST00000456328"; exon_number "1"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; trans
cript_name "DDX11L1-002"; transcript_source "havana"; exon_id "ENSE00002234944";


And all the genes in the text file have an expression value of 0. What exactly is going wrong in this step? Is the conversion to a SAM file the problem? I'm completely lost so any help would me much appreciated.

Thanks!

EDIT:: For anybody who sees this post cause they have a similar problem. Just note that the newer versions of htseq-count don't require sorted .bam files and, so following the editing of the .gtf file, all I needed to do was convert it to .sam using samtools view -h and then pipe this to htseq-count. Sorting the files prior to this conversion would cause errors.

• The obvious first thing to check is that your gtf and genome match. Did you get them from the same place? Specifically, do the chromosome names in the gtf match the ones in the bam? Mar 12 '19 at 16:37
• Good point..... I have just checked, and the gtf file contains no chromosomes actually. After running grep chr over the file it returns nothing. What is the reason for this? And how would I find one that works with the BAM files provided? Mar 12 '19 at 17:12
• "grep chr" won't find anything if the chromosome names are bare numbers. The question is, does your bam have bare numbers as chromosome names? Mar 12 '19 at 17:14
• Can you show us a few lines of your gtf file? I am guessing you have 1 instead of chr1. Mar 12 '19 at 17:18
• @h3ab74 you can change the names by simply running sed 's/^/chr/' file.gtf > file.fixed.gtf and then use the file.fixed.gtf for your analysis. However, and this is important: make sure that the bam files and the gtf are both referring to the same genome as Devon said, and also make sure there are no blank links in your gtf file. The command will only work correctly if all the lines in the file should start with chr. Mar 13 '19 at 20:14
• I ran a similar command to that one: sed /^[[:digit:]]/s/^/chr/ file.gtf > newfile.gtf. I checked out newfile.gtf and made sure the chromosome names were corrected; they're now chr1, chr2, etc. However, when I ran the original command with the new file I still got the same problem. Is it maybe an issue with my samtools view` command? Mar 14 '19 at 18:48