Hi,

I'm trying to use htseq-count to map trimmed transcriptomic reads from a pure culture of a bacteria to its annotated genome (GTF/GFF file: image 1). When I use hisat2 to align paired-end reads to a database I've built using hisat2-build and the .fna file for the bacteria, I get 97% alignment (image 2). However, when I then run htseq-count on the (sorted) .bam file, 94% of the mapped reads correspond to "__no_feature" (image 3). The names of the reference chromosomes are the same in the .sorted.bam file as the GFF file ("DM14_1"; image 4). Additionally, the names of the identifiers (--idattr) are the same as in the GFF file ("FCFFPC_XXXXX", which does show up a little in the counts file).

Why might I be getting such a high percentage of reads aligning to __no_feature, given that I am mapping to a well-annotated genome with 97% alignment for the reads and I don't think I have an inconsistency in sequence names? Thank you very much in advance!

HISAT2 library building code (where the .fna file is formatted as multifasta with headers "DM14_1" and "DM14_2"):

hisat2-build bacterial_genome.fna genome_database

HISAT2 alignment code:

hisat2 -x genome_database --no-unal -1 forward.gz -2 reverse.gz -S alignment_file.sam

Sorting and indexing code:

samtools sort alignment_file.sam -o sorted_bam_file.bam
samtools index -M sorted_bam_file.bam

Counting code:

htseq-count --type CDS --idattr ID sorted_bam_file.bam gff_file.gff >  output_count_file.sam

Versions:

• hitsat2: 2.2.1
• samtools: 1.16.1
• htseq: 2.0.2

• Note: the GFF file is in format GFF3, and I have removed the FASTA sequences from the end of it because htseq-count would throw an error when it got to them. However, I wouldn't think this should cause an issue.