Hello,

I have a genome with homologous pairs of chromosome. I am using hisat2 to map RNA-seq reads  to genome. I want to do variant calling with GATK, but before I want to extract the reads that mapped uniquely from BAM file(not SAM). Could you suggest how to do that? From other threads seems that flag NH:i:1 means mapped once, but how can I grep it from bam file? and is it a right way?

    HISEQ1:128:C0HUMACXX:6:2210:15441:67674	163	chr1A	1101543	60	59M	=	1104438	626	CAATCTTCTCGTGGTTTATTGGCGAATGGAGTTTCAGAAGATGCTCAAAATGGTTTCCA	@BCFFDDEFHDDFGEEEDHEHCBHFD<3D3C*?D:B?:B?B*98D?4?BB=CG)=FFHA	AS:i:0	ZS:i:0	XN:i:0	XM:i:0	XO:i:0	XG:i:0	NM:i:0	MD:Z:59	YS:i:0	YT:Z:CP	NH:i:1
    HISEQ1:130:C0G55ACXX:1:1203:2254:30173	163	chr1A	1101543	60	101M	=	1101599	155	CAATCTTCTCGTGGTTTATTGGCGAATGGAGTTTCAGAAGATGCTCAAAATGGTTTCCATGATGCTGATCGAACTGTCCGCAGAGGTGAGGGGCCATCTAA	CCCFFFFFHHHHHJGIJJJJJJJJJJJJJGIJIJJIJIJJJJJJJJJJJJJJJIIJJJJJJJJJJJJHHHHHFFFFEEEDDDDDDD@BDDDDDDDDDDDDD	AS:i:0	ZS:i:-15	XN:i:0	XM:i:0	XO:i:0	XG:i:0	NM:i:0	MD:Z:101	YS:i:-2	YT:Z:CP	NH:i:1

`samtools view -b -q 5 -o output.bam alignments.bam`

Having said that, there's no reason to bother doing this. GATK is smart enough to understand that low MAPQ alignments should be ignored.