SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Tophat uniquely mapped reads mrfox Bioinformatics 2 05-23-2013 05:58 AM
not uniquely mapped reads unidodo RNA Sequencing 2 04-22-2011 02:07 PM
Uniquely mapped reads with bowtie mapper Bioinformatics 2 11-22-2010 10:44 PM
BWA Uniquely Mapped Reads NF_seq Bioinformatics 0 09-06-2010 03:32 AM
cufflinks and non-uniquely mapped reads clariet Bioinformatics 1 05-08-2010 11:13 AM

Reply
 
Thread Tools
Old 12-23-2009, 10:33 AM   #1
subeet
Junior Member
 
Location: Philadelphia

Join Date: Dec 2009
Posts: 1
Default How to get uniquely mapped reads from Tophat

Hi,

I was wondering if anyone has used the "--max-multihits" option in Tophat. I specified "--max-multihits 1" in my command line with the hope of getting uniquely mapped reads; however, when I checked the run.log file, I noticed that the program still uses the default "--max-multihits 40". It seems that specifying the value of --max-multihits doesn't make a change in the analysis.

Has anyone experienced similar problems before?

Thanks,

Subeet
subeet is offline   Reply With Quote
Old 12-29-2009, 10:12 AM   #2
mgogol
Senior Member
 
Location: Kansas City

Join Date: Mar 2008
Posts: 193
Default

I've used -g 1 successfully, which is the same as --max-multihits according to the manual.
mgogol is offline   Reply With Quote
Old 08-28-2012, 02:00 AM   #3
wisense
Member
 
Location: Shanghai,China

Join Date: Sep 2011
Posts: 30
Default

If you want to get uniquely mapped reads from Tophat output bam file,you can use the command like this:
Code:
samtools view accepted_hits.bam | awk '$5==255{print $0}' > uniq_mapped_reads.sam
In the bam file ,when the 5 colum is 255,it means that this read is an uniquely mapped read,which also contain a string"NH:i:1" in the last field.You can see the statement of bam/sam format file for detail.
Also,you can use another command as follow,which will do the same as the first one:
Code:
samtools view accepted_hits.bam | grep -w "NH:i:1" > uniquely_mapped_reads.sam
wisense is offline   Reply With Quote
Old 08-29-2012, 12:02 AM   #4
simonandrews
Simon Andrews
 
Location: Babraham Inst, Cambridge, UK

Join Date: May 2009
Posts: 869
Default

We use -g 1 as part of our standard tophat pipeline. Seems to work OK and definitely doesn't produce the same results as running with the default options.
simonandrews is offline   Reply With Quote
Old 09-13-2012, 11:33 AM   #5
JQL
Member
 
Location: MO, USA

Join Date: Apr 2011
Posts: 82
Default

Quote:
Originally Posted by wisense View Post
If you want to get uniquely mapped reads from Tophat output bam file,you can use the command like this:
Code:
samtools view accepted_hits.bam | awk '$5==255{print $0}' > uniq_mapped_reads.sam
In the bam file ,when the 5 colum is 255,it means that this read is an uniquely mapped read,which also contain a string"NH:i:1" in the last field.You can see the statement of bam/sam format file for detail.
Also,you can use another command as follow,which will do the same as the first one:
Code:
samtools view accepted_hits.bam | grep -w "NH:i:1" > uniquely_mapped_reads.sam
Hi,

When I read the SAM file format, it states:
"MAPQ: ...A value 255 indicates that the mapping quality is not available."

Does that mean uniquely mapped?

thanks.
JQL is offline   Reply With Quote
Old 09-13-2012, 11:37 AM   #6
JQL
Member
 
Location: MO, USA

Join Date: Apr 2011
Posts: 82
Default

Quote:
Originally Posted by simonandrews View Post
We use -g 1 as part of our standard tophat pipeline. Seems to work OK and definitely doesn't produce the same results as running with the default options.
For RNA-seq study, is it too stringent to use -g 1? I recall I read somewhere, that one will miss gene famlily memebers.

what number is more optimal? the default is 20, i believe.
JQL is offline   Reply With Quote
Old 09-13-2012, 11:58 AM   #7
westerman
Rick Westerman
 
Location: Purdue University, Indiana, USA

Join Date: Jun 2008
Posts: 1,104
Default

Quote:
Originally Posted by JQL View Post
Hi,

When I read the SAM file format, it states:
"MAPQ: ...A value 255 indicates that the mapping quality is not available."

Does that mean uniquely mapped?

thanks.
For Tophat it does.

I am not sure who made the following observation on a Galaxy mailing list but it is what I refer to:

Quote:
Tophat/bowtie don't report mapping quality
values that are as meaningful as BWA, but there is some information in the
mapping quality values tophat reports. Tophat yields 4 distinct values for
its mapping quality values (you can do a "unique" count on the mapping
quality field of any SAM file from tophat to verify this):



255 = unique mapping

3 = maps to 2 locations in the target

2 = maps to 3 locations

1 = maps to 4-9 locations

0 = maps to 10 or more locations.
westerman is offline   Reply With Quote
Old 09-13-2012, 12:29 PM   #8
JQL
Member
 
Location: MO, USA

Join Date: Apr 2011
Posts: 82
Default

Quote:
Originally Posted by westerman View Post
For Tophat it does.

I am not sure who made the following observation on a Galaxy mailing list but it is what I refer to:
thanks. that info is helpful. It is makes senes to me now. Because it takes two mappings to calculate quality. So uniquely mapped reads do not have mapping quality to report, hence, 255.
JQL is offline   Reply With Quote
Old 11-27-2012, 05:44 PM   #9
Triple_W
Member
 
Location: China

Join Date: Jan 2012
Posts: 12
Smile

Quote:
Originally Posted by westerman View Post
For Tophat it does.

I am not sure who made the following observation on a Galaxy mailing list but it is what I refer to:
The information is helpful. Does tophat in Galaxy adopt such protocal?

However, I found .bam file from tophat.v2 used a different version:
255: unmapped reads
50: unique mapped reads
3: two mapped locations
1: 3 or 4 mapped locations
0: >4 mapped locations
Triple_W is offline   Reply With Quote
Old 11-27-2012, 09:33 PM   #10
wisense
Member
 
Location: Shanghai,China

Join Date: Sep 2011
Posts: 30
Default

Quote:
Originally Posted by Triple_W View Post
The information is helpful. Does tophat in Galaxy adopt such protocal?

However, I found .bam file from tophat.v2 used a different version:
255: unmapped reads
50: unique mapped reads
3: two mapped locations
1: 3 or 4 mapped locations
0: >4 mapped locations
You're right.I also find this change in tophat v2.0.4. In my opinion,maybe the best way to find the uniquely mapped reads is using the "NH:i:1" tag in the last field of every alignment.
wisense is offline   Reply With Quote
Old 11-28-2012, 05:56 AM   #11
westerman
Rick Westerman
 
Location: Purdue University, Indiana, USA

Join Date: Jun 2008
Posts: 1,104
Default

With version v2.0.6 I do not find this exact behavior. But perhaps it is because my reads get separated into 'accepted_hits.bam' and 'unmapped.bam'. The latter does, indeed, have all reads at mapq 255 but the former does not. The former has the mapq of 50, 3, 1 or 0. So when I look at the "interesting" reads -- i.e., those that map in accepted_hits.bam -- then just doing a samtools filtering by mapq is good enough.
westerman is offline   Reply With Quote
Reply

Tags
rna-seq, tophat

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 09:48 PM.


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