<p>If your BAM file is using UCSC chromosome names (chr1, chr2, etc.) and your GTF file is using Ensembl chromosome names (1, 2, etc.) then htseq-count won't know how to match your alignments to any of the genes. The simplest solution is to not use htseq-count, but instead to use featureCounts, which will handle the chromosome name difference transparently (it's one of the few programs that does) and is much faster anyway. In the future, download your GTF and fasta files from the same source, don't mix and match them.</p>

CC BY-SA 4.0
2019-03-13T16:29:51.337Z
Thanks for the reply! My BAM file does indeed use this chromosome name system.Thing is though, my lab wants to continue using htseq for consistency. Would this problem be fixed by editing the .gtf file chromosome names to chr1, chr2 etc. ?
h3ab74
CC BY-SA 4.0
2019-03-13T17:00:57.657Z
Yes, you can edit the GTF, though if you want real consistency you might see if you were using the exact same fasta and GTF files before.
Devon Ryan
CC BY-SA 4.0
2019-03-13T20:14:12.477Z
@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`.
terdon
CC BY-SA 4.0
2019-03-14T18:48:48.680Z
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?
h3ab74
CC BY-SA 4.0
2019-03-18T16:37:14.293Z
Thank you all for your helpful comments, I have made it work. And @DevonRyan, good point, I made sure to use the same .gtf files!
h3ab74
CC BY-SA 4.0
2019-03-13T09:04:52.473Z
2019-03-13T09:04:52.473Z
Devon Ryan
<p>I am trying to derive an expression matrix from BAM files using <code>htseq-count</code>. These are bulk RNASeq BAM's by the way.</p>
<p>I have read the <code>htseq-count</code> documentation as well as <code>samtools</code> and figured that the following command <em>should</em> work:</p>
<pre><code>samtools sort -on file.bam | samtools view -h - | htseq-count --type=exon --idattr=gene_id --additional-attr=gene_name - Homo_sapiens.GRCh38.77.gtf &gt; counts.txt
</code></pre>
<p>The <code>.gtf</code> file here is from ensembl.</p>
<p>I have tried many variations of the above command but have failed to generate a count matrix. The <em>error</em> message states:</p>
<p>&quot;0 SAM alignments processed&quot;</p>
<p>Here are the first few lines of my .gtf file:</p>
<pre><code>#!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 &quot;ENSG00000223972&quot;; gene_name &quot;DDX11L1&quot;;
gene_source &quot;ensembl_havana&quot;; gene_biotype &quot;pseudogene&quot;;
1       processed_transcript    transcript      11869   14409   .       +       .       gene_id &quot;ENSG00000223972&quot;; trans
cript_id &quot;ENST00000456328&quot;; gene_name &quot;DDX11L1&quot;; gene_source &quot;ensembl_havana&quot;; gene_biotype &quot;pseudogene&quot;; transcript_nam
e &quot;DDX11L1-002&quot;; transcript_source &quot;havana&quot;;
1       processed_transcript    exon    11869   12227   .       +       .       gene_id &quot;ENSG00000223972&quot;; transcript_id
 &quot;ENST00000456328&quot;; exon_number &quot;1&quot;; gene_name &quot;DDX11L1&quot;; gene_source &quot;ensembl_havana&quot;; gene_biotype &quot;pseudogene&quot;; trans
cript_name &quot;DDX11L1-002&quot;; transcript_source &quot;havana&quot;; exon_id &quot;ENSE00002234944&quot;;
</code></pre>
<p>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.</p>
<p>Thanks!</p>
<p>EDIT:: For anybody who sees this post cause they have a similar problem. Just note that the newer versions of <code>htseq-count</code> don't require sorted .bam files and, so following the editing of the <code>.gtf</code> file, all I needed to do was convert it to <code>.sam</code> using <code>samtools view -h</code> and then pipe this to <code>htseq-count</code>. Sorting the files prior to this conversion would cause errors.</p>

CC BY-SA 4.0
2019-03-12T16:37:19.763Z
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?
swbarnes2
CC BY-SA 4.0
2019-03-12T17:12:40.990Z
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?
h3ab74
CC BY-SA 4.0
2019-03-12T17:14:45.583Z
"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?
swbarnes2
CC BY-SA 4.0
2019-03-12T17:17:14.403Z
@swbarnes2 I converted your answer to a comment since it wasn't giving an answer but only requesting more information from the OP. Please only post answers if they are actually providing an _answer_. If you need more details, leave a comment instead.
terdon
CC BY-SA 4.0
2019-03-12T17:18:53.383Z
Can you show us a few lines of your gtf file? I am guessing you have `1` instead of `chr1`.
terdon
CC BY-SA 4.0
2019-03-12T17:19:48.593Z
My .bam file chromosomes are formatted as such: chr1, chr2,....
h3ab74
<rna-seq><sam><samtools><bash>
Running htseq-count over BAM files
# Running htseq-count over BAM files

(Asked by: username_0 on 2019-03-12)

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.

**Tags:** rna-seq, sam, samtools, bash
> **Comments:**
> **username_1**: 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?  
> **username_0**: 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?  
> **username_1**: "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?  
> **username_2**: @swbarnes2 I converted your answer to a comment since it wasn't giving an answer but only requesting more information from the OP. Please only post answers if they are actually providing an _answer_. If you need more details, leave a comment instead.  
> **username_2**: Can you show us a few lines of your gtf file? I am guessing you have `1` instead of `chr1`.  
> **username_0**: My .bam file chromosomes are formatted as such: chr1, chr2,....  

---

## 1 answers

### Answer by username_3 on 2019-03-13. Score: 2 | Accepted answer

If your BAM file is using UCSC chromosome names (chr1, chr2, etc.) and your GTF file is using Ensembl chromosome names (1, 2, etc.) then htseq-count won't know how to match your alignments to any of the genes. The simplest solution is to not use htseq-count, but instead to use featureCounts, which will handle the chromosome name difference transparently (it's one of the few programs that does) and is much faster anyway. In the future, download your GTF and fasta files from the same source, don't mix and match them.

> **Comments:**
> **username_0**: Thanks for the reply! My BAM file does indeed use this chromosome name system.Thing is though, my lab wants to continue using htseq for consistency. Would this problem be fixed by editing the .gtf file chromosome names to chr1, chr2 etc. ?  
> **username_3**: Yes, you can edit the GTF, though if you want real consistency you might see if you were using the exact same fasta and GTF files before.  
> **username_2**: @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`.  
> **username_0**: 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?  
> **username_0**: Thank you all for your helpful comments, I have made it work. And @DevonRyan, good point, I made sure to use the same .gtf files!  

---

