How to select uniquely and concordantly reads from hisat2 alingment for raw read count

Question

[Bioinfonext]

I mapped pair end RNAseq read to genome using HISAT2 and got statistics like this:

31026735 (100.00%) were paired; of these:
    3579612 (11.54%) aligned concordantly 0 times
    25541051 (82.32%) aligned concordantly exactly 1 time
    1906072 (6.14%) aligned concordantly >1 times
    ----
    3579612 pairs aligned concordantly 0 times; of these:
      290787 (8.12%) aligned discordantly 1 time
    ----
    3288825 pairs aligned 0 times concordantly or discordantly; of these:
      6577650 mates make up the pairs; of these:
        6068493 (92.26%) aligned 0 times
        475539 (7.23%) aligned exactly 1 time
        33618 (0.51%) aligned >1 times
90.22% overall alignment rate

Now how I can select uniquely and concordantly reads for raw read count using samtools?

Answer 1

[Matteo Schiavinato]

Ok, let's go through it. The -f 0x2 bitwise flag selection keeps proper pairs, which are read pairs mapped in correct orientation (I suppose illumina reads facing each other) and within insert size (-X and -I) -> concordantly.

Now, if you want the uniquely mapped, you should basically exclude all those reads that have a secondary alignment. From the HISAT2 manual:

ZS:i:    <N>     Alignment score for the best-scoring alignment found other than the alignment reported. Can be negative. Only present if the SAM record is for an aligned read and more than one alignment was found for the read. Note that, when the read is part of a concordantly-aligned pair, this score could be greater than [AS:i].

Basically, if the line contains the ZS: field, it has a secondary alignment, i.e. it is not unique. Filter those out, for example using awk:

samtools view -h -f 0x2 file.bam | awk 'substr($1, 0, 1)=="@" || $0 !~ /ZS:/' | samtools view -h -b > filtered_file.bam

The awk part selects only lines that either start with @ (like the header) or don't contain the ZS: tag, like the reads with secondary alignments (i.e. non-uniquely mapping).

[Bioinfonext]

Thanks Macspider! I am also interested to use HTseq for read count which only count uniquely mapped read. Could anyone suggest how to proceed with hisat2 sam file to do read count using HTseq?

[Matteo Schiavinato]

Documentation is your friend! https://htseq.readthedocs.io/en/release_0.9.1/tour.html

[MCH]

This is very wrong. You did not understand the meaning of the option ZS. It only reports the alignment score of the secondary alignment which can be less than the primary one. As long as the primary one has a higher score, that alignment is considered uniquely mapped. The only way is to comare AS and ZS together; if they are equal it means that the alignment is actually a multimapper.

Answer 2

[al.bodrug]

I think like this you're gonna end up with unpaired reads where one of the mates has a ZS tag and the other doesn't and is the only one kept. So they won't be concordant.