SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
SFF Read names johan 454 Pyrosequencing 8 04-19-2012 08:54 AM
what're performaces of spliced read mappers? jiyan Bioinformatics 1 10-17-2011 08:06 AM
Mosaik - Read Names Trimmed sichan Bioinformatics 0 01-26-2011 01:20 PM
Tophat - Difference in mapping percentages of the read pairs adarshjose RNA Sequencing 0 12-17-2010 08:30 AM
SpliceMap 3.3.3.x released: Tool for spliced read alignment john_mu Bioinformatics 0 09-07-2010 03:34 PM

Reply
 
Thread Tools
Old 09-09-2010, 07:35 AM   #1
mendl7
Member
 
Location: USA

Join Date: Aug 2010
Posts: 10
Default Obtain Spliced Read Names/Mapping with TopHat??

Hi all. Well, I've finally made my first post.

Here's my question. I have run Tophat using 27M 76nt Illumina reads mapped against a small (120kb) genome looking for spliced reads. The run appears to complete successfully and I have the typical output files such as junctions.bed which has 582 junctions.

However, I would really like to look at some of the spliced read alignments. Even when I use the command line option to save temp files I can't seem to find the reads that TopHat believes are spliced. I can see some fastq files but, of course, the read names have been changed. Is there a way to look at finer scale results such as this? Thanks for any help you can give.
mendl7 is offline   Reply With Quote
Old 09-10-2010, 03:05 PM   #2
mendl7
Member
 
Location: USA

Join Date: Aug 2010
Posts: 10
Default

bumpbumpbump
mendl7 is offline   Reply With Quote
Old 09-10-2010, 05:47 PM   #3
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

Quote:
Originally Posted by mendl7 View Post
bumpbumpbump
Please do not bump, your question will be answered in time. Keep searching the forums.
nilshomer is offline   Reply With Quote
Old 09-13-2010, 05:20 PM   #4
malachig
Senior Member
 
Location: WashU

Join Date: Aug 2010
Posts: 117
Default

You should have both a SAM file and a wig file in addition to the junctions.bed you mention ... From the TopHat manual:

TopHat Output
1. accepted_hits.sam. A list of read alignments in SAM format. SAM is a compact short read alignment format that is increasingly being adopted. The formal specification is here.

2. coverage.wig. A UCSC BedGraph wigglegram track, showing the depth of coverage at each position, including the spliced read alignments.

3. junctions.bed. A UCSC BED track of junctions reported by TopHat. Each junction consists of two connected BED blocks, where each block is as long as the maximal overhang of any read spanning the junction. The score is the number of alignments spanning the junction.

Have you tried converting the SAM file to BAM format, indexing it, and then viewing in IGV?

TopHat manual
http://tophat.cbcb.umd.edu/manual.html

IGV
http://www.broadinstitute.org/igv/AlignmentData
malachig is offline   Reply With Quote
Old 09-15-2010, 09:23 AM   #5
mendl7
Member
 
Location: USA

Join Date: Aug 2010
Posts: 10
Default

The issue is that the accepted_hits.sam file do not contain reads that are spliced. I believe these are the reads that are simply mapped by Bowtie in the beginning of the TopHat run to determine regions of expression that TopHat tries to link later in the algorithm by spliced reads. I say this because there are no results in accepted_hits.sam that contain the XO: flag which gives the number of gap openings.

The coverage.wig file does contain results from the spliceds read but only as the depth of coverage or number of reads that mapped to the position.

If I've got something confused I would greatly appreciate the help.


What I am seeking is simply a list of reads that TopHat has called spliced and the reference positions that the read positions were mapped to.
mendl7 is offline   Reply With Quote
Old 09-15-2010, 10:35 AM   #6
Uwe Appelt
Member
 
Location: Heidelberg, Germany

Join Date: Oct 2009
Posts: 27
Default

@mendl7:

accepted_hits.sam indeed does contain split-read-alignments. At least for me, but i might miss your point. Does the following deliver want you're looking for?

Code:
fgrep --invert-match "50M" accepted_hits.sam > spliced_hits.sam
(you probably have to adjust the search string from 50M to whatever your read length was - you can check it by "head -n 100 accepted_hits.sam" [see column 6])

For my data this delivers for example:

Code:
HWUSI-EAS618_0001:3:63:1663:315#0       147     1       14788   255     42M140N8M       =       14783   0       CTCCGGGCCCCTCACCAGCCCCAGGTCCTTTCCCAGAGATGCCCTTGCGC      =1791>/79<==5'>>9<@/8>6AA=6=A<>@AA8AA?B@B?AB@BBA==        NM:i:0  XS:A:-  NH:i:1
HWUSI-EAS618_0001:3:11:1021:354#0       147     1       14791   255     39M140N11M      =       14720   0       CGGGCCCCTCACCAGCCCCAGGTCCTTTCCCAGAGATGCCCTTGCGCCTC      =>13;?=?>>::AA=A=A?;@?A>BBA@ABAB@AA@BBBABA@A@AB<@B        NM:i:0  XS:A:-  NH:i:1
HWUSI-EAS618_0001:3:34:352:749#0        147     1       14791   255     39M140N11M      =       14786   0       CGGGCCCCTCACCAGCCCCAGGTCCTTTCCCAGAGATGCCCTTGCGCCTC      A;>A4BAABB?B@ABB@BA@@A@AB@BAA?A@AB@BBBB@9BBBB@BBAB        NM:i:0  XS:A:-  NH:i:1
HWUSI-EAS618_0001:3:113:114:1346#0      163     1       14794   0       36M140N14M      =       14818   0       GTCCCTCACCAGCCCCAGGTCCTTTCCCAGAGATGCCCTTGCGCCTCATG      B=A@>BB@>B>5@B@A@?A?BAABBAABBA@ABBA===BAA@B@=B?A>6        NM:i:1  XS:A:-  NH:i:6  CC:Z:15 CP:i:102516182
HWUSI-EAS618_0001:3:91:1318:920#0       73      1       14796   0       34M140N16M      *       0       0       CCCTCACCAGCCCCAGGTCCTTTCCCAGAGATGCCCTTGCGCCTCATGAC      B@@BA@A=AAAAABB@B=BABBAA=A<BABABA<<=BA@=A7>A??BA:7        NM:i:0  XS:A:-  NH:i:6  CC:Z:15 CP:i:102516180
HWUSI-EAS618_0001:3:51:933:611#0        83      1       14809   255     21M140N29M      =       14743   0       CAGGTCCTTTCCCAGAGATGCCCTTGCGCCTCATGACCAGCTTGTTGAAG      CB@9@BBAA@BAC?@BBBBBBBBABBC=CABCCBBABBC@BCBBACBBCB        NM:i:0  XS:A:-  NH:i:1
HWUSI-EAS618_0001:3:117:1392:1226#0     147     1       14814   255     16M140N34M      =       14764   0       CCTTTCCCAGAGATGCCCTTGCGCCTCATGACCAGCTTGTTGAAGAGATC      7=<>?==49A9@AA?AAA@AAB?BBB?ABB;B@B;B?A?AB>BA@A=BAA        NM:i:0  XS:A:-  NH:i:1
HWUSI-EAS618_0001:3:113:114:1346#0      83      1       14818   0       12M140N38M      =       14794   0       TCCCAGAGATGCCCTTGCGCCTCATGACCAGCTTGTTGAAGAGATCCGAC      ;:)?B<A=>A6;-37A==>?A=;A@A<=??=/65A>A=BCAC<B>A>?BB        NM:i:0  XS:A:-  NH:i:6  CC:Z:15 CP:i:102516158
HWUSI-EAS618_0001:3:100:997:670#0       137     1       15000   0       39M757N11M      *       0       0       GATCCGACATCAAGTGCCCACCTTGGCTCGTGGCTCTCACTTGCTCCTGC      73A?A;=A@?AB=<9=BA<;;@AA?7@AB@BB9B@AA@?BBB@B@9B@@@        NM:i:0  XS:A:-  NH:i:8  CC:Z:12 CP:i:89817
HWUSI-EAS618_0001:3:101:1284:137#0      83      1       15000   255     39M757N11M      =       15101   0       GATCCGACATCAAGTCCCCACCTTGGCTCGTGGCTCTCACTTGCTCCTGC      ###########==55;=>744A107:=;?>AAA@B?>B@B8@BBABCAAB        NM:i:1  XS:A:-  NH:i:1
HWUSI-EAS618_0001:3:110:1622:914#0      83      1       15000   255     39M757N11M      =       14920   0       GATCCGACATCAAGTGCCCACCTTGGCTCGTGGCTCTCACTTGCTCCTGC      #:::==>?A;BAA=>=?==5?=?AA?BB@?BB=BA@@A?ACBBAA?B@BB        NM:i:0  XS:A:-  NH:i:1
HWUSI-EAS618_0001:3:29:668:1845#0       137     1       15000   0       39M757N11M      *       0       0       GATCCGACATCAAGTGCCCACCTTGGCTCGTGGCTCTCACTTGCTCCTGC      709:93>A?@?>=>>>A??;AA@?=?8?@ABAABB>AA>B?@@BBAABAA        NM:i:0  XS:A:-  NH:i:8  CC:Z:12 CP:i:89817
HWUSI-EAS618_0001:3:52:794:1776#0       83      1       15000   1       39M757N11M      =       14743   0       GATCCGACATCAAGTGCCCACCTTGGCTCGTGGCTCTCACTTGCTCCTGC      =;B<BBAABB@BBABBABA;BBBBBBBBBBBBBBBB@>;<ACCC>C;BCB        NM:i:0  XS:A:-  NH:i:4  CC:Z:12 CP:i:89817
HWUSI-EAS618_0001:3:65:849:538#0        73      1       15000   0       39M757N11M      *       0       0       GATCCGACATCAAGTCCCCACCTTGGCTCGTGGCTCTCACTTGCTCCTGC      ?@?7?@@AA+A@BBBB<AB,BB<B<?@A9BA<?BCCBBBC@BBCCBC6@B        NM:i:1  XS:A:-  NH:i:8  CC:Z:12 CP:i:89817
HWUSI-EAS618_0001:3:69:1365:376#0       83      1       15000   255     39M757N11M      =       15085   0       GATCCGACATCAAGTGCCCACCTTGGCTCGTGGCTCTCACTTGCTCCTGC      ###=:4:=;;=;44=56=;49=7;=;?;?=AA?AAA?A=AAB@BABBABB        NM:i:0  XS:A:-  NH:i:1
HWUSI-EAS618_0001:3:78:894:857#0        73      1       15000   0       39M757N11M      *       0       0       GATCCGACATCAAGTGCCCACCTTGGCTCGTGGCTCTCACTTGCTCCAGA      BCBBCBBB>BBBA?<BAABBABBB;@:><>=<?<<<B#############        NM:i:2  XS:A:-  NH:i:8  CC:Z:12 CP:i:89817
HWUSI-EAS618_0001:3:79:517:552#0        147     1       15000   255     39M757N11M      =       15096   0       GATCCGACATCAAGTGCCCACCTTGGCTCGTGGCTCTCACTTGCTCCTGC      57@6A<:>??BA=84;A=A6BBBBA=BBB@B@BBBCBB=CBC?BBBBBAB        NM:i:0  XS:A:-  NH:i:1
HWUSI-EAS618_0001:3:97:1036:226#0       73      1       15000   0       39M757N11M      *       0       0       GATCCGACATCAAGTGCCCACCTTGGCTCGTGGCTCTCACTTGCTCCTGC      A9@BBBBB@@BBBBBBBB;4;BB4B@B5B@BB<BCCBCACCCBC@CCCBB        NM:i:0  XS:A:-  NH:i:8  CC:Z:12 CP:i:89817
Cheers
Uwe

Last edited by Uwe Appelt; 09-15-2010 at 10:39 AM.
Uwe Appelt is offline   Reply With Quote
Old 09-18-2010, 11:22 AM   #7
mendl7
Member
 
Location: USA

Join Date: Aug 2010
Posts: 10
Default

yep, that's exactly what I needed. I had FALSELY assumed that XO flag would be used. Thank you very much for your assistance Uwe.
mendl7 is offline   Reply With Quote
Old 09-20-2010, 11:53 AM   #8
mendl7
Member
 
Location: USA

Join Date: Aug 2010
Posts: 10
Default

Just a quick update in case this is useful to anyone. I trim my reads to get rid of obviously lower quality 3' bases. Thus, fgrep --invert-match is not a useful filter since my read lengths are variable.
Instead I now use "awk '$6 ~/N/{print $_;}' accepted_hits.sam >spliced-hits.sam" which obviously looks for the presence of "N", inserted bases, in the CIGAR string column of the tophat SAM output.
mendl7 is offline   Reply With Quote
Reply

Tags
tophat rna-seq read names

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 02:16 AM.


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