how to find "unmapped mates" from tophat output

[capricy]

Hello, there,

I looked at the tophat output and found that the "unmapped mate" of the mapped reads were not in the unmapped.bam file. I wonder if there is an easy way to find those reads?

Another question about the tophat:

If I just want to map the reads of the second end from paired-end sequencing results, how should I provide it to the tophat command?

For example, if I have paired end reads: read1, read2,

I will do

tophat <genome_index_base> <read1> <read2>

Then if I just want to map read2, may I use tophat <genome_index_base> <read2>

Thank you very much for the help

[WhatsOEver]

I looked at the tophat output and found that the "unmapped mate" of the mapped reads were not in the unmapped.bam file. I wonder if there is an easy way to find those reads?

Can you provide an example of a read in sam format that was not placed in the unmapped.bam though it was flagged as unmapped?

You could try to do it via samtools - see this post http://seqanswers.com/forums/showthread.php?t=12283

1) An unmapped read whose mate is mapped.
samtools view -u -f 4 -F264 alignments.bam > temp1.bam
2) A mapped read who's mate is unmapped
samtools view -u -f 8 -F 260 alignments.bam > temp2.bam
3) Both reads of the pair are unmapped
samtools view -u -f 12 -F 256 alignments.bam > temp3.bam

Then if I just want to map read2, may I use tophat <genome_index_base> <read2>

There is no PE/SE parameter in tophat. If you provide one fastq, it runs in PE mode, if you provide one in SE mode. So, yes you can

[capricy]

@WhatsOEver, thank you very much for reply.

Do you mean if I provide one fastq, tophat will run SE mode?

My understanding is, in PE mode, tophat somehow will get the reverse complementary sequences from read2. If I only provide read2, how does it know whether or not it should reverse-compliment them first?

Thanks

[dpryan]

It'll try to map to both strands regardless (at least unless you tell it not to).

[capricy]

If I have two fastq files, read1 and read2, which contain many unpaired reads after quality filtering, what is the best way to use tophat?


Still treat them as paired end read?
tophat <genome_index_base> <read1> <read2>

or tread them just as single end read

tophat <genome_index_base> <read1>, <read2>

Thanks for the help

[dpryan]

The best way is to use a trimmer that won't leave reads unpaired like that (e.g., trimmomatic or trim_galore). You would then map the pairs as pairs and the orphans as single-ended.

[capricy]

Why do I need to trim sequences? My issue is that some unpaired reads present in the two fastq file. For example, one read in read1, but its paired read is not in read2....

[dpryan]

Right and how did that happen? Likely because you or someone else trimmed or filtered them incorrectly. That's the root of your problem.

[WhatsOEver]

Why do I need to trim sequences? My issue is that some unpaired reads present in the two fastq file. For example, one read in read1, but its paired read is not in read2....

So you quality filtered the reads somehow (would be interesting to know how if it wasn't a trimmer) prior to tophat and you know that some reads were discarded because of this?! These reads will, of course, not be in the unmapped.bam afterwards...

Btw: The only thing considered in PE mode is that the reads have the opposite orientation. Read 1 can be on the "forward" strand as can be read 2. Unless a stranded-protocol was used.

[capricy]

Well, what if the unpaired reads are not resulting from the trimming step/quality screening, but a removal on purpose?

The sure thing is that read1 and read2 have the opposite orientation....

[dpryan]

Then the removal was done wrong. If one read in a pair is removed, then the other should be removed as well (and possibly written to a different file so it can be aligned single-ended). Otherwise you just end up desynchronizing the files.

[WhatsOEver]

@capricy: Do you actually still have a question or is it answered now?

If not, maybe its easier to express your question in an example:

1) Let's say you have 2x2 paired reads read1_1, read 1_2 and read2_1, read2_2.
2) Now, you remove read2_1 (for whatever reason).
3) The remaining 3 reads are mapped to your reference.
4) You would like to get read2_2 from the resulting sam file? Doesn't my previously posted command work for you?

2) A mapped read who's mate is unmapped
samtools view -u -f 8 -F 260 alignments.bam > temp2.bam

5) ???