SEQanswers

Go Back   SEQanswers > Applications Forums > RNA Sequencing



Similar Threads
Thread Thread Starter Forum Replies Last Post
Calculate number of multi-mapped reads? KAP Bioinformatics 13 02-17-2017 06:07 AM
How to extract mapped and unmapped raw reads from bwa's sam file ? vaibhavvsk Bioinformatics 11 02-07-2013 09:01 AM
Extract perfectly mapped reads from SAM/BAM file Graham Etherington Bioinformatics 2 07-21-2011 07:27 AM
samtools flagstats bug or it does not support multi-reads? xinwu Bioinformatics 0 12-22-2010 09:33 PM
extract subset (mapped reads) from csfasta and .qual files KevinLam SOLiD 1 01-18-2010 12:38 AM

Reply
 
Thread Tools
Old 06-10-2011, 06:53 AM   #1
mavishou
Junior Member
 
Location: China

Join Date: Jun 2011
Posts: 1
Default How to extract multi-mapped reads by samtools?

Hi guys,
Iím working with some RNA seq data and get some problems.

I mapped the paired-end reads on tophat and got the result in sam/bam format. Now Iím using samtools to analyze the result.

I want to extract the reads that map to more than one place in the genome, and this is my command line:
Samtools view Ėh Ėf 0x100 in.bam > out.sam
There are no output alignmens in the out.sam except the head, which means there are no multi-mapped reads

However, Iíve run my own program in perl and find that thereíre lots of reads whose IDs appear more than twice in the sam file, which means .these read mapped more than one place in the genome. (Two paired-end reads have the same ID, so most of the time, a certain ID appearing twice means itís a unique mapping)

I want to know whether my command line is right? I mean whether 0x100 is the right parameter to extract multi-mapped reads.
mavishou is offline   Reply With Quote
Old 06-10-2011, 12:24 PM   #2
westerman
Rick Westerman
 
Location: Purdue University, Indiana, USA

Join Date: Jun 2008
Posts: 1,104
Default

Quote:
Originally Posted by mavishou View Post
Hi guys,

I want to know whether my command line is right? I mean whether 0x100 is the right parameter to extract multi-mapped reads.
It could be. Depends on the program that created the BAM file. In other words perhaps the program thought that *all* matches were primary and thus never set the 0x100 ("the alignment is not primary") flag.

I suggest that, since you know the names of some of the reads that have multiple mappings, you look through the BAM file for those reads and see what flags are being set for those reads.
westerman is offline   Reply With Quote
Old 07-22-2011, 02:01 AM   #3
Marie_Noir
Junior Member
 
Location: Berlin

Join Date: Jul 2011
Posts: 8
Default

I have the same issue to work on. I know that there are multi-mapped reads in my alignment

Quote:
grep -i XA:Z:chr15 S1_unsorted.sam | wc -l
118242
but samtools view -f100 -c gives me 0

Checking the lines with XA:Tag

118242

Quote:
grep -i XA:Z:chr15 S1_unsorted.sam | head -n 10
61WG6AAXX_1:1:1205:11675 0 chr15 28901805 23 76M * 0 0 GCCTCCTGAGTGGCTGGGATTACGGTGTGTGCCACCACGCCTGGCTACTTTTGTCTTTTTAGTAGAGCTGG
GGGTT eeeceadd^deeceedd^T\\cc`ccddLTc\cTccTc\ceeec\baYYYbb^KKacYc[R^TTaL]L]Y```\Lc XT:A:U NM:i:3 X0:i:1 X1:i:1 XM:i:3 XO:i:0 XG:i:0 MD:Z:XA:Z:chr15,-20643137,76M,4
;
shows that actually the FLAG is not set. Question is, why not??!
Marie_Noir is offline   Reply With Quote
Old 06-06-2013, 07:17 AM   #4
wzchang
Junior Member
 
Location: Maryland

Join Date: Jun 2013
Posts: 1
Smile

You need give the right FLAG number 0x100 or 256, Both should work.

Quote:
Originally Posted by Marie_Noir View Post
I have the same issue to work on. I know that there are multi-mapped reads in my alignment



but samtools view -f100 -c gives me 0

Checking the lines with XA:Tag

118242



shows that actually the FLAG is not set. Question is, why not??!
wzchang is offline   Reply With Quote
Old 12-04-2016, 12:15 PM   #5
ashmc
Junior Member
 
Location: Rhode Island

Join Date: Nov 2016
Posts: 1
Default

Hello,

How would I get all the other locations in the genome where a particular read maps? If I use the 256 flag, it seems to only give a secondary alignment. I would like to rank the reads such that I can find the top 10 reads (top 10 that map to the most places in the genome).

Thanks
ashmc is offline   Reply With Quote
Old 12-05-2016, 05:27 AM   #6
BnaT
Member
 
Location: Barcelona

Join Date: Nov 2014
Posts: 10
Default

Quote:
Originally Posted by ashmc View Post
Hello,

How would I get all the other locations in the genome where a particular read maps? If I use the 256 flag, it seems to only give a secondary alignment. I would like to rank the reads such that I can find the top 10 reads (top 10 that map to the most places in the genome).

Thanks
If you only need to check a single read or a few, you could use BLAST. Otherwise, you could try bwa-mem with -a option to get the secondary alignments. Also you could try GEM-mapper which in default outputs all secondary alignments given a maximum threshold of mismatch (e.g 4% of the read length).
BnaT is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 12:16 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO