Question: How to get regions with mapped reads
gravatar for nanoide
4 months ago by
nanoide10 wrote:

Hello all,

So I would like to get the sequence of the genomic regions with mapped reads after an RNA-seq analysis.
First, I have mapped the reads into a reference genome using hisat2 and --rna-strandness F.
Does anyone have a suggestion on how I could extract the regions of the genome with mapped reads?
My plan was to use bedtools genomecov on the bam file to find the regions with coverage > 0 and then use bedtools getfasta using the bed with the regions and the genome in fasta format.
However, I'm not sure on how to keep strand information.

I think I have to include -strand in genomecov and -s in getfasta, but I have found that genomecov gives the same results with and without the -strand flag. I was expecting including -strand I would get strand information in the bed file but it is not the case. This is the command I have tried:

bedtools genomecov -ibam input.bam -strand -bga | awk '$4>0' >> output.bed

Does anyone have any comment?

Thank you for your time

ADD COMMENTlink modified 4 months ago by RamRS20k • written 4 months ago by nanoide10

A previous post that would be of interest: A: navigate IGV to the points of interest

ADD REPLYlink modified 4 months ago • written 4 months ago by genomax62k

Thank you very much, I think bamToBed worked for my purpose here

ADD REPLYlink written 4 months ago by nanoide10

Please use the formatting bar (especially the code option) to present your post better. I've done it for you this time.

ADD REPLYlink written 4 months ago by RamRS20k

Thank you for your help, I will do it next time

ADD REPLYlink written 4 months ago by nanoide10
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1397 users visited in the last hour