# Counting the reads that map to each gene¶

Now that we know which reads go with which gene, we’ll use htseq-count.

First, load the PySAM and HTSeq software packages:

module load PySAM/0.6


And next, run HTSeq:

htseq-count --format=bam --stranded=no --order=pos \
tophat_salivary_repl1/accepted_hits.bam \
~/RNAseq-model/Homo_sapiens/Ensembl/GRCh37/Annotation/Genes/genes.gtf > salivary_repl1_counts.txt


When this is done, type:

less salivary_repl1_counts.txt


(again, use ‘q’ to exit). These are your gene counts.

Note, these are raw gene counts - the number of reads that map to each feature (gene, in this case). They are not normalized by length of gene. According to this post on seqanswers, both DEseq and edgeR want exactly this kind of information!

Questions:

• what are the ‘ENS*...’ names?
• what do these parameters mean?
• what parameters does HTSeq take?
• why are we using so many programs?