Tweet
Summary: How can I best get tophat to use only high-confidence junctions?
Details:
I'm working with mRNA in yeast, and I know there are about 350 junctions, which are pretty cleanly identified.
Now, I could come through and create pseudo-chromosomes and align to them with bowtie, since I know where the junctions are, but I'm trying to use tophat for this instead, particularly since I want to eventually extend my analysis to mammalian cells where things are going to get a lot less straightforward.
Now, my goal here is just to align my mRNA reads to the genome. I'm not actually trying to identify novel junctions.
I'm seeing three things in my output that surprise me, here.
1) I get ~ 14,000 junctions in my junctions.bed file.
2) Despite having a --min-anchor-length of 12, I see .bed entries like:
chrXIII 26606 33139 JUNC00000004 1 + 26606 33139 255,0,0 2 10,9 0,6524
My understanding of the meaning of blockSizes (the second to last field with "10,9", here) is that this is the anchor block sizes on the exons, and {10,9} < 12. I assume I'm just misunderstanding this, so I'd appreciate guidance on how I should be interpreting this.
3) I have not set --splice-mismatches, and it defaults to 0. Yet I see alignments in the output that look like:
SCS:1:100:1003:923#0 0 chrXII 457284 1 32M3775N8M * 0 0 AGGATTGGGTAATTTGCGCGCCTGCTGCCTTCCTTTCGTA B@CBCCC@@3@BCCC>BC=8A@@9@A;@3@CC@=@CBB>C NM:i:3 XS:A:- NH:i:3 CC:Z:= CP:i:457284
Note the NM:i:3, in particular. Doesn't this mean I'm getting an alignment with 3 mismatches in the exons, which I understand to be the "anchor" area in which --splice-mismatches is setting a max of 0 mismatches? What am I missing here?
***************
Ultimately, this is all just exposition. Really what I want is to make sure that I'm getting as much signal as possible in my read alignments as possible with very very little noise. My limited (!) understanding of my current output suggests to me that I've set things up to make the recall/precision tradeoff much more heavily in favor of recall, and I want to skew back towards the precision part of the curve.
How can I best achieve that?
Thanks!
-John
Details:
I'm working with mRNA in yeast, and I know there are about 350 junctions, which are pretty cleanly identified.
Now, I could come through and create pseudo-chromosomes and align to them with bowtie, since I know where the junctions are, but I'm trying to use tophat for this instead, particularly since I want to eventually extend my analysis to mammalian cells where things are going to get a lot less straightforward.
Now, my goal here is just to align my mRNA reads to the genome. I'm not actually trying to identify novel junctions.
I'm seeing three things in my output that surprise me, here.
1) I get ~ 14,000 junctions in my junctions.bed file.
2) Despite having a --min-anchor-length of 12, I see .bed entries like:
chrXIII 26606 33139 JUNC00000004 1 + 26606 33139 255,0,0 2 10,9 0,6524
My understanding of the meaning of blockSizes (the second to last field with "10,9", here) is that this is the anchor block sizes on the exons, and {10,9} < 12. I assume I'm just misunderstanding this, so I'd appreciate guidance on how I should be interpreting this.
3) I have not set --splice-mismatches, and it defaults to 0. Yet I see alignments in the output that look like:
SCS:1:100:1003:923#0 0 chrXII 457284 1 32M3775N8M * 0 0 AGGATTGGGTAATTTGCGCGCCTGCTGCCTTCCTTTCGTA B@CBCCC@@3@BCCC>BC=8A@@9@A;@3@CC@=@CBB>C NM:i:3 XS:A:- NH:i:3 CC:Z:= CP:i:457284
Note the NM:i:3, in particular. Doesn't this mean I'm getting an alignment with 3 mismatches in the exons, which I understand to be the "anchor" area in which --splice-mismatches is setting a max of 0 mismatches? What am I missing here?
***************
Ultimately, this is all just exposition. Really what I want is to make sure that I'm getting as much signal as possible in my read alignments as possible with very very little noise. My limited (!) understanding of my current output suggests to me that I've set things up to make the recall/precision tradeoff much more heavily in favor of recall, and I want to skew back towards the precision part of the curve.
How can I best achieve that?
Thanks!
-John