Question: How to get regions with mapped reads
0
gravatar for nanoide
2.1 years ago by
nanoide60
nanoide60 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 2.1 years ago by RamRS30k • written 2.1 years ago by nanoide60
2

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

ADD REPLYlink modified 2.1 years ago • written 2.1 years ago by genomax92k

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

ADD REPLYlink written 2.1 years ago by nanoide60
1

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

ADD REPLYlink written 2.1 years ago by RamRS30k

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

ADD REPLYlink written 2.1 years ago by nanoide60
Please log in to add an answer.

Help
Access

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