Htseq Count

Question

[Bjorn]

Hello,

I have ran htseq-count numerous times and continue to get the same error. That NONE of my genes are counted as seen here.

ZXDC    0
ZYG11B  0
ZYX 0
ZZEF1   0
ZZZ3    0
__no_feature    70257177
__ambiguous 0
__too_low_aQual 1509790
__not_aligned   3970775
__alignment_not_unique  4277765

However, I have a very high allignment rate when i use Hisat2

43268953 reads; of these:
  43268953 (100.00%) were paired; of these:
    5758279 (13.31%) aligned concordantly 0 times
    35403271 (81.82%) aligned concordantly exactly 1 time
    2107403 (4.87%) aligned concordantly >1 times
    ----
    5758279 pairs aligned concordantly 0 times; of these:
      151183 (2.63%) aligned discordantly 1 time
    ----
    5607096 pairs aligned 0 times concordantly or discordantly; of these:
      11214192 mates make up the pairs; of these:
        8563285 (76.36%) aligned 0 times
        2522477 (22.49%) aligned exactly 1 time
        128430 (1.15%) aligned >1 times
90.10% overall alignment rate

This is driving me crazy. What is wrong here?

hisat2code:

hisat2 -t --rna-strandness RF --summary-file /Volumes/cachannel/ZebraFinchBrain/CB-4a_genomemapping/summaryfilecb4a.txt -p 24 -x /Volumes/cachannel/ZebraFinchBrain/RefFolder/Reference -1 /Volumes/cachannel/ZebraFinchBrain/CB-4a_Fastqc_Report/CB-4a_R1_001_val_1.fq.gz -2 /Volumes/cachannel/ZebraFinchBrain/CB-4a_Fastqc_Report/CB-4a_R2_001_val_2.fq.gz | samtools view -@ 24 -Shu - | samtools sort -@ 24 -n -o /Volumes/cachannel/ZebraFinchBrain/CB-4a_genomemapping/alignmentcb4a.bam

Htseq code:

samtools sort /Volumes/cachannel/ZebraFinchBrain/CB-4a_genomemapping/alignmentcb4a.bam -o /Volumes/cachannel/ZebraFinchBrain/CB-4a_genomemapping/sorted_alignmentcb4a.bam

then need to index bam:

samtools index /Volumes/cachannel/ZebraFinchBrain/CB-4a_genomemapping/sorted_alignmentcb4a.bam

Then run:

htseq-count -r name -t gene -i gene -s yes -f bam /Volumes/cachannel/ZebraFinchBrain/CB-4a_genomemapping/sorted_alignmentcb4a.bam /Volumes/cachannel/ZebraFinchBrain/GCF_003957565.2/ncbi_dataset/data/GCF_003957565.2/genomic.gff > /Volumes/cachannel/ZebraFinchBrain/HTSEQ_withautomate/output_counts.txt

This is the only warning I get: Warning: Read J00138:157:HVCVCBBXX:4:1127:5051:41862 claims to have an aligned mate which could not be found in an adjacent line.

Answer 1

[GenoMax]

To sort alignments by name you need to use samtools sort -n .....
If you don't want to re-sort add option -r pos (instead of -r name) to htseq-count command since your file is sorted by position at the moment.

[Bjorn]

Thank you. This still gets me zero reads counted. I am getting an error that my reads have missing mates. Is there a setting to eliminate this issue?