Low hits in express

[TLongi1]

Hi,

I'm working on RNA-seq data (Illumina Hi-Seq). After Trimmomatic, CLC (de novo assembly) and MIRA assembly and adding the Debris, I run Tophat to map my raw reads against the transcriptome.

Here's the code and the alignment summary:

tophat -p 64 -o Path/to/Output /path/to/bowtie/database /path/to/Paired_forward.fq.gz /path/to/Paired_reverse.fq.gz

The Output I get is:

Left reads
Input: 11157893
Mapped: 10404612 (93.2% of Input)
of These: 8239736 (79.2%) have multiple alignments (1167061 have >20)

Right reads
Input: 11157893
Mapped: 10404612 (93.2% of Input)
of These: 8239736 (79.2%) have multiple alignments (1167061 have >20)

93.2% Overall read mapping rate.

Aligned pairs: 10404612
of These: 8239769 (79.2%) have multiple alignments
10404497 (100.0%) are discordant alignments
0.0% concordant pair alignment rate.

This already does look not too good for me.
After sorting with samtools (samtools -n), express gives me about 10-30 hits per sample.

Does anybody have any idea what might be going wrong here?

Thanks a lot!

[GenoMax]

Having 80% of your data multi-map (when 93+% of the data is mapping) seems to indicate that there is some issue with either the data you started with or the assembly that was generated.

Was the quality of the sequences you started with good (reasonably clean FastQC results)? Is this data for an organism with no reference genome (or is a closely related genome available)? How did the assembly look (in terms of number of contigs, length of contigs, genes you expected to find)?

[Brian Bushnell]

Having 80% of your data multi-map (when 93+% of the data is mapping) seems to indicate that there is some issue with either the data you started with or the assembly that was generated.

I think it's OK in this case, as he's mapping RNA-seq data to a transcriptome. I'm more concerned with the fact that almost 100% of the alignments are discordant, which indicates the reads don't go together, or pair ordering was broken.

First, you need to re-pair the reads, or perhaps even validate that they are from the same library. You could alternatively reprocess them from the raw reads. Be sure to process read 1 and read 2 at the same time to keep them together. In fact, it might help if you displayed all the command lines you used up for processing up to mapping.

Secondly, for mapping RNA to a transcriptome, you should not use Tophat2, which is designed for spliced alignments to a genome; it's better to use Bowtie2 or something else.

EDIT: Looking at the statistics again, I think the mapping command was wrong, and the same file was specified twice. Note that the exact same number of read 1 was mapped as read 2, with the exact same number having multiple alignments, and the same number with mapq>20; that can't be a coincidence.

[TLongi1]

Brian: Oh yes, very good observation! I completely missed that. I will check it right away. Thanks

Max: The quality of the sequences was good. Yes, I sequenced the transcriptome of an ant species and there are no genomes of closely related species out there. That's why I had to do a denovo assembly. I agree that there seems to be a problem with the data and hope (also it would be a bit shameful..) that it is what Brian suggested

[TLongi1]

Indeed, I accidently used the Reverse reads twice rather than 1x reverse and 1x forward. This completely solved the problem. Thanks guys!

[GenoMax]

Ouch!

What did the new stats look like (compared to one's posted above)?