Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-03-13T16:29:51.337Z
Answers.Text: 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. ?
Answers.UserDisplayName: h3ab74
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-03-13T17:00:57.657Z
Answers.Text: 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.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-03-13T20:14:12.477Z
Answers.Text: @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`.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-03-14T18:48:48.680Z
Answers.Text: 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?
Answers.UserDisplayName: h3ab74
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-03-18T16:37:14.293Z
Answers.Text: 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!
Answers.UserDisplayName: h3ab74
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-03-13T09:04:52.473Z
Answers.LastActivityDate: 2019-03-13T09:04:52.473Z
Answers.UserDisplayName: Devon Ryan
Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-03-12T16:37:19.763Z
Comments.Text: 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?
Comments.UserDisplayName: swbarnes2
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-03-12T17:12:40.990Z
Comments.Text: 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?
Comments.UserDisplayName: h3ab74
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-03-12T17:14:45.583Z
Comments.Text: "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?
Comments.UserDisplayName: swbarnes2
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-03-12T17:17:14.403Z
Comments.Text: @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.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-03-12T17:18:53.383Z
Comments.Text: Can you show us a few lines of your gtf file? I am guessing you have `1` instead of `chr1`.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-03-12T17:19:48.593Z
Comments.Text: My .bam file chromosomes are formatted as such: chr1, chr2,....
Comments.UserDisplayName: h3ab74
Tags: <rna-seq><sam><samtools><bash>
Title: Running htseq-count over BAM files
ThreadText: # 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**: @username_1 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**: @username_0 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 username_3 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 @username_3, good point, I made sure to use the same .gtf files!  
  
  ---
  
  