How to extract uniquely mapped paired end reads from bam file

serenaliao · 08-20-2014, 11:42 AM

Dear all,

I believe this questions has been partially brought up several times. But I cannot find a definite answer from all the previous discussion.

I am now working on bam files of paired-end RNA-Seq samples. And I want to extract uniquely mapped paired-end reads (no multiple alignments). Basically, from the summary of tophat below, I want to get the aligned pairs and remove the reads with multiple alignments and disconcordant alignments, which result in 43226526*(1-21.2%-12.6%)=28615960 pe reads.

Left reads:
Input: 160074322
Mapped: 63468641 (39.6% of input)
of these: 15015591 (23.7%) have multiple alignments (13075 have >20)
Right reads:
Input: 160074322
Mapped: 54678109 (34.2% of input)
of these: 12050250 (22.0%) have multiple alignments (7459 have >20)
36.9% overall read alignment rate.

Aligned pairs: 43226526
of these: 9158530 (21.2%) have multiple alignments
and: 5451746 (12.6%) are discordant alignments
23.6% concordant pair alignment rate.

I have tried several command in samtools, including
samtools view -c -f 1 -F 12 accepted_hits.bam
samtools view -c -hf 0x2 accepted_hits.bam.

But all failed to do the job. Could anyone give any suggestions?

Thank you!

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Uniquely mapped reads from BAM file (BWA)	a_mt	Bioinformatics	2	04-29-2014 11:46 PM
Extract paired end reads from sff file.	ojy	Bioinformatics	4	12-13-2012 05:07 AM
Extract mapped reads at certain position from paired end RNA-seq data	RedMary	Bioinformatics	2	09-04-2012 07:16 AM
How to extract uniquely mapped reads?	wisense	RNA Sequencing	0	05-16-2012 01:04 AM
Extract perfectly mapped reads from SAM/BAM file	Graham Etherington	Bioinformatics	2	07-21-2011 08:27 AM

08-20-2014, 11:42 AM	#1
serenaliao Member Location: Pittsburgh Join Date: Jan 2013 Posts: 22	How to extract uniquely mapped paired end reads from bam file Dear all, I believe this questions has been partially brought up several times. But I cannot find a definite answer from all the previous discussion. I am now working on bam files of paired-end RNA-Seq samples. And I want to extract uniquely mapped paired-end reads (no multiple alignments). Basically, from the summary of tophat below, I want to get the aligned pairs and remove the reads with multiple alignments and disconcordant alignments, which result in 43226526*(1-21.2%-12.6%)=28615960 pe reads. Left reads: Input: 160074322 Mapped: 63468641 (39.6% of input) of these: 15015591 (23.7%) have multiple alignments (13075 have >20) Right reads: Input: 160074322 Mapped: 54678109 (34.2% of input) of these: 12050250 (22.0%) have multiple alignments (7459 have >20) 36.9% overall read alignment rate. Aligned pairs: 43226526 of these: 9158530 (21.2%) have multiple alignments and: 5451746 (12.6%) are discordant alignments 23.6% concordant pair alignment rate. I have tried several command in samtools, including samtools view -c -f 1 -F 12 accepted_hits.bam samtools view -c -hf 0x2 accepted_hits.bam. But all failed to do the job. Could anyone give any suggestions? Thank you!