Biostar Beta. Not for public use.
Question: Extracting uniquely mapped reads from bams
Entering edit mode

Hi everyone!

I generated multiple bam files with hisat2 using following commands:

for single-end reads:

hisat2 --max-intronlen 20000 -q -x hisat2_index/rheMac8_index -U myfile.fastq -S myfile.sam

for paired-end reads:

hisat2 --max-intronlen 15000 --score-min L,0,-0.1 --no-discordant -p 5 -q -x hisat2_index/rheMac8_index -1 myfile_1.fastq -2 myfile_2.fastq -S myfile.sam

Now, I would like to extract uniquely mapped (aligned exactly 1 time) reads from my bam files to be able to differentiate between reads that map to a gene and its pseudogene. I tried a few things with awk to select reads with mapq value of 60, but could not succeed.

I also changed the --score-min formula slightly to make the mapping more stringent and differentiate between reads that map to the original gene and to the pseudogene. I was wondering if anyone would have any tips about options that might help to differentiate between original vs. pseduogene mapped reads.

I would be more than happy if anyone could help me with these issues.

Best, Gökberk

ADD COMMENTlink 9 months ago gokberk • 30
Entering edit mode

I think this biostars thread might be helpful in your case.

ADD REPLYlink 9 months ago
Nitin Narwade
• 400
Entering edit mode

Also see How to select uniquely and concordantly reads from hisat2 alingment for raw read count and please use google and the search function. This has been adressed before.

ADD REPLYlink 9 months ago
Entering edit mode

Thanks for your responses Nitin and ATpoint. I've already seen the thread ATpoint mentioned, however the suggested code there gave me an empty bam file for some reason. That's why I opened a new issue. And on top of that, I was hoping to find some suggestions for mapping RNA-seq reads to a gene and its pseudogenes. Cheers.

edit: I've just tried the samtools view -bq 1 file.bam > unique.bam command in the thread Nitin has pointed. However, after I do that I still have some reads with mapq score other than 60 (as far as I know, uniquely mapped reads have a mapq score of 60).

ADD REPLYlink 9 months ago
• 30

Login before adding your answer.

Similar Posts
Loading Similar Posts
Powered by the version 2.0