How to extract multi-mapped reads by samtools?

[mavishou]

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.

[westerman]

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.

[Marie_Noir]

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

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

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??!

[wzchang]

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

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??!

[ashmc]

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

[BnaT]

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).