Thank you very much. I will try this out.
• 0 views
•
link
written
by
KVC_bioinfo
extracting mapped reads from the bam file
Hello all,
I am trying to extract the mapped reads from a bam file for a specific region.
samtools view -h in.bam "chr1:regionstart-regionend" > out.bam
this gives me all the reads falling into that region. (which contains read that mapped into that region plus reads unmapped into that region)
I am looking for reads that mapped to that specific region. Which options can I include?
Or samtools in unable to do it?
Thank you in advance.
• 9,769 views
•
link
updated
by
Devon Ryan
•
written
by
KVC_bioinfo
•
written
by
KVC_bioinfo
2 answers
using http://lindenb.github.io/jvarkit/VcfFilterJdk.html after samtools to get the reads strictly in the regions:
samtools view -b -F 4 in.bam "chr1:regionstart-regionend" |\
java -jar dist/vcffilterjdk -e 'return record.getStart()>=regionstart && record.getEnd()<=regionend;'
• 0 views
•
link
written
by
Pierre Lindenbaum
This works well. Thank you very much.
• 0 views
•
link
written
by
KVC_bioinfo
samtools view -b -F 4 in.bam "chr1:regionstart-regionend" > out.bam
I assume by "unmapped in a region" you mean "having coordinates in a region but also being unmapped". So, just filter unmapped reads.
• 0 views
•
link
written
by
Devon Ryan
Suppose I have a long read which maps to few exons of a gene A and it does not map to rest of the exons of gene A. Now when I extract the reads that map to say exon 1, I want to have only reads that mapped to a region of exon one. When I use the command you gave, that extracts the reads falling in that regions. however, few reads in that region do not map to the region specified but they map somewhere else on the gene.
I could see that on IGV few genes were not mapping to the region I specified but still present in the bam file I got from above command.
• 0 views
•
link
written
by
KVC_bioinfo
few reads in that region do not map to the region specified but they map somewhere else on the gene.
it's not clear if you're talking about a read or about a (paired) read and its' mate.
• 0 views
•
link
written
by
Pierre Lindenbaum
my original reads are single end reads coming from nanopore sequencing.
• 0 views
•
link
written
by
KVC_bioinfo
Log in to answer this question.
More posts like this
-
samtools view not removing unwanted headers
written by bpf24 •
-
Separating mapped and unmapped reads from BAM file
written by marongiu.luigi
-
Extracting reads from given region
written by KVC_bioinfo
-
Alignment in IGV: extracting specific reads.
written by KVC_bioinfo
-
extracting reads from sam/bam file
written by KVC_bioinfo
-
extracting reads from bam file
written by KVC_bioinfo
-
Trying to create short test bam not working.
written by ariel.balter
-
Getting read counts per loci
written by JameB11 •
-
samtools view with specified regions, different output with and without -F 4
written by tonja.r
-
mapped/unmapped mate with samtools
written by Flavia •
How are you distinguishing between reads that "fall" in a position, from reads that "map" to it?
That is basically I am not able to do it. I am interested in extracting reads mapped in a specific region.
To most of us, there is no difference in what you've done and what you're asking for, so clarification is required. Please elaborate on what the difference is between what you've done and what you're asking for by editing your original post.
Thank you. I added more explanation.
please, stop asking questions and acknowledge the answers (green tick on the left) or comment the answers:
Hello Pierre,
All the posts you pulled here does not answer the question I asked in this post. I have always acknowledged the answer.
no: eg.: Understanding SAM file format answer is correct but not commented or not flagged as answered.
Got it. Changed it to answer. Thank you.
Perhaps you're unaware of how to acknowledge answers. There should be a check mark (in addition to up/down arrows) that you as the OP can toggle on/off.