Counting Reads¶
The next step after mapping reads is to count the number of reads that fall within each annotated gene in the genome, so lets set up a count directory.
Shell Variables¶
[1]:
# Source the config script
source bioinf_intro_config.sh
mkdir -p $COUNT_OUT
ls $CUROUT
count_out igv qc_output stuff_for_igv_shorter_intron.tgz trimmed_fastqs
genome myinfo star_out stuff_for_igv.tgz
Counting Reads¶
[2]:
htseq-count --help
usage: htseq-count [options] alignment_file gff_file
This script takes one or more alignment files in SAM/BAM format and a feature
file in GFF format and calculates for each feature the number of reads mapping
to it. See http://htseq.readthedocs.io/en/master/count.html for details.
positional arguments:
samfilenames Path to the SAM/BAM files containing the mapped reads.
If '-' is selected, read from standard input
featuresfilename Path to the file containing the features
optional arguments:
-h, --help show this help message and exit
-f {sam,bam}, --format {sam,bam}
type of <alignment_file> data, either 'sam' or 'bam'
(default: sam)
-r {pos,name}, --order {pos,name}
'pos' or 'name'. Sorting order of <alignment_file>
(default: name). Paired-end sequencing data must be
sorted either by position or by read name, and the
sorting order must be specified. Ignored for single-
end data.
--max-reads-in-buffer MAX_BUFFER_SIZE
When <alignment_file> is paired end sorted by
position, allow only so many reads to stay in memory
until the mates are found (raising this number will
use more memory). Has no effect for single end or
paired end sorted by name
-s {yes,no,reverse}, --stranded {yes,no,reverse}
whether the data is from a strand-specific assay.
Specify 'yes', 'no', or 'reverse' (default: yes).
'reverse' means 'yes' with reversed strand
interpretation
-a MINAQUAL, --minaqual MINAQUAL
skip all reads with alignment quality lower than the
given minimum value (default: 10)
-t FEATURETYPE, --type FEATURETYPE
feature type (3rd column in GFF file) to be used, all
features of other type are ignored (default, suitable
for Ensembl GTF files: exon)
-i IDATTR, --idattr IDATTR
GFF attribute to be used as feature ID (default,
suitable for Ensembl GTF files: gene_id)
--additional-attr ADDITIONAL_ATTR
Additional feature attributes (default: none, suitable
for Ensembl GTF files: gene_name). Use multiple times
for each different attribute
-m {union,intersection-strict,intersection-nonempty}, --mode {union,intersection-strict,intersection-nonempty}
mode to handle reads overlapping more than one feature
(choices: union, intersection-strict, intersection-
nonempty; default: union)
--nonunique {none,all}
Whether to score reads that are not uniquely aligned
or ambiguously assigned to features
--secondary-alignments {score,ignore}
Whether to score secondary alignments (0x100 flag)
--supplementary-alignments {score,ignore}
Whether to score supplementary alignments (0x800 flag)
-o SAMOUTS, --samout SAMOUTS
write out all SAM alignment records into SAM files
(one per input file needed), annotating each line with
its feature assignment (as an optional field with tag
'XF')
-q, --quiet suppress progress report
Written by Simon Anders (sanders@fs.tum.de), European Molecular Biology
Laboratory (EMBL). (c) 2010. Released under the terms of the GNU General
Public License v3. Part of the 'HTSeq' framework, version 0.11.2.
We will use htseq-count to do the counting, but first we need to make some decisions, because the htseq-count defaults do not work with some annotation files. Here are the most important commandline options that we need to consider: * –format=: Format of the input data. Possible values are sam (for text SAM files) and bam (for binary BAM files). Default is sam. * –stranded=: whether the data is from a strand-specific
assay (default: yes). For stranded=no, a read is considered overlapping with a feature regardless of whether it is mapped to the same or the opposite strand as the feature. For stranded=yes and single-end reads, the read has to be mapped to the same strand as the feature. For paired-end reads, the first read has to be on the same strand and the second read on the opposite strand. For stranded=reverse, these rules are reversed. * –type=: feature type (3rd column in GFF file) to be evaluated,
all features of other type are ignored (default, suitable for RNA-Seq analysis using an Ensembl GTF file: exon) * –idattr=: GFF attribute to be used as feature ID. Several GFF lines with the same feature ID will be considered as parts of the same feature. The feature ID is used to identity the counts in the output table. The default, suitable for RNA-Seq analysis using an Ensembl GTF file, is gene_id.
And here is how we will set those options: * –format=bam: Since Tophat generated BAM files for us * –stranded=reverse: The dUTP method that we used for generating a strand-specific library produces reads that are anti-sense, htseq-count considers this to be “reverse”.
We need to look at the GFF file to understand what exactly the --type and --idattr options are, and why we are setting them this way.
[3]:
head -20 $GENOME_DIR/$GTF
#!genome-build CNA3
#!genome-version CNA3
#!genome-date 2015-11
#!genome-build-accession GCA_000149245.3
#!genebuild-last-updated 2015-11
1 ena gene 100 5645 . - . gene_id "CNAG_04548"; gene_source "ena"; gene_biotype "protein_coding";
1 ena transcript 100 5645 . - . gene_id "CNAG_04548"; transcript_id "AFR92135"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding";
1 ena exon 5494 5645 . - . gene_id "CNAG_04548"; transcript_id "AFR92135"; exon_number "1"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding"; exon_id "AFR92135-1";
1 ena CDS 5494 5645 . - 0 gene_id "CNAG_04548"; transcript_id "AFR92135"; exon_number "1"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding"; protein_id "AFR92135"; protein_version "1";
1 ena start_codon 5643 5645 . - 0 gene_id "CNAG_04548"; transcript_id "AFR92135"; exon_number "1"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding";
1 ena exon 5322 5422 . - . gene_id "CNAG_04548"; transcript_id "AFR92135"; exon_number "2"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding"; exon_id "AFR92135-2";
1 ena CDS 5322 5422 . - 1 gene_id "CNAG_04548"; transcript_id "AFR92135"; exon_number "2"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding"; protein_id "AFR92135"; protein_version "1";
1 ena exon 3958 5263 . - . gene_id "CNAG_04548"; transcript_id "AFR92135"; exon_number "3"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding"; exon_id "AFR92135-3";
1 ena CDS 3958 5263 . - 2 gene_id "CNAG_04548"; transcript_id "AFR92135"; exon_number "3"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding"; protein_id "AFR92135"; protein_version "1";
1 ena exon 3206 3890 . - . gene_id "CNAG_04548"; transcript_id "AFR92135"; exon_number "4"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding"; exon_id "AFR92135-4";
1 ena CDS 3206 3890 . - 1 gene_id "CNAG_04548"; transcript_id "AFR92135"; exon_number "4"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding"; protein_id "AFR92135"; protein_version "1";
1 ena exon 2846 3126 . - . gene_id "CNAG_04548"; transcript_id "AFR92135"; exon_number "5"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding"; exon_id "AFR92135-5";
1 ena CDS 2846 3126 . - 0 gene_id "CNAG_04548"; transcript_id "AFR92135"; exon_number "5"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding"; protein_id "AFR92135"; protein_version "1";
1 ena exon 2322 2782 . - . gene_id "CNAG_04548"; transcript_id "AFR92135"; exon_number "6"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding"; exon_id "AFR92135-6";
1 ena CDS 2322 2782 . - 1 gene_id "CNAG_04548"; transcript_id "AFR92135"; exon_number "6"; gene_source "ena"; gene_biotype "protein_coding"; transcript_source "ena"; transcript_biotype "protein_coding"; protein_id "AFR92135"; protein_version "1";
[4]:
head -50 $GENOME_DIR/$GTF | cut -c -55
#!genome-build CNA3
#!genome-version CNA3
#!genome-date 2015-11
#!genome-build-accession GCA_000149245.3
#!genebuild-last-updated 2015-11
1 ena gene 100 5645 . - . gene_id "CNAG_04548"; gene_so
1 ena transcript 100 5645 . - . gene_id "CNAG_04548"; t
1 ena exon 5494 5645 . - . gene_id "CNAG_04548"; transc
1 ena CDS 5494 5645 . - 0 gene_id "CNAG_04548"; transcr
1 ena start_codon 5643 5645 . - 0 gene_id "CNAG_04548";
1 ena exon 5322 5422 . - . gene_id "CNAG_04548"; transc
1 ena CDS 5322 5422 . - 1 gene_id "CNAG_04548"; transcr
1 ena exon 3958 5263 . - . gene_id "CNAG_04548"; transc
1 ena CDS 3958 5263 . - 2 gene_id "CNAG_04548"; transcr
1 ena exon 3206 3890 . - . gene_id "CNAG_04548"; transc
1 ena CDS 3206 3890 . - 1 gene_id "CNAG_04548"; transcr
1 ena exon 2846 3126 . - . gene_id "CNAG_04548"; transc
1 ena CDS 2846 3126 . - 0 gene_id "CNAG_04548"; transcr
1 ena exon 2322 2782 . - . gene_id "CNAG_04548"; transc
1 ena CDS 2322 2782 . - 1 gene_id "CNAG_04548"; transcr
1 ena exon 1823 2274 . - . gene_id "CNAG_04548"; transc
1 ena CDS 1823 2274 . - 2 gene_id "CNAG_04548"; transcr
1 ena exon 1556 1767 . - . gene_id "CNAG_04548"; transc
1 ena CDS 1556 1767 . - 0 gene_id "CNAG_04548"; transcr
1 ena exon 100 1497 . - . gene_id "CNAG_04548"; transcr
1 ena CDS 168 1497 . - 1 gene_id "CNAG_04548"; transcri
1 ena three_prime_utr 100 167 . - . gene_id "CNAG_04548
1 ena gene 5928 7982 . - . gene_id "CNAG_07303"; gene_s
1 ena transcript 5928 7982 . - . gene_id "CNAG_07303";
1 ena exon 7685 7982 . - . gene_id "CNAG_07303"; transc
1 ena CDS 7685 7769 . - 0 gene_id "CNAG_07303"; transcr
1 ena start_codon 7767 7769 . - 0 gene_id "CNAG_07303";
1 ena exon 5928 7626 . - . gene_id "CNAG_07303"; transc
1 ena CDS 6209 7626 . - 2 gene_id "CNAG_07303"; transcr
1 ena five_prime_utr 7770 7982 . - . gene_id "CNAG_0730
1 ena three_prime_utr 5928 6208 . - . gene_id "CNAG_073
1 ena transcript 6209 7769 . - . gene_id "CNAG_07303";
1 ena exon 7685 7769 . - . gene_id "CNAG_07303"; transc
1 ena CDS 7685 7769 . - 0 gene_id "CNAG_07303"; transcr
1 ena start_codon 7767 7769 . - 0 gene_id "CNAG_07303";
1 ena exon 6209 7626 . - . gene_id "CNAG_07303"; transc
1 ena CDS 6212 7626 . - 2 gene_id "CNAG_07303"; transcr
1 ena stop_codon 6209 6211 . - 0 gene_id "CNAG_07303";
1 ena gene 8766 9603 . - . gene_id "CNAG_07304"; gene_s
1 ena transcript 8766 9603 . - . gene_id "CNAG_07304";
1 ena exon 9311 9603 . - . gene_id "CNAG_07304"; transc
1 ena CDS 9311 9480 . - 0 gene_id "CNAG_07304"; transcr
1 ena start_codon 9478 9480 . - 0 gene_id "CNAG_07304";
1 ena exon 8766 9246 . - . gene_id "CNAG_07304"; transc
1 ena CDS 9129 9246 . - 1 gene_id "CNAG_07304"; transcr
Running htseq-count¶
So now we are ready! We run htseq-count using htseq-count ALIGNMENT_FILE GFF_FILE. Here is our command for our test sample:
–format=bam: Since Tophat generated BAM files for us
–stranded=reverse: The dUTP method that we used for generating a strand-specific library produces reads that are anti-sense, htseq-count considers this to be “reverse”.
[5]:
ls ${STAR_OUT}
21_2019_P_M1_S21_L001_R1_short_introns_Aligned.sortedByCoord.out.bam
21_2019_P_M1_S21_L001_R1_short_introns_Aligned.sortedByCoord.out.bam.bai
21_2019_P_M1_S21_L001_R1_short_introns_Log.final.out
21_2019_P_M1_S21_L001_R1_short_introns_Log.out
21_2019_P_M1_S21_L001_R1_short_introns_Log.progress.out
21_2019_P_M1_S21_L001_R1_short_introns_ReadsPerGene.out.tab
21_2019_P_M1_S21_L001_R1_short_introns_SJ.out.tab
21_2019_P_M1_S21_L002_R1_Aligned.out.bam
21_2019_P_M1_S21_L002_R1_Log.final.out
21_2019_P_M1_S21_L002_R1_Log.out
21_2019_P_M1_S21_L002_R1_Log.progress.out
21_2019_P_M1_S21_L002_R1_ReadsPerGene.out.tab
21_2019_P_M1_S21_L002_R1_short_introns_Aligned.sortedByCoord.out.bam
21_2019_P_M1_S21_L002_R1_short_introns_Aligned.sortedByCoord.out.bam.bai
21_2019_P_M1_S21_L002_R1_short_introns_Log.final.out
21_2019_P_M1_S21_L002_R1_short_introns_Log.out
21_2019_P_M1_S21_L002_R1_short_introns_Log.progress.out
21_2019_P_M1_S21_L002_R1_short_introns_ReadsPerGene.out.tab
21_2019_P_M1_S21_L002_R1_short_introns_SJ.out.tab
21_2019_P_M1_S21_L002_R1_SJ.out.tab
21_2019_P_M1_S21_L002_R1__STARtmp
genome_Log.out
multiqc_data
multiqc_report.html
[6]:
htseq-count --quiet \
--format=bam \
--stranded=reverse \
${STAR_OUT}/21_2019_P_M1_S21_L002_R1_Aligned.out.bam \
$GENOME_DIR/$GTF > ${COUNT_OUT}/21_2019_P_M1_S21_L002_R1.tsv
Let’s take a quick peek at the results
[7]:
head ${COUNT_OUT}/21_2019_P_M1_S21_L002_R1.tsv
CNAG_00001 0
CNAG_00002 51
CNAG_00003 23
CNAG_00004 78
CNAG_00005 5
CNAG_00006 546
CNAG_00007 188
CNAG_00008 119
CNAG_00009 29
CNAG_00010 146
There’s also some useful information at the end of the file:
[8]:
tail ${COUNT_OUT}/21_2019_P_M1_S21_L002_R1.tsv
ENSRNA049551862 0
ENSRNA049551899 0
ENSRNA049551942 0
ENSRNA049551964 0
ENSRNA049551993 0
__no_feature 18783
__ambiguous 316
__too_low_aQual 0
__not_aligned 46060
__alignment_not_unique 35006