Biostar Beta. Not for public use.
Extracting uniquely mapped reads from bams
0
Entering edit mode
20 months ago
gokberk • 30
@gokberk54066

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

rna-seq hisat2 • 139 views
ADD COMMENTlink
2
Entering edit mode

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

ADD REPLYlink
2
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
0
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

Login before adding your answer.

Similar Posts
Loading Similar Posts
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.3