3
$\begingroup$

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.

Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ 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? $\endgroup$
    – swbarnes2
    Mar 12, 2019 at 16:37
  • 1
    $\begingroup$ 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? $\endgroup$
    – h3ab74
    Mar 12, 2019 at 17:12
  • 2
    $\begingroup$ "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? $\endgroup$
    – swbarnes2
    Mar 12, 2019 at 17:14
  • 1
    $\begingroup$ @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. $\endgroup$
    – terdon ♦
    Mar 12, 2019 at 17:17
  • 1
    $\begingroup$ Can you show us a few lines of your gtf file? I am guessing you have 1 instead of chr1. $\endgroup$
    – terdon ♦
    Mar 12, 2019 at 17:18

1 Answer 1

Reset to default
2
$\begingroup$

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.

Improve this answer
$\endgroup$
5
  • $\begingroup$ 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. ? $\endgroup$
    – h3ab74
    Mar 13, 2019 at 16:29
  • 1
    $\begingroup$ 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. $\endgroup$
    – Devon Ryan
    Mar 13, 2019 at 17:00
  • $\begingroup$ @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. $\endgroup$
    – terdon ♦
    Mar 13, 2019 at 20:14
  • $\begingroup$ 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? $\endgroup$
    – h3ab74
    Mar 14, 2019 at 18:48
  • $\begingroup$ 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! $\endgroup$
    – h3ab74
    Mar 18, 2019 at 16:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.